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

Magnetic reconnection can power bright, rapid flares originating from the inner magnetosphere of accreting black holes. We conduct extremely high-resolution (5376 × 2304 × 2304 cells) general-relativistic magnetohydrodynamics simulations, capturing plasmoid-mediated reconnection in a 3D magnetically arrested disk for the first time. We show that an equatorial, plasmoid-unstable current sheet forms in a transient, nonaxisymmetric, low-density magnetosphere within the inner few Schwarzschild radii. Magnetic flux bundles escape from the event horizon through reconnection at the universal plasmoid-mediated rate in this current sheet. The reconnection feeds on the highly magnetized plasma in the jets and heats the plasma that ends up trapped in flux bundles to temperatures proportional to the jet’s magnetization. The escaped flux bundles can complete a full orbit as low-density hot spots, consistent with Sgr A* observations by the GRAVITY interferometer. Reconnection near the horizon produces sufficiently energetic plasma to explain flares from accreting black holes, such as the TeV emission observed from M87. The drop in the mass accretion rate during the flare and the resulting low-density magnetosphere make it easier for very-high-energy photons produced by reconnection-accelerated particles to escape. The extreme-resolution results in a converged plasmoid-mediated reconnection rate that directly determines the timescales and properties of the flare.


INTRODUCTION
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 (Albert et al. 2007;Aharonian et al. 2007Aharonian 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 41 erg/s (Abramowski et al. 2012), which is non-negligible compared to the total jet power of 10 42 −10 44 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, 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 quasi-periodic 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 plasmoidunstable current sheet of 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 electronpositron 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 (Sironi & Spitkovsky 2014;Guo et al. 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)1 .Resolving plasmoidmediated 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 time-scale is governed by the flux decay which is directly set by the reconnection rate (Bransgrove et al. 2021), which 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 highestresolution global 3D GRMHD simulations to-date to show that plasmoid-mediated magnetic reconnection in transient, non-axisymmetric 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 2 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.

NUMERICAL SETUP
Reconnecting current sheets are plasmoid-unstable for Lundquist numbers (Bhattacharjee et al. 2009;Uzdensky et al. 2010) 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 η num ∝ ∆x p , where p ≈ 2 depending on the details of the second order accurate algorithm.Thus, the constraint on S (Eq. 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 3 internal and 4 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 0.5 • − 7.5 • of each pole.In all of these runs we fix the radial domain to [1.2, 2000]r g and we use a minimum 10000r 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 φ ∝ max ρ/ρ max (r/r in ) 3 sin 3 θ exp (−r/400) − 0.2, 0 , normalized to the gas-to-magnetic-pressure ratio β = 2p/b 2 = 100.We replenish gas density ρ in low-density regions to maintain σ max = 25 where the magnetization σ = b 2 /(4πρc 2 ) is defined using the magnetic field strength b co-moving with the fluid, and fluid-frame restmass density ρ.We adopt an equation of state for a relativistic ideal gas with an adiabatic index of γ = 13/9, in between a fully relativistic gas γ = 4/3 and a fully non-relativistic gas γ = 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 2 with ion mass m i and Boltzmann's constant k B such that T > 1 indicates relativistic ion temperatures.

RECONNECTION-POWERED FLARES
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 2 , and rest-mass density ρ around the time of one such flares 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, non-axisymmetric (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 mass accretion rate.A macroscopic equatorial current sheet forms in the magnetosphere, extending from the horizon to the disk at x = r sin θ cos φ ≈ −5r g at z = r cos θ ≈ 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 σ max = 25 in our simulations, whereas in reality it is much more strongly magnetized (σ σ max ) pair plasma.Reconnection occurs locally in X-points where a field line breaks and reconnects to other field line (see insets in Figures 1D and 1E).In these X-points, reconnection heats the plasma up to T ∼ σ max = 25 (left panels) after which it is expelled from the layer at Lorentz factors up to Γ ∝ √ σ max = 5 (Lyubarsky 2005, see also Appendix B for an exploration of different σ max in 2D).
The flux is expelled through reconnection into the lowdensity region in between the large low-mode-number RTI modes accreting spirals.Electrons and positrons accelerated to non-thermal 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.
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 of 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 horizontal magnetic field into 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 time scale 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 mix 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 anti-parallel field lines (top and bottom rows in Figures 1 and 2).
Figure 3A 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 (3B) there are no structures at T > 1.During flare, the (green) field lines in the current sheet (i.e., seeded in the T > 1 region in 3A) 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 3B) no such distinction is visible and all field lines (green and blue, which are seeded at the same points as in panel 3A) are part of the disk.The extreme resolution allows to capture multiple plasmoids identified as 3D helical field line structures in the sheet (Figure 3C) 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 1D.
Figure 4A-D shows 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 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, 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  −gdθdφ through the inner 5r g (4G), where g is the metric determinant, u µ is the fluid 4velocity, * 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 Figure 4E) 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 mass accretion rate (Figure 4G) 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.
For the high resolution run (red line, Figure 4E), 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 line) 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.Mini-flares (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 mini-flares are typically not accompanied by a clear drop in ṁ5rg (Figure 4G), while large flares are showing a clear drop in ṁ5rg 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 (Figure 4A,B).Mini-flares 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 ex- pressed 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 4F.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 4H we then measure the flow speeds from left and right of the current sheet as O(10 5 ) with explicit resistivity η = 5 × 10 −5 .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 two to three 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 time scale.
In Figure 5 we show 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 mini-flare 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 mini-flares could potentially result in smaller very high energy flares and shorter variability time scales (Abramowski et al. 2012).

RADIATIVE PROPERTIES OF THE RECONNECTION LAYER
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. 2020Crinquand et al. , 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.

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 σ 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 σ up = B2 up /(4πnm e c 2 ) = 2U B /(nm e c 2 ), where n is the number density of electrons with mass m e , B up is magnetic field strength upstream from the current sheet and magnetic energy density U B = B 2 up /8π 2 .The plasma particles are impulsively accelerated by non-ideal electric fields at the X-points (Sironi & Spitkovsky 2014).When they encounter plasmoids, they experience strong synchrotron losses.To parametrize the effect of the radiation backreaction, we define the particle Lorentz factor γ sync rad 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): where σ T = (8π/3)r 2 e is the Thomson cross-section, r e = e 2 /(m e c 2 ) is the classical electron radius, and e is the electron's charge.We then find (γ sync rad ) 2 = 3v rec B cl /(2cB up ), where B cl = m 2 e c 4 /e 3 6 × 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 − 150G 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 42 − 10 44 erg/s (Prieto et al. 2016) to the Blandford-Znajek jet power L BZ = κΩ 2 BH φ2 BH /(4π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 9 M for M87 * , and ) is the horizon radius (Blandford & Znajek 1977;Tchekhovskoy et al. 2011), resulting in a range B horizon ∼ 20 − 200 G at the horizon.By normalizing to a fiducial B up = 100 G in this range, we then obtain (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 Γ ∼ √ σ 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 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 chargestarved gap in the funnel region λ 10 3 (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 14 (M/6 • 10 9 M )(B up /100G), such that σ up ∼ 10 14 (M/6 • 10 9 M )(B up /100G)/2λ.For these parameters, γ sync rad σ up such 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 2 ) (Beloborodov 2017).
Using Eq. 2 and the ω B /Ω BH relation, we can rewrite B = v rec wω B /(c 2 (γ sync rad ) 2 ) and obtain (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 41 erg/s low-energy photons with energy density U soft rad ∼ 0.003 erg cm −3 in the inner 10r g results in γ IC rad ∼ γ sync rad U B /U soft rad ∼ 10 9 (Broderick & Tchekhovskoy 2015; EHT MWL Science Working Group et al. 2021), which is well above γ sync rad .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 rec rad ∼ 0.1U 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 burnoff limit E ph ∼ (γ sync rad ) 2 ω B ∼ 200MeV (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(E ph ) = m e c 2 γ sync rad ∼ 0.511MeV•γ sync rad ∼ few TeV.Additionally, particles can be accelerated beyond γ > γ sync rad 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 2 = 0.5MeV, e ± pairs are created if there are enough photon-photon collisions with seed photons with low energy E s ∼ (m e c 2 ) 2 /E ph .Highenergy photons of energy E ph,TeV produced in the magnetospheric region around the current sheet will interact most efficiently with seed photons of energy E s ∼ (1TeV/E ph,TeV ) eV.Given the uncertainties about the density of a 1eV 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 41 − 10 43 , is emitted in very high-energy photons, a γ-ray flux of 10 39 − 10 41 erg/s can be emitted as a flare.
Our extreme resolution GRMHD simulation shows transient flaring periods where the mass accretion rate drops (and, thus, luminosity of seed photons) signifi-cantly, 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 of 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 230GHz images (Chatterjee et al., in prep.).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 timescale (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 3 , 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, that can emit synchrotron photons with energies ranging from the radio to optical wavelengths (see Appendix A).

Sgr A * flares powered by radiative reconnection
Sgr A * shows daily near-infrared and X-ray flares from the inner 10r g , on average every 6 and 12 hours, lasting for 30-80 minutes, respectively (Baganoff et al. 2001;Eckart et al. 2006;Gravity Collaboration et al. 2018;Witzel et al. 2021;Murchikova & Witzel 2021).The flare periods in our simulation last for ∼ 100r g /c ∼ 30 minutes, and the subsequent quiescent period for ∼ 2000r g /c ∼ 10 hours for Sgr A * .The resulting hot spot orbits for ∼ 500r g /c ∼ 150 minutes in the inner 20r g until it diffuses due to mixing instabilities.The magnetic field strength in quiescence is well constrained in the range of 10−50 G in the inner 10r g for Sgr A * with black hole mass 4•10 6 M (Dodds-Eden et al. 2009).Using Eq. 3, this results in γ sync rad ≈ 9 • 10 6 (B up /10G) −1/2 , limiting the energy of accelerated particles by syn-chrotron cooling for a typical magnetization σ up ∼ 10 10 (M/4 • 10 6 M )(B up /10G)/2λ γ sync rad .Using Eq. 4, the compactness is B ∼ 10 −5 (w/1r g )(M/4 • 10 6 M )(B up /10G) 2 .Synchrotron photons emitted by the particles accelerated to the highest energies in the reconnection layer, up to γ sync rad ∼ 10 7 , should extend in the hard X-ray range.The energy of particles accumulated in the orbiting hot spot will be constrained by the synchrotron cooling time which has to be larger than the lightcrossing time of the current sheet, ct sync /w ∼ 1/( B γ) ≥ 1, or γ 1/ B ∼ γ cool = 10 4 for the hot spot at ∼ 10r g .These particles are likely to emit in the (near-)infrared range, (γ cool ) 2 ω B ∼ 1eV(B/10G).Thus, reconnection near the event horizon can power a multi-wavelength flare solely by synchrotron emission from reconnection-accelerated particles.Mini-flares are a potentially viable route to produce only near-infrared emission without strong enough X-rays to be detected as flares, as they do not produce a long-lasting extended current sheet, which would be the source of highest energy particles.Mini-flares could be distinguishable from large flares with concurrent multi-wavelength observations of Sgr A * (e.g., Ponti et al. 2017).The characteristic power of the X-ray emission can be estimated from the total dissipated power in reconnection, ∼ 0.1L BZ ∼ 10 35 (B horizon /10G) 2 erg/s.Thus, reconnection in the magnetospheric current sheet provides enough energy to power the observed X-ray flares from Sgr A * with typical luminosities in a range 10 34 − 10 35 erg/s (Neilsen et al. 2015).

CONCLUSIONS
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 non-axisymmetric 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 − 10r 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.Reconnectionheated 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 time scales 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 -Xpoints 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 sense to their black hole; Porth et al. 2021).However, due to limited numerical resolution they did do not capture plasmoidmediated 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 reconnectionheated 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, that 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;Chatterjee et al. 2021;Dexter et al. 2020;Yoon et al. 2020, and references therein) or two-temperature GRMHD approaches (Ressler et al. 2015;Chael et al. 2019) therefore cannot describe the non-thermal emission from these events which involves reconnection in high-σ collisionless pair plasma regime, the transport and cooling of non-thermal 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 an 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 non-ideal effects beyond scalar resistivity (e.g., Ripperda et al. 2019b) into GRMHD simulations, such as electron inertia and anisotropic electron pressure tensor effects in the 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. 2021Crinquand et al. , 2020) ) are crucial for probing the non-thermal 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 (Chatterjee et al., in prep.).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 non-thermal 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 with concurrent Event Horizon Telescope and multi-wavelength 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(Zhdankin et al. , 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 6 (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 non-ideal 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 plasmoidunstable current sheet close 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  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 E ph is τ γγ ∼ 0.1σ T U s rad w/E s , where U s rad is the energy density of photons at energy E s ∼ (m e c 2 ) 2 /E ph such that the photon-photon pair production opacity is maximal for E s E ph ∼ (m e c 2 ) 2 .Most of the reconnected power, L rad ∼ 0.1L jet , is radiated around the burnoff limit, ∼ 200MeV, and the peak can be quite broad (Hakobyan et al. 2019).Estimating U s rad ∼ L rad /4πw 2 c, we get τ γγ ∼ 0.1σ T L rad /(4πcwE s ).For w ∼ r g , and for typical photon energies E ph ∼ E s ∼ MeV, we find τ γγ ∼ 10 −2 σ T L jet /(4πcr g E s ) ∼ 10 −4 for a jet power L jet ∼ 10 43 erg/s.The rate of pair creation is then Ṅpair ∼ τ γγ • (0.1L jet /E ph ) ∼ 10 44 s −1 .We can compare this estimate to the Goldreich Julian number flux ṄGJ = I GJ /e, where I GJ (cL jet ) 1/2 is the Goldreich Julian current, such that ṄGJ ∼ (cL jet ) 1/2 /e ∼ 10 36 s −1 .The resulting multiplicity λ = Ṅpair / ṄGJ ∼ 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 γ 2 ∼ 10 4 GHz(γ/400) 2 , ranging from the radio to optical wavelengths.

APPENDIX B. INFLUENCE OF MASS LOADING ON PLASMA HEATING DUE TO RECONNECTION
We performed two additional 2D GRMHD simulations to show that reconnection heats the plasma to T ∼ σ max and Γ ∼ √ σ max , for the floor magnetizations in the jet σ max = 25, 60, 100. Figure 6 shows both the Lorentz factor γ (top row) and temperature T (bottom row) for the three values of σ max .The plasma in the current sheath is indeed heated to T ≈ σ max and Γ ≈ √ σ max and the reconnection exhaust deposits hot plasma in the jet sheath up to large distances.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 During mini-flares, 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 mini-flare 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.

Figure 1 .
Figure 1.Plasmoid-mediated reconnection, which takes place at sufficiently high resolutions in MHD, is seen in a 3D GRMHD simulation for the first time.Resolving the dynamics of X-points and plasmoids in the current sheet can be the key to understanding the source of black hole non-thermal emission, e.g., high-energy flares.Dimensionless temperature T = p/ρ, plasma-β, and density ρ (from left to right) in the meridional plane before (top row), during (middle row) in the inner 10rg and after (bottom row) the large magnetic flux eruption in the inner 40rg.During the magnetic flux eruption, the accretion disk is ejected and the broad accretion inflow is reduced to a thin plasmoid-unstable current sheet, indicated by X-points and magnetic nulls shown by the antiparallel in-plane field lines (in green, see inset in panel D) and the high β (inset panel E).The hot (T ∼ σmax) exhaust of the reconnection layer heats the jet sheath.Reconnection transforms the horizontal field in the current sheet to vertical field that is ejected in the form of hot coherent flux tubes (panel G) at low β and density (panels H,I).

Figure 2 .
Figure2.Our extreme resolution simulation reveals small-scale structure and interface instabilities of magnetic flux bundles escaping from the black hole, in an equatorial slice through the system.Dimensionless temperature T = p/ρ, plasma-β, and density ρ (from left to right) in the equatorial plane before a large magnetic flux eruption (top row), during the magnetic flux eruption (middle row) in the inner 10rg and after the magnetic flux eruption (bottom row) in the inner 40rg.Gaps of low β and density form during the pre-eruption quiescence while many azimuthal RTI modes accrete.During the magnetic flux eruption a single large T > 1 spiral forms with a gap where the sheet moved out of the equatorial plane.Magnetic flux escapes through the spiral current sheet, while accretion continues over a small angle φ < 2π at x ≈ 2rg and y ≈ −1 to y ≈ −2.In the bottom row the inner 10rg is in quiescent accretion state, and a hot flux tube that is ejected from the reconnection layer is in orbit at x ≈ 10rg to x ≈ 30rg and y ≈ −10rg to y ≈ 20rg.The low β flux tube shows clear signatures of instabilities at its boundaries mixing low density plasma into the disk.

Figure 3 .
Figure 3. Volume rendering of the temperature T = p/ρ shows plasmoids and hot current sheets.Extreme resolution allows the current sheets to become thinner and hotter than typically seen in GRMHD simulations.(Panel A:) During a large flare a relativistically hot T > 1 spiral current sheet forms.Accretion occurs over a small azimuthal angle φ < 2π in the T < 1 (white) regions.The green field lines, seeded in the current sheet (T > 1), remain in the current sheet and are mostly attached to the black hole.Blue field lines are seeded in the disk, where some disk field lines are accreting onto the black hole in the T < 1 region.(Panel B:) In the quiescent state T ≤ 1 everywhere, and both green and blue field lines (with the same seeds as in panel A) are in the disk, accreting onto the black hole.The inset (C) shows a zoom into the inner rg in the flare state with multiple escaping flux loops (green field lines).In the small black box we highlight an escaping flux tube with vertical field as the result of reconnection (green) and an infalling flux tube (purple).We also show a plasmoid, indicated by the helical field line (green) in the second small black box.

Figure 4 .
Figure 4.The equatorial current sheet that forms during the magnetic flux eruption is unresolved at low and standard resolutions (panels A,B) such that magnetic field lines (green lines) diffuse through the current sheet and do not reconnect, due to the high numerical resistivity.At high and extreme resolutions (C,D), the field lines are antiparallel in the current sheet, and they reconnect in well-defined X-points.Smaller current sheets are resolved in the accretion disk at high and extreme resolutions, potentially heating the plasma through reconnection.Panel E shows the magnetic flux on the horizon for the four numerical resolutions.The extreme and high resolution runs show two and three large flare periods, respectively, indicated by flux decay at a rate ∝ e −t/500 governed by the reconnection rate (dashed black lines).A mini-flare is indicated by the small flux drop at t ≈ 6800rg/c in the extreme resolution run.The standard and low resolution runs show a faster flux decay ∝ e −t/350 governed by the enhanced reconnection rate due to an increased numerical resistivity.Flares in the extreme resolution run are accompanied by clear drops in the mass accretion rate (panel G), due to the expulsion of the disk over a large azimuthal angle.Panel F shows a cut through the equatorial current sheet at x ≈ 1.5rg during the flare state (indicated by the red dashed line in panels A-D).Both the (nearly) radial field B L and the (nearly) toroidal field B M components (see definition in the text) change sign in the equatorial current sheet, while B N is (close to) zero, indicating zero-guide field reconnection.Panel H shows the flow speeds left and right of the current sheet.After correcting for the bulk flow, we measure the reconnection rate to be vrec ≈ 0.01c.We confirm this measurement at 10 radial cuts during separate flare periods.
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 (Figure4F).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 A w/η num = (v rec /c) −2 = 10 4 .Reconnection thus occurs in the asymptotic plasmoid-mediated regime where S ≥ S crit = 10 4(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 A = σ up /(σ up + 1)c ∼ c, for σ up = 25, and numerical resistivity η num .The reconnection rate is consistent with 2D resistive GRMHD simulations of plasmoiddominated reconnection in MAD flows (RBP20) for Lundquist numbers S = L sheet c/η

Figure 5 .
Figure 5. Smaller flux eruptions show shorter current sheets, potentially powering mini-flares that are not accompanied by a large-scale evacuation of the accretion disk.Meridional (top row) and equatorial (bottom row) cuts of temperature T = p/ρ (left), plasma-β (middle) and density ρ (right) during the mini-flare at t = 6852rg/c.The magnetic flux is expelled through a smaller (w 3rg) current sheet, close to the horizon, in a short time 100rg/c.The accretion disk is not expelled over a large azimuthal angle, yet the flare is accompanied by a significant drop in mass accretion rate (see Figure 4G) and clear gaps in the density (F).Multiple small current sheets are visible in the accretion disk at x ≥ 3rg indicated by the high plasma-β (B).

Figure 6 .
Figure 6.Lorentz factor Γ (top rows) and temperature T = p/ρ for the extreme resolution 3D run with density floor σmax = 25 (left), and two supplementary 2D runs with density floors σmax = 60 (middle) and σmax = 100 (right).Reconnection heats up the plasma in the equatorial current sheet to T ∼ σmax and Γ ∼ √ σmax.The hot reconnection exhaust heats up the jet sheath up to temperatures proportional to the magnetization in the jet.We observe limb-brightened jets up to at least 20rg.

Figure 7 .Figure 8 .
Figure 7. Top row: Profiles of the magnetic field in minimum variance coordinates in a current sheet in the high (A), standard (B), and low (C) resolution runs, as in Figure 4. Bottom row: Profile of the inflow speed into the current sheet, showing a reconnection rate of 0.01c for high resolution (D), and enhanced reconnection rates of 0.02c (standard resolution, E) and 0.025c (low resolution, F) as a result of a larger numerical resistivity and a broader current sheet (top panels).

Figure 9 .
Figure 9. Dimensionless temperature T = p/ρ (left), plasma-β (middle), and density ρ (right) in the meridional plane (first and third rows) and equatorial plane (second and fourth rows) during a magnetic flux eruption in the inner 10rg (first and second row) and during quiescence in the inner 40rg (third and fourth row).
and flaring hot spots from Sgr A * .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 Ac-cretion (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 Dutch Research Council (NWO) VICI award, grant Nr. *