© The Authors 2026

1. Introduction

Accreting black holes are frequently associated with the launching of powerful relativistic jets, originating from the innermost regions of the accretion flow and possibly powered by the black hole ergosphere (see Blandford et al. 2019 for a review and references therein). In high-luminosity systems, the accretion rate is large enough for the flow to cool efficiently, resulting in a geometrically thin, optically thick Keplerian disk (Shakura & Sunyaev 1973). In contrast, at low accretion rates, radiative cooling becomes inefficient, and the accretion flow becomes hot, geometrically thick, and more tenuous (Narayan & Yi 1994). In fact, the plasma density is sufficiently low that the flow becomes effectively collisionless, with Coulomb collision times exceeding the accretion timescale (e.g., Quataert et al. 2002).

This regime is particularly relevant for the nearby supermassive black holes M87^★ and Sgr A^★, where recent horizon-scale observations provide some of the most stringent empirical constraints on accretion physics. A key conclusion inferred from polarimetric measurements of these two sources is the presence of an organized vertical magnetic field on large scales (Johnson et al. 2015; GRAVITY Collaboration 2018, 2020; Event Horizon Telescope Collaboration 2021, 2024; Wielgus et al. 2022). Such a dynamically important magnetic field can dramatically alter the nature of accretion by exerting pressure and stress on the flow. Although Sgr A^★ does not exhibit a detectable jet, it nevertheless displays puzzling activity characterized by frequent nonthermal flares originating from the immediate vicinity of the black hole horizon (Baganoff et al. 2001; Genzel et al. 2003), which indicates ongoing particle acceleration (Markoff et al. 2001; Dodds-Eden et al. 2010).

In Sgr A^★, the accretion flow most likely results from the collective contributions of massive star winds orbiting the black hole (Paumard et al. 2006). Dedicated studies show that the plasma penetrates the black hole vicinity with a low net angular momentum (Cuadra et al. 2005; Ressler et al. 2018), which in turn suggests that accretion may be close to spherical. In the presence of a large-scale vertical magnetic field, the classical Bondi accretion model is significantly altered (Bisnovatyi-Kogan & Ruzmaikin 1974; Igumenshchev & Narayan 2002; Ressler et al. 2021). The frozen-in magnetic field lines are dragged inward with the plasma as it is gravitationally pulled towards the black hole event horizon. As accretion proceeds, field lines become increasingly stretched in the radial direction, leading to a partial storage of the gravitational potential energy in the field.

This process continues until the field becomes dynamically important near the hole. At this stage, the restoring force from the magnetic tension pushes the accretion flow away from the horizon. This phenomenon is accompanied by a large-scale reconnection event that facilitates the global reorganization of the field and yields efficient non-thermal particle acceleration (Bisnovatyi-Kogan & Ruzmaikin 1974, 1976; Meszaros 1975). This is reminiscent of flux eruptions reported in general relativistic magnetohydrodynamic (GRMHD) simulations of magnetically arrested disks (e.g., Tchekhovskoy et al. 2011; Dexter et al. 2020; Porth et al. 2021; Ripperda et al. 2022). This scenario provides a promising explanation for the recurrent flares observed in Sgr A^★ mentioned earlier. Nevertheless, this picture remains incomplete because the questions of dissipation and particle acceleration have not yet been fully settled. Fluid models have been widely used to elaborate the above picture, but they are unable to capture the microphysics of collisionless plasmas and particle acceleration. A kinetic plasma approach must therefore be employed.

Recent global general relativistic particle-in-cell (GRPIC) simulations of magnetized spherical accretion around a maximally spinning black hole have uncovered a clear connection between flux eruptions and particle acceleration (Galishnikova et al. 2023; Vos et al. 2025). Building on these previous works, we present here a new series of two-dimensional (2D) GRPIC simulations of collisionless spherical accretion for a Schwarzschild black hole immersed in an initially vertical laminar field (Sect. 2). This configuration, by design, prevents the formation of a spin-powered jet (Blandford & Znajek 1977), in contrast to previous GRPIC studies, and thus focuses on the accretion process. Solutions are integrated over an unprecedentedly long timescale so that a quasi-steady state can be reached and multiple flux eruptions can be observed in the presence of a pair plasma (Sect. 3) or an electron-ion plasma (Sect. 4). We further dissect each stage of the accretion process and propose a simplified model to construct a quantitative physical framework based on the simulation results (Sect. 5). Results are then discussed in the context of Sgr A^★ flares and isolated stellar-mass black holes accreting from the interstellar medium (Sect. 6).

2. Methods

In this work, we used the GRPIC code Zeltron (Parfrey et al. 2019), where the fields and particle properties are derived in the 3+1 formalism of Komissarov (2004) using spherical Kerr-Schild coordinates. In this formalism, B and D are the magnetic and electric fields as measured by fiducial observers (FIDOs) and H and E are auxiliary, metric-induced magnetic and electric fields. We modeled zero-angular-momentum accretion of matter by a Schwarzschild black hole of radius r_H = 2r_g, where r_g = 𝒢M_BH/c² is the gravitational radius and M_BH is the black hole mass. In this article, we employed Gaussian-cgs units for the electromagnetic fields and used r_g and t_g = r_g/c as explicit units of space and time, but we otherwise set 𝒢 = c = M_BH = 1.

We employed a 2D spherical grid of size N_r × N_θ. The grid is logarithmically spaced in the radial direction and uniform in θ, and it covers the domain [r_min, r_max]×[0.013π, 0.987π], where r_min < r_H. We enforced axial symmetry at the θ boundaries. Waves and particles are absorbed in a layer between r_pml = 30r_g and r_max = r_pml/0.9 ≃ 33r_g. The magnetic field lines in this layer are matched to the initial field.

The initial electromagnetic configuration followed the Wald (1974) solution with zero spin, corresponding to a uniform magnetic field at infinity of intensity B₀. We assumed that the black hole is initially surrounded by a stationary homogeneous cloud of thermal plasma with total number density n₀ and temperature T₀ as measured by FIDOs. The plasma is composed of electrons of mass m_e and charge −e and ions of mass m_i and charge +e, both with equal number densities n₀/2. To maintain a supply of fresh infalling plasma, we injected electrons and ions in a thin spherical shell between r_inj = r_pml − r_g and r_pml, enforcing a local number density floor of n ≥ n₀. This injected plasma is thermal with the same temperature T₀ as the initial cloud. We enforced roughly N_PPC = 10 particles per cell per species both throughout the domain at the simulation onset and at later times in the injection layer. We checked that the results are not quantitatively affected up to N_PPC = 50. Our initial and boundary conditions are summarized in Fig. 1.

Fig. 1. Sketch of the initial GRPIC setup representing zero-angular-momentum accretion onto a central black hole (labeled BH) in the center. On the left, we show in cyan the initial cloud of plasma, at a density n0 and temperature T0, and in pink the injection area of fresh plasma. On the right, the initial uniform magnetic field lines are shown in blue, and the limit of the outer matching layer in red.

Accretion can only occur if the ion thermal velocity, $v_{\mathrm{th}}=\sqrt{k_B{T}_0/m_i}$, where k_B is the Boltzmann constant, is lower than the gravitational escape velocity, $v_{\mathrm{esc}}=\sqrt{2r_g/r}$. This imposes a constraint on the normalized temperature, ϑ₀ ≡ k_BT₀/m_i, of ϑ₀ ≤ 2r_g/r_pml. To respect this limit, we set ϑ₀ = 1/30, which also means that the Bondi radius R_Bondi = 2𝒢M_BH/v_th² lies outside of our simulation box. In practice, we set the species-specific normalized temperatures, ϑ_0, e = k_BT_0, e/m_e and ϑ_0, i = k_BT_0, i/m_i, both equal to ϑ₀ = k_BT₀/m_i. For m_i ≠ m_e, this renders electrons initially colder than ions, but only reduces the overall plasma temperature from T₀ by at most a factor of two.

We targeted a low initial magnetization σ₀ = B₀²/4πρ₀ < 1, where ρ₀ = (m_i + m_e)n₀/2 is the plasma mass density. If σ₀ were too high in the injection zone, matter would not accrete due to the rigidity of magnetic field lines. We varied σ₀ for different simulations, but it is always of order 0.1. This also means that the plasma beta-parameter is β_plasma = 8πn₀k_BT/B² ≳ 1.

The main parameter we explored was the ion-to-electron mass ratio, m_i/m_e. A realistic mass ratio of m_i/m_e = 1836 imposes a broad separation between ion and electron scales and is therefore too computationally expensive to achieve. Instead, we employed more modest values of m_i/m_e, distinguishing between the accretion dynamics of a pair plasma, m_i/m_e = 1, and that of a progressively more realistic electron-ion plasma, m_i/m_e = 16 and 256.

The grid spacing in our simulations, parameterized through N_r and N_θ, was designed to resolve the plasma microscales, the smallest of which is the Debye length, ${\lambda }_D=\sqrt{k_B{T}_0/4\pi n_0{e}^2}$. Thanks to our logarithmically stretched radial grid, it suffices to resolve λ_D at the outer boundary of the box. The compression of the grid toward inner radii then compensates the plasma compression, keeping the kinetic scales resolved everywhere. We set N_r = N_θ in all of our simulations, adjusting N_r between 512 and 2048 according to the resolution requirements imposed by the values of m_i/m_e and σ₀ of a given run. In addition, we set the initial magnetic field strength B₀ such that the nominal ion gyroradius is r_L = m_i/eB₀ = r_g.

3. Accretion of plasma and flare dynamics

3.1. Global dynamics

Here we present the behavior of a representative pair-plasma simulation (m_i = m_e) with σ₀ = 0.05 and N_r = 1024, illustrated in Fig. 2. The top panel of Fig. 2 shows the magnetic flux threading the black hole, Φ_H, the accretion rate, Ṁ, and the bulk electromagnetic dissipation rate, $\dot{\mathcal{E}}$, defined as

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{\mathrm{H}}& =\frac{1}{2}{\displaystyle \int \int \sqrt{h}}|B^r|\text{d}\theta \text{d}\varphi ,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \dot{M}& =-{\displaystyle \int \int \sqrt{-g}}{N}^r\text{d}\theta \text{d}\varphi ,\text{and}\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \dot{\mathcal{E}}& ={\displaystyle \int \int \int \sqrt{h}}\mathbf{J}·\mathbf{E}\text{d}r\text{d}\theta \text{d}\varphi ,\hfill \end{array}$$(3)

Fig. 2. Top half: Time evolution of: the horizon magnetic flux, ΦH; the horizon accretion rate, Ṁ, normalized by ${\dot{M}}_0=4\pi m_i{n}_0{r}_{\mathrm{inj}}^2{v}_{\mathrm{th}}$; the normalized horizon magnetic flux, ${\mathrm{\Phi }}_{\mathrm{H}}/\sqrt{\dot{M}}$; and the total dissipation rate, $\dot{\mathcal{E}}$, normalized by Ṁ0. Bottom half: Simulation snapshots at t/tg = 2689, 3286, 3596,  and 3648. Each snapshot shows the number density n (top left), the electromagnetic energy dissipation J ⋅ E (top right), the magnetization σ (bottom left) and the local average particle Lorentz factor ⟨γ − 1⟩ (bottom right). Black lines show magnetic field lines.

where J is the electric current, h is the determinant of the 3-metric h_ij (Komissarov 2004), g is the determinant of the 4-metric g_μν, and N^r represents the rest-mass flux averaged over all particle species. Integrals (1) and (2) are evaluated at the horizon, r = r_H. The behavior in all of the above quantities is quasi-periodic, with each cycle comprising three distinct phases:

1. Φ_H and Ṁ both increase linearly in time;

2. The ratio ${\mathrm{\Phi }}_{\mathrm{H}}/\sqrt{\dot{M}}$ reaches a saturation value of roughly 60, remaining approximately constant as magnetic flux continues to accrue on the horizon;

3. A critical point is reached, and the system undergoes an eruption in which the magnetic flux decreases by a factor of ∼4 over a timescale of ∼100t_g – a very short time compared to the ∼10³t_g cycle period.

The eruption is accompanied by a sharp rise in $\dot{\mathcal{E}}$, indicating the dissipation of a substantial fraction of the magnetic energy stored during the preceding flux accumulation phases. After the ejection, the cycle restarts, with slow flux accumulation culminating in further abrupt eruptions. All of our simulations exhibit the three phases outlined above, and we refer to these phases – one, two, and three – throughout this work. The first cycles in Fig. 2 exhibit deviations from the others both in duration and peak magnetic flux storage. This may reflect memory effects associated with the initial conditions. Starting from the fourth cycle, however, the system settles into a more quasi-stationary pattern, with successive cycles exhibiting increasingly uniform behavior.

Four snapshots of the simulation are shown in Fig. 2. Three of these snapshots represent the three phases of the simulation cycle; one represents the moment where Φ_H peaks, from which the eruption (phase three) is triggered. Each snapshot displays spatial maps of the plasma number density n, the magnetization σ = B²/4πρ, the averaged particle Lorentz factor ⟨γ⟩, and the electromagnetic energy dissipation J ⋅ E. In the spatial maps of Fig. 2, the magnetic field divides the box into two topologically distinct zones. One zone is threaded by field lines attached to the black hole. These field lines form a highly magnetized funnel region where plasma accretes by flowing along the magnetic field. The other zone is threaded by field lines not attached to the black hole. These field lines are held in place by a dense and weakly magnetized plasma that builds up along the equator. Unlike in the funnel region, the dense equatorial plasma accretes mostly perpendicular to the field lines, carrying them toward the black hole and, thus, building up the horizon-threading magnetic flux and corresponding funnel zone. As the funnel broadens, oppositely oriented field lines are pushed closer to the equator, building up a radial field reversal. This reversal is supported by an azimuthal electric current that develops in the equatorial plasma.

Starting in phase two, the equatorial current layer thins to the point where intermittent magnetic reconnection events appear, contributing to the regulation of the magnetic flux threading the black hole (second snapshot of Fig. 2). These events are too minor in phase two to repulse the heavy plasma channeled in through the equator. However, during phase three, a much more powerful reconnection event (third snapshot of Fig. 2) completely arrests and expels the accretion flow (fourth snapshot of Fig. 2), preying on the energy stored in the accumulated magnetic field to do so. This energy reservoir depletes once a significant fraction of Φ_H is consumed. At that point, the electromagnetic force can no longer resist the force of gravity, and material infall resumes.

To quantify the energetic efficiency of this process, we characterize key time-averaged values of Ṁ and $\dot{\mathcal{E}}$. We measure:

$$\begin{array}{cc}\hfill {⟨\dot{M}⟩}_{\mathrm{cycle}}& =\frac{\int \dot{M}\text{d}t}{\int \text{d}t}=0.14{\dot{M}}_0,\hfill \\ \hfill {⟨\dot{\mathcal{E}}⟩}_{\mathrm{cycle}}& =\frac{\int \mathrm{\Theta }\left(\dot{\mathcal{E}}\right)\dot{\mathcal{E}}\text{d}t}{\int \text{d}t}=0.0044{\dot{M}}_0,\mathrm{and}\hfill \\ \hfill {⟨\dot{\mathcal{E}}⟩}_{\mathrm{flare}}& =\frac{\int \mathrm{\Theta }\left(\dot{\mathcal{E}}\right)\dot{\mathcal{E}}\text{d}t}{\int \mathrm{\Theta }\left(\dot{\mathcal{E}}\right)\text{d}t}=0.021{\dot{M}}_0,\hfill \end{array}$$(4)

where, as in Fig. 2, ${\dot{M}}_0=4\pi m_i{n}_0{r}_{\mathrm{inj}}^2{v}_{\mathrm{th}}$. In the above, integrals are taken over the fourth and fifth accretion-expulsion cycles (last two cycles in the time series of Fig. 2) and the Heaviside function, Θ(x), selects only positive values of $\dot{\mathcal{E}}$ (red regions in the $\dot{\mathcal{E}}$-panel of Fig. 2). Positive $\dot{\mathcal{E}}$ corresponds to electromagnetic energy transferred to particles and thus represents energy made available for emission as observable radiation. Instead, negative $\dot{\mathcal{E}}$ represents the work of the particles (and ultimately gravity) increasing the electromagnetic energy density within our simulation box. We can use these measurements to define the black hole efficiency factors

$$\begin{array}{c}\hfill {\eta }_{\mathrm{cycle}}\equiv \frac{{⟨\dot{\mathcal{E}}⟩}_{\mathrm{cycle}}}{{⟨\dot{M}⟩}_{\mathrm{cycle}}}\simeq 0.03\mathrm{and}{\eta }_{\mathrm{flare}}\equiv \frac{{⟨\dot{\mathcal{E}}⟩}_{\mathrm{flare}}}{{⟨\dot{M}⟩}_{\mathrm{cycle}}}\simeq 0.2.\end{array}$$(5)

These efficiency factors are useful in different observational situations. For a dim source, one may wish to integrate over many flare cycles to obtain a significant detection. In that case, η_cycle⟨Ṁ⟩_cycle can be used to estimate the exposure time needed to accumulate a required number of photons. However, for a bright source where observations of single flares are possible, η_flare⟨Ṁ⟩_cycle provides an estimate for the on-flare luminosity. We discuss use cases of both η_cycle and η_flare in Section 6.

3.2. Particle acceleration in the different phases

In addition to regulating accretion, reconnection results in particle acceleration visible in the J ⋅ E and ⟨γ⟩−1 maps of Fig. 2. The intermittent reconnection events of phase two correspond to small flashes of particle acceleration that occur near the horizon, although the accelerated particles are ultimately swallowed by the black hole. During the erupting phase, a substantial fraction of the energy built up in the magnetosphere is abruptly transferred to particles, resulting in much more intense equatorial particle acceleration over a broader range of radii. Some of the accelerated particles escape along the wall of the funnel region (fourth snapshot of Fig. 2).

These observations are corroborated by the evolution of the particle energy distribution, shown in Fig. 3. The figure shows instantaneous particle distributions as well as distributions obtained by time-averaging over phases one, two, and three of one cycle of the simulation. All distributions show a thermal bump with a maximum at $\gamma -1\sim 0.1\sim \sqrt{{\vartheta }_0}$, corresponding to the temperature of injected particles. Particle acceleration manifests as high-energy nonthermal tails in the distributions.

Fig. 3. Top panel: Particle energy distributions at different stages of the eruption cycle. Transparent lines show instantaneous distributions; opaque lines represent distributions time-averaged over each phase. Phases one, two, and three are color-coded, respectively, as green, blue, and red. Bottom panel: Time evolution of ΦH over one eruption cycle. Shaded regions (green, blue, red) indicate the time intervals during which particle distributions are presented in the top panel.

In phase one (green), the particle spectrum appears rather soft with little variability. The absence of magnetic reconnection from this phase results in little to no particle acceleration. In contrast, during the reconnection-regulated phase (blue), intermittent reconnection accelerates particles episodically, resulting in a somewhat extended and highly variable tail in the distributions. The erupting phase (red) shows the strongest particle acceleration, with more particles being accelerated overall and with a higher maximum particle energy. A large number of magnetic reconnection studies in 2D slab geometry show that, in pair plasmas, the particle energy distribution cuts off (or at least turns over) around γ ∼ σ, corresponding to an approximately equal partitioning of the available magnetic energy density among the accelerated particles: γnm_ec² ∼ B²/8π (see Sironi et al. 2025 and references therein). During the eruption phase, the particle distribution turns over at around γ ∼ 10, hinting that the magnetization witnessed by particles in the reconnection layer is σ ∼ 10, consistent with the values in the funnel region shown in Fig. 2.

3.3. Robustness of the dynamics

All of the above phases and particle acceleration results (with the exception of the separation between species that develops for m_i > m_e, described later on) hold for all of the simulations presented in this work. However, we also conducted a broad suite of supplementary test simulations, pushing the numerical parameters beyond the confines of those presented here. Below, we comment on the robustness of the results described so far in light of these test campaigns.

• Effect of σ₀: We explored initial magnetizations ranging from σ₀ = 0.3 to σ₀ = 0.03. We observed that the magnetization on the event horizon, σ_H, reaches ∼10 in all cases. In addition, the maximum value of Φ_H (reached immediately before eruption) decreases with increasing σ₀. Less magnetic flux accumulation is required to reach σ_H = 10 for a higher σ₀. However, if the simulation box is not big enough, σ_H may not reach 10, as described below.

• Effect of box size: Reducing r_max by 2/3 resulted in a decreased maximum value of Φ_H with respect to the simulation presented above, with σ_H no longer saturating near 10. On the contrary, increasing r_max by a factor of 5/3 altered neither the maximum Φ_H nor the saturated σ_H value. We interpret this as meaning that our runs are converged, possessing enough initial magnetic flux in the box to supply the flux Φ_H necessary to achieve a horizon magnetization of σ_H ∼ 10.

• Effect of initial temperature: We did not alter ϑ₀ = k_BT₀/m_i much, keeping it close to the Bondi limit, 2r_g/r_pml (see Sect. 2). However, we noticed that plasma instabilities such as mirroring arose in the funnel region for simulations with σ < 0.1 and β_plasma > 1 (a regime barely accessed by our runs). This feature has also been observed in the context of black hole accretion by Galishnikova et al. (2023).

In addition, we note that our simulations consider an unrealistically small Bondi radius, which lies at the edge of our simulation domain (∼30r_g). Modeling realistic Bondi accretion is not computationally accessible with kinetic simulations. However, many of the essential features we observe, including magnetic flux eruptions, have also been witnessed in GRMHD simulations with a realistic material supply from larger scales (Ressler et al. 2021). In light of this broad agreement with prior GRMHD work, and also thanks to the robustness of our results to variations in numerical parameters discussed here, we believe that the reduced-scale dynamics observed in our kinetic simulations may scale up to larger systems.

4. Effect of scale separation between the species

4.1. Scale separation features

In this second part of our work, we study the effect of the mass ratio, m_i/m_e, on the dynamics. We present simulations with m_i/m_e ∈ [1, 16, 256], increasing the initial magnetization to σ₀ = 0.3 to make this range numerically feasible. All cases exhibit the same three-phase cyclic dynamics as detailed for pair plasmas in Section 3. The typical eruption timescale is still ∼100t_g and σ_H still peaks at ∼10. However, due to the increased σ₀ (see Section 3.3), the maximum value of Φ_H is somewhat reduced, reaching ∼40 − 50. Overall, the marked similarity to the pair-plasma case suggests that pair plasmas may be used to model eruptive magnetospheric dynamics at substantially lower computational cost.

However, a nonzero inter-species scale separation does impact certain microscopic details of the evolution. First, due to the higher ion inertia, ions show a somewhat different distribution of equatorial current than electrons. As shown in the top panels of Fig. 4, the azimuthal current carried by the electrons is more intense and organized on smaller spatial scales, whereas the ion current is more diffuse. The electron-ion scale separation is also highlighted at reconnection X-points. Electrons approach X-points more closely than ions do before becoming demagnetized and crossing into the outflow region. This decoupling between species generates a poloidal electric current in the outflow from X-points, inducing a quadrupolar magnetic field H_ϕ (Werner et al. 2018), as shown in the bottom panel of Fig. 4.

Fig. 4. Top: map of the azimuthal electric current carried by the electrons (Jϕ−). Center: map of the azimuthal electric current carried by the ions (Jϕ+). Bottom: map of the auxiliary azimuthal magnetic field Hϕ. Magnetic field lines are represented in solid black lines for all panels. The mass ratio is mi/me = 256 in the depicted simulation.

4.2. Impact on particle acceleration

Perhaps the starkest difference between electrons and ions is with respect to particle energization. As previously explained, magnetic reconnection is the key mechanism that accelerates particles during flux eruptions. If the liberated energy were divided evenly among all particles, then the electrons would reach much higher Lorentz factors than the ions. Even though the species do not share energy precisely equally in reality (Werner et al. 2018; Werner & Uzdensky 2024; Comisso 2024), this is a small effect compared to the mass ratio in deciding the characteristic Lorentz factors that particles attain. For example, Fig. 5 displays a map of the average kinetic energy for each species during a magnetic flux eruption. The strong discrepancy between electron and ion energies in the funnel region is merely an initialization artifact; electrons are injected colder than ions (Section 2). Despite the initial temperature gap, the species achieve comparable average energies once processed by reconnection, indicating characteristic Lorentz factors of accelerated electrons a factor of m_i/m_e higher than those of accelerated ions.

Fig. 5. Map of the local average particle energy for electrons (left) and ions (right) for the simulation with mi/me = 256. The energies are normalized to mec2 and grey lines represent magnetic field lines.

Here again, we see that the spatial distribution of accelerated ions is more diffuse than that of high-energy electrons. High-energy electrons reside both in the equatorial reconnection layer and along the wall of the funnel region. This suggests both locations as possible sites of high-energy emission during flux eruption events.

We also compare the particle energy distributions of the species as a function of mass ratio in Fig. 6. The figure shows snapshots taken during the third erupting phase for all the mass ratios considered. The ion energy distributions remain virtually identical when increasing the mass ratio; the only difference between ion distributions is a horizontal shift consistent with their higher rest-mass energy. This again demonstrates the robustness of the overall dynamics with respect to the mass ratio: the species dominating the plasma inertia behaves similarly in all cases.

Fig. 6. Time-averaged energy distributions of the particles during the erupting phase for all electron-ion simulations (including the reference mi = me run). Ions are represented as dashed lines and electrons as solid lines.

However, qualitative differences emerge for the electrons as m_i/m_e is increased. The extent of the nonthermal tail of the electron distribution is enhanced, covering almost two decades for m_i/m_e = 256 instead of barely one for m_i = m_e. The power-law slope also appears to harden for higher mass ratios, reaching an index of dlog N/dlog γ ∼ −1.7 at the highest mass ratio. We attribute both effects to the increased magnetization felt by the electrons near the event horizon, σ_H, e ∼ (m_i/m_e)σ_H ∼ 10³. With more magnetic energy available per unit of rest-mass energy, the tail of the electron distribution turns over at a higher Lorentz factor ∼σ_H, e. The hardening of the power-law slope with increasing magnetization probably owes to a similar mechanism as has been previously observed in 2D simulations of pair-plasma (Guo et al. 2014; Werner et al. 2016) and electron-ion plasma (Werner et al. 2018) reconnection in slab geometry.

This robust trend towards stronger nonthermal electron acceleration with increasing mass ratio can be extrapolated to the full proton-electron value, m_i/m_e = 1836. We expect that, in this case, the electron energy distribution should turn over or cut off at γ ∼ σ_H, e ∼ (m_i/m_e)σ_H ∼ 10⁴. However, we expect the power-law slope to plateau in the ultrarelativistic limit (Guo et al. 2014; Werner et al. 2016, 2018; Ball et al. 2018) and, thus, to remain roughly constant between m_i/m_e = 256 (corresponding to σ_H, e ∼ 10³) and m_i/m_e = 1836 (for which σ_H, e ∼ 10⁴).

5. Modeling the cyclic activity

In what follows, we present several theoretical aspects of the three-phase cycle outlined in Section 3.1. We discuss major physics points of phases one, two, and three in Sections 5.1, 5.2, and 5.3, respectively. Then, in Section 5.4, we argue that the transitions between phases are governed by the progressive thinning of the equatorial current layer.

5.1. Phase I: Ideal accretion

This is the first of two phases where the plasma accretes onto the black hole, dragging magnetic field lines with it. The hallmark of this phase is that it respects the ideal flux-freezing condition of magnetohydrodynamics (MHD). The accretion of flux is thus tied to that of matter, as we demonstrate below.

To show how material and magnetic flux transport are linked, we here derive and solve the magnetic flux transport equation. To arrive at this equation, we need three ingredients: (i) the ideal MHD condition, which can be expressed as

$$\begin{array}{c}\hfill \mathbf{E}+\mathbf{V}\times \mathbf{B}=0\end{array}$$(6)

in the 3 + 1 formalism of Komissarov (2004), where V is the bulk fluid 3-velocity; (ii) the magnetic flux threading a spherical cap of radius r and half opening angle θ,

$$\begin{array}{c}\hfill \mathrm{\Phi }(r,\theta ,t)=2\pi {\displaystyle {\int }_0^{\theta }}\sqrt{h}{B}^r\mathrm{d}\theta ;\end{array}$$(7)

and (iii) the Maxwell-Faraday law

$$\begin{array}{c}\hfill {\partial }_t\mathbf{B}=-\mathrm{\nabla }\times \mathbf{E}.\end{array}$$(8)

Integrating Eq. (8) over the same spherical cap as Eq. (7) yields the time-evolution equation for Φ:

$$\begin{array}{c}\hfill {\partial }_t\mathrm{\Phi }=-2\pi {\displaystyle {\int }_0^{\theta }}{\partial }_{\theta }{E}_{\varphi }\mathrm{d}\theta =-2\pi E_{\varphi }.\end{array}$$(9)

Finally, plugging in the ϕ-component of Eq. (6),

$$\begin{array}{c}\hfill E_{\varphi }+\sqrt{h}(V^r{B}^{\theta }-V^{\theta }{B}^r)=0,\end{array}$$(10)

and using

$$\begin{array}{c}\hfill 2\pi \sqrt{h}{B}^r={\partial }_{\theta }\mathrm{\Phi }\mathrm{and}2\pi \sqrt{h}{B}^{\theta }=-{\partial }_r\mathrm{\Phi },\end{array}$$(11)

gives the equation for magnetic flux transport,

$$\begin{array}{c}\hfill {\partial }_t\mathrm{\Phi }=2\pi \sqrt{h}(V^r{B}^{\theta }-V^{\theta }{B}^r)=-V^r{\partial }_r\mathrm{\Phi }-V^{\theta }{\partial }_{\theta }\mathrm{\Phi }.\end{array}$$(12)

We now specialize to the equatorial plane, θ = π/2, for which the top-down symmetry of our problem dictates V^θ = 0. This results in the simplified flux-advection equation

$$\begin{array}{c}\hfill ({\partial }_t+V^r{\partial }_r)\mathrm{\Phi }=0.\end{array}$$(13)

The equatorial material accretion velocity, V^r, dictates magnetic flux deposition onto the black hole. Conversely, if the flux transport, Φ(r, π/2, t), is known, then the equatorial accretion velocity follows. Given that we have measured empirically that Φ(r_H, π/2, t) = Φ_H grows linearly in time, we investigate which velocity profiles, V^r, are consistent with ${\dot{\mathrm{\Phi }}}_{\mathrm{H}}=\mathrm{constant}$. We consider self-similar, time-independent profiles V^r = −𝒞/r^ξ, with proportionality constant 𝒞 > 0. Then, we solve Eq. (13) by the method of characteristics, searching for the ξ that reproduces ∂_tΦ_H = const.

For the ansatz, V^r = −𝒞/r^ξ, the solution to Eq. (13) is

$$\begin{array}{c}\hfill \mathrm{\Phi }(r,\pi /2,t)=\mathrm{\Phi }(r_0,\pi /2,0)=\pi B_0{r}_0^2,\end{array}$$(14)

where

$$\begin{array}{c}\hfill r_0={[\mathcal{C}(1+\xi )t+r^{1+\xi }]}^{1/(1+\xi )},\end{array}$$(15)

and we have used our simulation initial condition, Φ(r₀, π/2, 0) = πB₀r₀². We search to match the solution at the horizon, r = r_H, to a linearly growing function of t. For r = r_H and $t\gg r_{\mathrm{H}}^{1+\xi }/[\mathcal{C}(1+\xi )]$, we have

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{H}}=\mathrm{\Phi }(r_{\mathrm{H}},\pi /2,t)\simeq \pi B_0{[\mathcal{C}(1+\xi )t]}^{2/(1+\xi )}.\end{array}$$(16)

Thus, for a self-similar, time-independent velocity profile V^r, the flux threading the horizon, Φ_H, grows linearly in time if and only if ξ = 1. Rewriting the constant 𝒞 in terms of ${\dot{\mathrm{\Phi }}}_{\mathrm{H}}$ then gives

$$\begin{array}{c}\hfill V^r=-\frac{{\dot{\mathrm{\Phi }}}_{\mathrm{H}}}{2\pi B_{0}r}·\end{array}$$(17)

Next, we examine how well the profile (Eq. (17)) matches the equatorial accretion velocity of our simulations. We present in Fig. 7 the equatorial profiles of V^r measured from our simulation presented in Section 3. To compare with Eq. (17), we plug in ${\dot{\mathrm{\Phi }}}_{\mathrm{H}}\simeq 0.20$, measured during phase one in Fig. 3.

Fig. 7. Profiles of Vr measured from the simulation presented in Section 3. Solid lines denote time averages over the intervals defined in Fig. 3; envelopes denote one-sigma percentiles for each averaging interval. The orange dashed curved corresponds to Eq. (17) where ${\dot{\mathrm{\Phi }}}_{\mathrm{H}}$ is measured during phase one on Fig. 3. The dashed black line denotes the free-fall velocity of a particle starting from rest at infinity (Eq. (18)).

The measured V^r profiles agree very well with the theoretical prediction (Eq. (17)) far from the black hole. Close to the black hole, there is some deviation, suggesting that some amount of magnetic diffusivity may be present near the horizon even before reconnection starts to act. Importantly, however, the accretion velocities are all much slower than that of a particle that freely falls from rest at infinity,

$$\begin{array}{c}\hfill V_{\mathrm{ff}}^r=-{\left(\frac{r}{r_{\mathrm{H}}}\right)}^{-1/2}\frac{1-{(r/r_{\mathrm{H}})}^{-1}}{1-{(r/r_{\mathrm{H}})}^{-3/2}},\end{array}$$(18)

indicating that magnetic tension is strong enough to significantly slow gravitational infall even in phase one.

5.2. Phase II: Reconnection-regulated accretion

Matter and magnetic fields continue to accrete in this phase. The crucial difference from the preceding phase is the appearance of intermittent near-horizon reconnection events. These are traced by the high-frequency noise in the time series of Φ_H in this phase (Fig. 3). Because reconnection allows particles to locally slip across magnetic field lines, less flux accrues on the black hole per unit of mass accreted. In fact, the average equatorial accretion velocity V^r, as shown in Fig. 7, is not that different in this phase from phase one. Nevertheless, reconnection-mediated magnetic diffusion reduces the average growth rate of Φ_H by a factor of ≃2 with respect to the first phase, as shown in Fig. 3.

Another distinguishing aspect of this phase is that the magnetic field near the horizon strongly regulates accretion. Namely, the magnetic tension across the equator approximately matches the force of gravity on the heavy equatorial plasma. This fixes the ratio of ${\mathrm{\Phi }}_{\mathrm{H}}/\sqrt{\dot{M}}$ to roughly 60, as we derive below.

In the 3 + 1 formalism of Komissarov (2004), FIDOs measure a gravitational pull on the plasma of

$$\begin{array}{c}\hfill f_r^G=-\rho \frac{{\partial }_r\alpha }{\alpha },\end{array}$$(19)

where $\alpha =1/\sqrt{1+2r_g/r}$ is the lapse function, ρ is the FIDO-measured mass density, and we assume a nonrelativistic plasma. The expression for the Lorentz force measured by FIDOs is

$$\begin{array}{c}\hfill f_r^L=-\sqrt{h}\frac{J^{\varphi }}{\alpha }{B}^{\theta },\end{array}$$(20)

where

$$\begin{array}{c}\hfill J^{\varphi }=\frac{1}{4\pi }{(\mathrm{\nabla }\times H)}^{\varphi }=-\frac{1}{4\pi \sqrt{h}}{\partial }_{\theta }{H}_r.\end{array}$$(21)

The auxiliary magnetic field, H_r, is related to B_r by H_r = αg_rrB^r = B^r/α. In the vicinity of an X-point, B^r switches sign from B_up^r to −B_up^r over an angular width Δθ = 2δ/r_X. Here, r_X is the radial position of the X-point and δ is the local half-thickness of the current sheet. Plugging this into Eq. (21) gives

$$\begin{array}{c}\hfill J^{\varphi }=\frac{1}{4\pi \alpha \sqrt{h}}\frac{r_X{B}_{\mathrm{up}}^r}{\delta }·\end{array}$$(22)

At the X-point, a vertical magnetic field component B_X^θ is created by the reconnection of the upstream magnetic field B_up^r. In the orthonormal basis ($B^{\widehat{i}}=\sqrt{g_{\mathit{ii}}}{B}^i$), we expect that the reconnection rate, β_rec, relates the r and θ components of the magnetic field by

$$\begin{array}{c}\hfill {\beta }_{\mathrm{rec}}\sim \frac{B_X^{\widehat{\theta }}}{B_{\mathrm{up}}^{\widehat{r}}}=\sqrt{\frac{g_{\theta \theta }}{g_{\mathit{rr}}}}\frac{B_X^{\theta }}{B_{\mathrm{up}}^r}=r_X\alpha \frac{B_X^{\theta }}{B_{\mathrm{up}}^r},\end{array}$$(23)

Using Eqs. (20), (22), and (23), we recover the Lorentz force estimate,

$$\begin{array}{c}\hfill f_r^L\sim \frac{{\beta }_{\mathrm{rec}}}{4\pi {\alpha }^3\delta }{\left(B_{\mathrm{up}}^r\right)}^2.\end{array}$$(24)

Balancing the gravitational and Lorentz forces then gives

$$\begin{array}{c}\hfill \frac{{\beta }_{\mathrm{rec}}}{4\pi {\alpha }^3\delta }{\left(B_{\mathrm{up}}^r\right)}^2\sim \rho \frac{{\partial }_r\alpha }{\alpha }·\end{array}$$(25)

We approximate B_up^r ∼ αΦ_H/2πr_X² and assume that a significant part of the accretion flow passes through the current layer (which we verified in our simulations) so that Ṁ(r = r_H) ∼ Ṁ(r = r_X) ∼ ρV_acc 4πr_Xδ/α. In this manner, we replace B^r with Φ_H and ρ with Ṁ to get

$$\begin{array}{c}\hfill \frac{V_{\mathrm{acc}}{\beta }_{\mathrm{rec}}}{\alpha }\frac{{\mathrm{\Phi }}_{\mathrm{H}}^2}{4{\pi }^2{r}_X^3}=\dot{M}{\partial }_r\alpha .\end{array}$$(26)

Rearranging, we arrive at

$$\begin{array}{c}\hfill \frac{{\mathrm{\Phi }}_{\mathrm{H}}}{\sqrt{\dot{M}}}=2\pi r_X^{3/2}\sqrt{\frac{\alpha {\partial }_r\alpha }{V_{\mathrm{acc}}{\beta }_{\mathrm{rec}}}}·\end{array}$$(27)

During phase two, the X-point location is observed in the simulations to lie on average at r_X ∼ 3r_g, and the accretion velocity is measured to be V_acc ∼ 0.1. Plugging these in and assuming β_rec ∼ 0.1, we find that

$$\begin{array}{c}\hfill \frac{{\mathrm{\Phi }}_{\mathrm{H}}}{\sqrt{\dot{M}}}\sim 60.\end{array}$$(28)

This derivation was made using strong assumptions on the dynamics and the accretion flow, and should be explored in the context of more realistic accretion conditions. However, the fact that our simplified setup of accretion and this analytical scaling match the larger-scale GRMHD estimates (e.g. Igumenshchev & Narayan 2002; Ressler et al. 2021) indicates that it might be robust upon more realistic accretion flow conditions.

Conversely, the above arguments cannot predict σ_H. The latter depends on the mass density within the funnel, whereas the model in this section constrains only the mass density near the equator. We expect to conduct a more comprehensive theoretical and numerical exploration of σ_H in future work and, in particular, to examine its dependence on the black hole spin.

5.3. Phase III: Flux eruption

The delicate balance between magnetic tension and gravity that characterizes the second phase sets the stage for the eruptive phase. This is because reconnection tends to convert horizontal field (B^r) into vertical field (B^θ). Thus, when a sudden, strong reconnection event occurs, it rapidly builds up B^θ, disrupting the balance between magnetic tension and gravity and triggering an eruption. During the eruption, a macroscopic reconnection zone forms at the equator, removing flux from the black hole and expelling it outward. The rate at which magnetic flux decays off of the horizon is governed by the collisionless relativistic reconnection rate, β_rec ∼ 0.1, and follows an exponential decay law (Fig. 3). Following models developed by Crinquand et al. (2021) and Bransgrove et al. (2021) we derive this decay law below.

FIDOs just upstream of the reconnection layer (i.e., slightly above or below the equator) measure a ratio of the electric-to-magnetic field equal to the reconnection rate:

$$\begin{array}{c}\hfill \frac{D^{\widehat{\varphi }}}{B^{\widehat{r}}}\sim {\beta }_{\mathrm{rec}},\end{array}$$(29)

where $D^{\widehat{i}}=\sqrt{g_{\mathit{ii}}}{D}^i$ and β_rec is the reconnection rate expectation in flat spacetime. Assuming a split monopolar magnetic geometry near the event horizon, and placing the FIDO above the X-point, we have

$$\begin{array}{c}\hfill D^{\varphi }={\beta }_{\mathrm{rec}}\sqrt{\frac{g_{\mathit{rr}}}{g_{\varphi \varphi }}}{B}^r={\beta }_{\mathrm{rec}}\sqrt{\frac{g_{\mathit{rr}}}{g_{\varphi \varphi }}}\frac{\alpha {\mathrm{\Phi }}_{\mathrm{H}}}{2\pi r_X^2}·\end{array}$$(30)

Neglecting any B^θ, E_ϕ is related to D^ϕ by

$$\begin{array}{c}\hfill E_{\varphi }\sim \alpha g_{\varphi \varphi }{D}^{\varphi }.\end{array}$$(31)

Given that ∂_tΦ = −2πE_ϕ (Eq. (9)), g_rr = 1/α², and g_ϕϕ = r²sin²θ, we get

$$\begin{array}{c}\hfill {\partial }_t\mathrm{\Phi }=-{\beta }_{\mathrm{rec}}\frac{\alpha sin\theta }{r_X}{\mathrm{\Phi }}_{\mathrm{H}}.\end{array}$$(32)

We assume that flux loss is transferred from r_X to r_H much faster than it reconnects at r_X, such that we may write Φ(r_X, π/2)≃Φ(r_H, π/2) = Φ_H. This gives,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{H}}}{\mathrm{d}t}=-\frac{\alpha {\beta }_{\mathrm{rec}}}{r_X}{\mathrm{\Phi }}_{\mathrm{H}},\end{array}$$(33)

an exponential decay law that matches the behavior in Fig. 3. In that figure, we fit a decay rate of αβ_rec/r_X ∼ 0.013t_g⁻¹, which implies, for r_X ∼ 3r_g, that β_rec ∼ 0.05. This is a factor of two lower than the expectation for relativistic magnetic reconnection in flat spacetime. We emphasize, however, that this estimate of β_rec should be seen as an average over the whole eruption, during which the upstream plasma properties (such as the magnetization) change. We believe that the evolving upstream, as well as the spherical geometry and the presence of a vertical magnetic field, may account for the slightly slower β_rec measured here than in typical slab-geometry simulations.

5.4. Phase transitions dictated by current layer stability

Sections 5.1–5.3 underscore the role of reconnection in distinguishing each phase. As the layer thins, it becomes progressively more prone to tearing (Zelenyi & Krasnoselskikh 1979; Zenitani & Hoshino 2008), the formation of X-points and plasmoids, and, thus, to the rapid reconnection of magnetic flux (Shibata & Tanuma 2001; Loureiro et al. 2007; Uzdensky et al. 2010). We argue in this section that the transitions between phases correspond to critical tearing instability thresholds of the current sheet.

To begin with, we measure the half-thickness, δ, of the equatorial current layer as a function of radius, r, and time, t. These measurements are presented in Fig. 8. In the figure, we normalize δ(r) by r − r_H, the latter being a proxy for the length of the current sheet. Figure 8 suggests that phase one corresponds to a half opening angle, δ(r)/(r − r_H), greater than 0.1. In phases two and three, the half-opening angle is roughly 0.1 and 0.05, respectively.

Fig. 8. Traces of δ(r)/(r − rH) at specific times within (transparent lines), and time-averaged over (solid dashed lines), each phase.

These critical opening angles correspond to tearing instability thresholds of the equatorial current sheet. The rate of collisionless plasmoid-mediated relativistic reconnection is β_rec ≃ 0.1. This corresponds to the aspect ratio of elementary current layers between small-scale plasmoids. On average, once the distance between two plasmoids exceeds ${\beta }_{\mathrm{rec}}^{-1}\simeq 10$ times the width of the current sheet between them, that current sheet tears and spins off another plasmoid (Shibata & Tanuma 2001; Loureiro et al. 2005, 2007; Uzdensky et al. 2010; Cerutti et al. 2014).

When δ/(r − r_H) = 0.1, the magnetic field above and below the equatorial current sheet opens at the same angle as the separatrix field lines near a single reconnection X-point. This is conducive to the formation of single X-points, which sporadically appear just beyond the horizon in this regime. They are, however, short-lived, falling quickly into the black hole. Thus, in phase two, the current sheet episodically undergoes single-X-point reconnection (Ji & Daughton 2011), with the X-point formation site being forced by the global field geometry to lie close to the event horizon.

Once δ/(r − r_H) decreases down to 0.05, a new, multiple X-point regime is attained. This is because the equatorial current layer is now so thin that its full width, 2δ, becomes 10 times smaller than its full length, r − r_H. The layer thus becomes statistically likely to form X-points far from the horizon, developing a rapidly reconnecting plasmoid chain (Shibata & Tanuma 2001; Uzdensky et al. 2010; Ji & Daughton 2011). Whereas the exhaust from a single X-point is characterized by a relatively small vertical magnetic field, $B^{\widehat{\theta }}\sim {\beta }_{\mathrm{rec}}{B}^{\widehat{r}}$, in this new fully nonlinear stage of the tearing instability, characterized by circular plasmoids, the vertical field becomes comparable to the upstream, unreconnected field: $B^{\widehat{\theta }}\sim B^{\widehat{r}}$. This amplification upsets the equilibrium between magnetic tension and gravity that previously held during phase two, thus triggering the eruption phase.

These ideas are summarized in Fig. 9. The transition from phase one to phase two corresponds to the moment when the half-opening angle at the horizon reaches θ_open ≃ β_rec = 0.1, triggering single X-point reconnection. The subsequent transition to phase three occurs once the layer further thins to the point, θ_open = β_rec/2 = 0.05, marking the transition to the multiple X-point, plasmoid-dominated reconnection regime.

Fig. 9. Conceptual three-phase accretion cycle. Transitions between phases coincide with critical opening angles of the magnetic field lines about the equator. Blue and green field lines are accreted during phases one and two, respectively; black field lines never reach the black hole.

6. Astrophysical applications

Generally, we only expect strong nonthermal radiation signatures from the third (eruptive) phase of each flux accumulation cycle. Reconnection does not act in the first phase, and so nonthermal particle acceleration and concomitant nonthermal radiation should be virtually absent. We also expect the second phase to produce very little high-energy emission. Although intermittent reconnection occurs in this phase, it happens very close to the horizon with the emitting particles primarily on infalling trajectories. This presents unfavorable conditions for high-energy radiation to escape the magnetosphere. In contrast to these first two phases, the erupting phase exhibits the strongest reconnection, with a significant fraction of the accelerated particles lying farther away from the black hole and moving outward. We therefore expect most of the observable high-energy nonthermal emission to arise during this phase.

Next, we outline specific astrophysical implications of our results with respect to two different types of low-luminosity black holes. We focus primarily on infrared (IR) and X-ray flares observed from Sgr A^★. Secondarily, we speculate on radiation that may be produced by a mechanism similar to our model from isolated stellar-mass black holes interacting with the Galactic interstellar medium (ISM).

6.1. Sgr A^★ IR and X-ray flares

We organize our discussion around the following observed characteristics of IR and X-ray flaring in Sgr A^★:

• Duration: The eruption phase lasts ∼100t_g in all our simulations – a robust timescale tied to the universal rate at which magnetic reconnection removes flux from the black hole (Section 5.3). Adopting the mass of Sgr A^★, ∼4 × 10⁶ M_⊙, this corresponds to a flare duration of ∼30 min, broadly consistent with the timescales of the X-ray and IR Sgr A^★ flares (Genzel et al. 2003; Ghez et al. 2004; Porquet et al. 2008; GRAVITY Collaboration 2021; von Fellenberg et al. 2025).

• Recurrence timescale: Our simulations show a cycle period of about 10³t_g. Recently, Jacquemin-Ide et al. (2025) argued, in the context of GRMHD simulations, that this timescale arises from the balance of magnetic flux advection and diffusion near the black hole, a balance reminiscent of how reconnection-induced diffusivity regulates flux accumulation in phase two of our simulations (Section 5.2). For Sgr A^★, 10³t_g corresponds to a flare period of ∼6 h, of order the characteristic near-IR recurrence timescale (Meyer et al. 2009; Witzel et al. 2018). However, this identification is still speculative at this stage; more theoretical work is needed to determine the dependence, if any, of the recurrence time on other aspects such as the dimensionality or the mass, angular momentum, and flux supply from larger scales.

• Luminosity: Given typical horizon-scale magnetic field strengths inferred for Sgr A^★, generally in the range 10 − 100 G (Loeb & Waxman 2007; Event Horizon Telescope Collaboration 2024; von Fellenberg et al. 2025), we can use the flux saturation condition, $\varphi ={\mathrm{\Phi }}_{\mathrm{H}}/\sqrt{{\dot{M}}_{\mathrm{H}}}\simeq 60$, to compute a characteristic flare luminosity in the context of our model. Since this condition applies in phase two during which mass accretes steadily, we estimate the mass accretion rate using its cycle-averaged value, Ṁ_H ≃ ⟨Ṁ⟨_cycle. Then, to tie this rate to the on-flare electromagnetic heating of particles, ${⟨\dot{\mathcal{E}}⟩}_{\mathrm{flare}}$, we use the appropriate efficiency, ${\eta }_{\mathrm{flare}}={⟨\dot{\mathcal{E}}⟩}_{\mathrm{flare}}/{⟨\dot{M}⟩}_{\mathrm{cycle}}\sim 0.1$, measured in Section 3.1 but rounded here to the nearest order of magnitude. Finally, we assume that particles radiate efficiently, so that the observed bolometric luminosity of the flare, L, is close to ${⟨\dot{\mathcal{E}}⟩}_{\mathrm{flare}}$. Putting all this together, the flux saturation condition becomes $\varphi ={\mathrm{\Phi }}_{\mathrm{H}}/\sqrt{L/{\eta }_{\mathrm{flare}}}$. Rearranging, and using Φ_H ≃ 2πB_Hr_H², where r_H = 2r_g for a non-spinning black hole, we obtain

  $$\begin{array}{c}\hfill L\sim {10}^{35}\mathrm{erg}{\mathrm{s}}^{-1}\left(\frac{{\eta }_{\mathrm{flare}}}{0.1}\right){\left(\frac{M_{\mathrm{BH}}}{4\times {10}^6{M}_{\odot }}\right)}^2{\left(\frac{B_{\mathrm{H}}}{30\mathrm{G}}\right)}^2{\left(\frac{\varphi }{60}\right)}^{-2}.\end{array}$$(34)

  This luminosity is in the observable range of Sgr A^★ X-ray flares, which occur once per day with a typical luminosity of L ∼ 10³⁴ erg s⁻¹ but in extreme cases can reach almost L ∼ 10³⁶ erg s⁻¹ (Genzel et al. 2010; Ponti et al. 2015). The above estimate would decrease for an imperfect radiative efficiency and probably also for an eruption with finite azimuthal extent – possible only in three dimensions (3D).

• Photon energies: As detailed in Sect. 4.2, we predict the electron energy distribution to extend up to Lorentz factors of γ ∼ σ_Hm_i/m_e ∼ 10 m_i/m_e before cutting off. Using the proton-electron mass ratio, we expect electron Lorentz factors up to at least 10⁴. This translates to a synchrotron photon energy of

  $$\begin{array}{c}\hfill {\epsilon }_c=\frac{3}{2}{\gamma }^2\frac{eħB_{\mathrm{H}}}{m_{e}c}\sim 50\mathrm{eV}\left(\frac{B_{\mathrm{H}}}{30\mathrm{G}}\right){\left(\frac{\gamma }{{10}^4}\right)}^2.\end{array}$$(35)

  Below ∼50 eV, reconnection results in a hard electron power-law index, explaining the rising νF_ν spectra typical of IR Sgr A^★ flares (Ponti et al. 2017; GRAVITY Collaboration 2021). Above ∼50 eV, the electron distribution may or may not continue. Our 2D simulations potentially impose an artificial cutoff here; 3D studies of reconnection in slab geometry suggest that accelerated particles can extend to much higher energies, owing to the uniquely 3D possibility that particles escape plasmoids to undergo additional acceleration (Zhang et al. 2021; Chernoglazov et al. 2023). Unless σ_H reaches values much higher than in our simulations, such acceleration is necessary in order for particles to reach X-ray emitting energies. We should also note that estimates of B_H based on one-zone models may be lower than the actual near-horizon value. Indeed, since our simulations show a steep radial dependence of B^r ∝ r⁻², the true local value of B_H may exceed observation-based inferences by up to an order of magnitude.

• Astrometry and polarization: Our axisymmetric, 2D simulations enforce eruptions to occur simultaneously at all azimuthal angles, ϕ, which is in tension with Sgr A^★ observations. In reality, the Sgr A^★ flares show orbital motion – both in real space and in linear polarization – of a localized hotspot around the black hole (GRAVITY Collaboration 2018), suggesting an azimuthally restricted emitting region. Our model could be refined to capture these features through two main improvements. First, extending our setup to full 3D would allow for non-axisymmetric eruptions. Second, in order to reproduce orbital motion, our models would need to include angular momentum either embedded in the spacetime (in the form of black hole spin) or into the inflowing plasma.

In summary, our simulations reproduce the characteristic duration (∼30 min) and luminosity of Sgr A^★ flares. The electron energies and recurrence times that emerge from our models reproduce near-IR flares, but more theoretical work is needed to more robustly constrain the recurrence time. Future 3D simulations may be able to show how particles are accelerated all the way up to X-ray-emitting energies while enabling study of non-axisymmetric eruptions and the formation of orbiting hotspots.

6.2. Stellar-mass black holes accreting the Galactic ISM

The Milky Way is expected to host up to ∼10⁹ stellar-mass black holes (Shapiro & Teukolsky 1983; Agol & Kamionkowski 2002; Olejak et al. 2020), one of which has recently been detected by gravitational microlensing (Sahu et al. 2022; Lam et al. 2022; Lam & Lu 2023). These objects should accrete their surrounding interstellar matter and magnetic field. We present here a discussion of the characteristic luminosities, variability timescales, and peak photon energies expected from such sources were they to accrete via the accretion-eruption mechanism studied here.

As shown in our simulations, the typical duration of an eruption-powered flare is ∼100t_g = 5 ms (M_BH/10 M_⊙). The recurrence of such flares is not well constrained by our model, as explained in Sect. 6.1, but for the present purposes we speculate that it matches our simulations (and the Sgr A^★ near-IR flares, rescaled to a smaller black hole mass). This gives an estimate for the accretion-eruption quasi-period of ∼1000t_g = 50 ms (M_BH/10 M_⊙).

Following the discussion of Sect. 6.1, the efficiency with which the black hole converts infalling rest-mass energy into plasma energization, averaged over the full accretion-eruption cycle, is ${\eta }_{\mathrm{cycle}}={⟨\dot{\mathcal{E}}⟩}_{\mathrm{cycle}}/{⟨\dot{M}⟩}_{\mathrm{cycle}}\sim 0.03$. If the plasma radiates efficiently, the cycle-averaged luminosity¹ is $L\sim {⟨\dot{\mathcal{E}}⟩}_{\mathrm{cycle}}=0.03{⟨M⟩}_{\mathrm{cycle}}$. Typical accretion rate estimates for isolated stellar-mass black holes feeding on the ISM lie in the range 10^{10 − 15} g s⁻¹ (Agol & Kamionkowski 2002; Matsumoto et al. 2018; Martinez et al. 2025), translating to L ∼ 10^{29 − 33} erg s⁻¹.

Next, we estimate the magnetic field strength at the horizon of a stellar-mass black hole accreting the ISM. Using the flux saturation condition from phase two of our model, $\varphi ={\mathrm{\Phi }}_{\mathrm{H}}/\sqrt{{\dot{M}}_{\mathrm{H}}}\simeq 60$, we estimate that, at the beginning of an eruption,

$$\begin{array}{c}\hfill B_{\mathrm{H}}\sim 3\times {10}^4\mathrm{G}{\left(\frac{M_{\mathrm{BH}}}{10M_{\odot }}\right)}^{-1}{\left(\frac{{\dot{M}}_{\mathrm{H}}}{{10}^{10}\text{g}{\text{s}}^{-1}}\right)}^{1/2}.\end{array}$$(36)

Such magnetic fields are much stronger than in the case of Sgr A^★, amplifying the expected photon energies. Electrons accelerated up to γ ∼ 10⁴, as predicted by our model, radiate synchrotron photons of energy

$$\begin{array}{c}\hfill {\epsilon }_c=\frac{3}{2}{\gamma }^2\frac{eħB_{\mathrm{H}}}{m_{e}c}\sim 20\mathrm{keV}\left(\frac{B_{\mathrm{H}}}{{10}^4\mathrm{G}}\right){\left(\frac{\gamma }{{10}^4}\right)}^2.\end{array}$$(37)

We stress, however (as discussed in Section 6.1), that higher electron energies might be possible in 3D. With a spectrum rising at least until ε_c, Equation (37) suggests that the spectral peak of isolated black holes lies not in the soft X-rays, but in the hard X-rays and, for the most vigorously accreting sources, beyond. Given current instrumental sensitivities, flares peaking in this range might not be easy to detect (see the recent attempt by Mereghetti et al. 2025). However, given the number of potential sources, their overall population might power a Galactic diffuse hard X-ray background.

7. Conclusions

In this work, we have performed self-consistent kinetic simulations that capture a cyclic eruptive activity driven by the accretion of zero-angular-momentum, magnetized plasma onto a Schwarzschild black hole. We show that magnetic reconnection plays a central role in regulating magnetic-flux accretion, triggering magnetospheric eruptions, and accelerating particles to ultrarelativistic energies. Building upon analytical models and empirical knowledge of plasma instabilities, we derive self-consistent constraints on the characteristic particle Lorentz factors and eruption timescales, finding encouraging agreement with observational expectations for Sgr A^★. Because our model is agnostic to the black-hole mass, similar processes may also occur in accreting, isolated stellar-mass black holes, potentially contributing to diffuse galactic high-energy emission.

However, our current setup lacks key features required to reproduce several observational properties of Sgr A^★’s flares, particularly the orbital motion of the emitting region. We expect that incorporating angular momentum – either as black-hole spin or in the inflowing plasma – together with fully 3D simulations, will enable us to capture these effects. Such a more realistic framework will also allow for self-consistent modeling of the associated synchrotron emission, including the production of light curves, polarization maps, and spectra, which we leave for future work.

Acknowledgments

The authors thank I. El Mellah, B. Crinquand, N. Scepi, and K. Kin for insightful discussions. We are thankful for the anonymous referee’s comments, which helped improve some of the discussions in the paper. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 863412). JM is supported by a grant from the Simons Foundation (MP-SCMPS-00001470). AS is supported by NSF through grant AST-2508744. Computing resources were provided by TGCC under the allocations A0170407669 made by GENCI.

References

All Figures

Fig. 1. Sketch of the initial GRPIC setup representing zero-angular-momentum accretion onto a central black hole (labeled BH) in the center. On the left, we show in cyan the initial cloud of plasma, at a density n0 and temperature T0, and in pink the injection area of fresh plasma. On the right, the initial uniform magnetic field lines are shown in blue, and the limit of the outer matching layer in red.

In the text

Fig. 2. Top half: Time evolution of: the horizon magnetic flux, ΦH; the horizon accretion rate, Ṁ, normalized by ${\dot{M}}_0=4\pi m_i{n}_0{r}_{\mathrm{inj}}^2{v}_{\mathrm{th}}$; the normalized horizon magnetic flux, ${\mathrm{\Phi }}_{\mathrm{H}}/\sqrt{\dot{M}}$; and the total dissipation rate, $\dot{\mathcal{E}}$, normalized by Ṁ0. Bottom half: Simulation snapshots at t/tg = 2689, 3286, 3596,  and 3648. Each snapshot shows the number density n (top left), the electromagnetic energy dissipation J ⋅ E (top right), the magnetization σ (bottom left) and the local average particle Lorentz factor ⟨γ − 1⟩ (bottom right). Black lines show magnetic field lines.

In the text

Fig. 3. Top panel: Particle energy distributions at different stages of the eruption cycle. Transparent lines show instantaneous distributions; opaque lines represent distributions time-averaged over each phase. Phases one, two, and three are color-coded, respectively, as green, blue, and red. Bottom panel: Time evolution of ΦH over one eruption cycle. Shaded regions (green, blue, red) indicate the time intervals during which particle distributions are presented in the top panel.

In the text

	Fig. 4. Top: map of the azimuthal electric current carried by the electrons (Jϕ−). Center: map of the azimuthal electric current carried by the ions (Jϕ+). Bottom: map of the auxiliary azimuthal magnetic field Hϕ. Magnetic field lines are represented in solid black lines for all panels. The mass ratio is mi/me = 256 in the depicted simulation.
In the text

	Fig. 5. Map of the local average particle energy for electrons (left) and ions (right) for the simulation with mi/me = 256. The energies are normalized to mec2 and grey lines represent magnetic field lines.
In the text

	Fig. 6. Time-averaged energy distributions of the particles during the erupting phase for all electron-ion simulations (including the reference mi = me run). Ions are represented as dashed lines and electrons as solid lines.
In the text

Fig. 7. Profiles of Vr measured from the simulation presented in Section 3. Solid lines denote time averages over the intervals defined in Fig. 3; envelopes denote one-sigma percentiles for each averaging interval. The orange dashed curved corresponds to Eq. (17) where ${\dot{\mathrm{\Phi }}}_{\mathrm{H}}$ is measured during phase one on Fig. 3. The dashed black line denotes the free-fall velocity of a particle starting from rest at infinity (Eq. (18)).

In the text

	Fig. 8. Traces of δ(r)/(r − rH) at specific times within (transparent lines), and time-averaged over (solid dashed lines), each phase.
In the text

	Fig. 9. Conceptual three-phase accretion cycle. Transitions between phases coincide with critical opening angles of the magnetic field lines about the equator. Blue and green field lines are accreted during phases one and two, respectively; black field lines never reach the black hole.
In the text