© The Authors 2026

1. Introduction

The accretion of matter onto black holes is central to most high-energy astrophysical environments as the engine that powers intense cosmic processes, such as the launching of relativistic jets of magnetized plasma (Blandford & Znajek 1977; Begelman et al. 1984). Because its role in shaping the dynamical behavior of these most extreme of environments in nature is crucial, a significant amount of literature over the span of more than four decades has been devoted to uncovering the details of matter accretion onto black holes.

The accreting material is typically thought to assume the configuration of a thin disk (Novikov & Thorne 1973; Shakura & Sunyaev 1973; Pringle et al. 1973) or of a geometrically thick torus (Fishbone & Moncrief 1976; Abramowicz et al. 1978) around the black hole. As it accretes, its magnetic field, which is frozen into the plasma, is advected toward the black hole, with magnetic flux in the vicinity of the event horizon steadily increasing until the accumulated magnetic field is strong enough to counteract the ram pressure of the infalling gas. Accretion then proceeds intermittently in the highly transient MAD state (Bisnovatyi-Kogan & Ruzmaikin 1974, 1976; Narayan et al. 2003), which has been confirmed by numerical simulations (Igumenshchev 2008) and agrees with the ordered magnetic fields observed in M87* (Event Horizon Telescope Collaboration 2021) and Sgr A* (Event Horizon Telescope Collaboration 2024).

The MAD state, however, appears not to be exclusive to the accretion of a Fishbone-Moncrief (FM) torus onto a spinning black hole. Rather, it has been observed to emerge in simulations that present considerable differences to the idealized FM setup that is traditionally employed in general relativistic magnetohydrodynamic (GRMHD) simulations of accretion flows, provided that sufficient net poloidal magnetic flux accumulates near the black hole over the relevant accretion timescale. Specifically, Ressler et al. (2021) simulated the spherical accretion of zero angular momentum matter with a uniform magnetic field onto spinning black holes. Their setup reached states presenting several similarities to the MAD state achieved by typical FM setups. A long-duration study of a comparable Bondi-like feeding configuration by Lalakos et al. (2024) showed that such a system can develop a more traditional MAD state with powerful jets when sufficient magnetic flux accumulates near the horizon, while subsequently exhibiting additional transitions in its accretion and jet behavior. Generalizing this setup to the spherical accretion of low angular momentum matter, Galishnikova et al. (2025) found that systems of this type also reach states with similarities to the typical MAD state over the simulated duration, but with weaker jets and stronger variability.

Moreover, using a multiscale GRMHD approach, Guo et al. (2025) observed the emergence of the typical MAD state in simulations of accreting FM tori of a much larger scale than typical ones, while Lalakos et al. (2025) used contiguous long-duration GRMHD accretion simulations of an initially uniform medium onto a rotating black hole and observed the emergence of the MAD state for all scale separations considered, determined by the value of the Bondi radius. Taken together, these results indicate that MAD formation is possible beyond the traditional FM torus setup, but that the physical state of the system depends on the feeding conditions and on whether it has evolved for sufficiently long times to reach its asymptotic behavior. In this respect, these studies may correspond to different evolutionary stages of similar systems, while the standard FM torus setup largely bypasses the initial flux-accumulation timescale. The timescale required to reach a time-averaged MAD state may also depend on the angular momentum content of the inflow.

In the MAD state, the excess magnetic flux in the vicinity of the black hole is episodically expelled back into the disk in what are known as flux eruption events (Igumenshchev 2008; Tchekhovskoy et al. 2011). During these events, coherent long-lived bundles of vertical magnetic field are generated close to the equatorial plane (Ripperda et al. 2020, 2022). These outwardly moving flux tubes are filled with low-density highly magnetized plasma, with their footpoints anchored in regions of the equatorial matter distribution where the vertical magnetic field is dynamically dominant (Porth et al. 2021). This whole process repeats multiple times during a typical simulation in the MAD state, resulting in the high temporal variability exhibited by the mass accretion rate and the magnetic flux through the black hole horizon in MAD simulations (Tchekhovskoy et al. 2011; Narayan et al. 2022; Chatterjee & Narayan 2022).

Recent studies have emphasized that hot spots formed via reconnection are central to understanding the flaring activity observed in systems such as Sgr A* (Younsi & Wu 2015; Nathanail et al. 2022; Lin et al. 2023; El Mellah et al. 2023; Dimitropoulos et al. 2025). Antonopoulou et al. (2025) showed that magnetic reconnection during flux eruption events in MADs produces hot spots that travel outward along newly formed magnetized flux tubes, giving rise to near-infrared flares whose properties match the observational characteristics of Sgr A* (GRAVITY Collaboration 2018).

Crucially, the same eruptive processes that drive hot-spot formation also create favorable conditions for particle acceleration, namely the formation of magnetic reconnection sites (Sironi & Spitkovsky 2014; Werner & Uzdensky 2017). Specifically, Vos et al. (2025) found that reconnection during flux eruptions can efficiently accelerate particles of both hadronic and leptonic nature, and trigger enhanced pair production. These energetic particles and their secondary products contribute to the total cosmic-ray reservoir of a galaxy, which shapes the diffuse X-ray, γ-ray, and neutrino emission (Amenomori et al. 2021; Kantzas & Calore 2024) and can affect the physical and chemical properties of the interstellar medium on a galactic scale (Koutsoumpou et al. 2025a,b).

Since reconnection is central to the dynamics of flux eruption events, resistivity of the accreting plasma can possess a critical role in shaping the evolution of these events, and consequently, in the dynamical evolution and variability of MADs. Nathanail et al. (2025) studied the effects of including physical resistivity in GRMHD simulations of MADs, finding that low-resistivity models are, like their ideal MHD counterparts, highly variable because their evolution is determined by frequent reconnection events, while a moderate resistivity can effectively diffuse the accumulated magnetic flux. Moreover, (Aktar et al. 2025) corroborated the results of (Nathanail et al. 2025) and additionally showed that high resistivity can efficiently diffuse the magnetic field of the accreting matter, suppressing plasmoid formation and turbulence related to the magnetorotational instability (MRI), as well as leading to significantly lower jet powers.

We investigated the generation of vertical flux tubes and the evolution of the non-axisymmetric features of the equatorial accretion flow during a flux eruption event by using the results of a standard 3D GRMHD simulation that reached the MAD state. In Sect. 2 we present the setup of our numerical simulation. In Sect. 3 we describe the mechanism via which vertical magnetic flux tubes are detached from the black hole horizon and subsequently transported outward, effectively reducing the accumulated magnetic flux. The description of this mechanism is supported by results from our numerical simulation. Section 4 is dedicated to the examination of the azimuthal structure of the equatorial plasma distribution, featuring quantitative results regarding the evolution of the non-axisymmetric features and dominant azimuthal modes of the equatorial flow, which are intrinsically connected to the formation of the vertical flux tube. We outline our conclusions in Sect. 5.

2. Numerical setup

2.1. Simulation parameters

We used the results of a 3D GRMHD simulation performed with the open-source module of the MPI-AMRVAC framework, BHAC (Porth et al. 2017). BHAC employs second-order shock-capturing finite-volume methods to solve the equations of ideal GRMHD,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mu }(\rho u^{\mu })=0,\end{array}$$(1)

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\nu }{T}^{\mu \nu }=0,\end{array}$$(2)

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mu }{(}^{★}{F}^{\mu \nu })=0,\end{array}$$(3)

where ρ is the fluid rest-mass density, u^μ is its four-velocity, T^μν is the energy-momentum tensor, and ^★F^μν is the Hodge dual of the electromagnetic tensor, F^μν. The code also uses constraint transport methods compatible with adaptive mesh refinement (AMR) (Londrillo & del Zanna 2004; Del Zanna et al. 2007) to satisfy the solenoidal condition, ∇ · B = 0, for the magnetic field (Olivares et al. 2019). BHAC has been widely used for numerical investigations of various relativistic environments and processes (Nathanail et al. 2019; Cruz-Osorio et al. 2022; Mpisketzis et al. 2024; Das et al. 2025) and has undergone extensive benchmarking against similar codes (Porth et al. 2019).

The three-dimensional simulation was performed in modified Kerr-Schild (KS) coordinates (McKinney & Gammie 2004) and featured an effective resolution of N_r × N_θ × N_ϕ = 384 × 192 × 192. As initial conditions, we used a standard axisymmetric Fishbone-Moncrief (FM) torus (Fishbone & Moncrief 1976) with a constant specific angular momentum of l = 6.76, in hydrodynamic equilibrium around a Kerr black hole with a dimensionless spin of α = J/M² = 0.94 and an outer horizon radius of r_H ≃ 1.34 r_g. The inner edge of the initial torus was set at r_in = 20 r_g from the central black hole, while the location of the density maximum was located at r_max = 40 r_g, where r_g = M is the gravitational radius of the central Kerr black hole with mass M, in units in which G = c = 1. Last, we used an ideal equation of state with a constant adiabatic index, $\widehat{\gamma }=4/3$.

The initial magnetic field was purely poloidal and determined by the vector potential,

$$\begin{array}{c}\hfill A_{\varphi }\propto max[\frac{\rho }{{\rho }_{\text{max}}}{\left(\frac{r}{r_{\text{in}}}\right)}^3{\mathrm{e}}^{-r/r_{\text{mag}}}{sin}^3\theta -A_{\varphi ,\text{cut}},0],\end{array}$$(4)

where A_ϕ, cut = 0.2, r_mag = 400 r_g, and ρ_max = 1 is the maximum value of the initial density distribution. This prescription for the initial magnetic vector potential corresponds to a purely poloidal magnetic field in a single nested-loop configuration and is standard for simulations that reach quasi-steady MAD state (Wong et al. 2021; Zhang et al. 2024). Furthermore, the strength of the initial magnetic field was normalized so that the ratio of the maximum gas pressure over the maximum magnetic pressure in the initial FM torus, ${\beta }_{\text{max}}=\frac{{(P_g)}_{\text{max}}}{{(P_{\text{mag}})}_{\text{max}}}=100$. We note that (P_g)_max and (P_mag)_max do not necessarily occur in the same location in the initial torus. Finally, numerical density floors were treated in the standard manner for GRMHD codes (Porth et al. 2017; Nathanail et al. 2020).

2.2. Identification of flux eruption events

We defined the mass accretion rate, Ṁ, as the integral of the mass flux, ρu^r, over the surface of the black hole event horizon, which is located at r ≃ r_H. This integral is

$$\begin{array}{c}\hfill \dot{M}(t)={\displaystyle {\int }_{\varphi }}{\displaystyle {\int }_{\theta }}\rho u^r\mathrm{d}{A}_{\theta \varphi }.\end{array}$$(5)

Similarly, the dimensional or absolute magnetic flux at the event horizon, Φ, is

$$\begin{array}{c}\hfill \mathrm{\Phi }(t)=\frac{\sqrt{4\pi }}{2}{\displaystyle {\int }_{\varphi }}{\int }_{\theta }\left|B^r\right|\mathrm{d}{A}_{\theta \varphi }.\end{array}$$(6)

Following convention, we multiplied the magnitude of the radial magnetic field, |B^r| with $\sqrt{4\pi }$, in order to express it in Gaussian units. Additionally, the factor 1/2 in Eq. (6) appears in order to specify that Φ is calculated over one hemisphere. The dimensionless normalized magnetic flux, ϕ_BH, can be calculated from Ṁ and Φ as

$$\begin{array}{c}\hfill {\varphi }_{\text{BH}}(t)=\frac{\mathrm{\Phi }(t)}{\sqrt{\dot{M}(t)}}.\end{array}$$(7)

In this system of units, the normalized horizon magnetic flux saturation value is ϕ_BH ≃ 50 (Tchekhovskoy et al. 2011). The mass accretion rate, Ṁ, and normalized magnetic flux, ϕ_BH, at the black hole horizon over the duration of our simulation are presented in Fig. 1.

Fig. 1. Time series of the mass accretion rate (top panel) and normalized magnetic flux (bottom panel) on the black hole outer horizon for the entirety of our simulation. The horizontal dashed black line in the plot of ϕBH corresponds to its time average ⟨ϕBH⟩t = 47.

After t ≃ 12 000 t_g, the simulation reached the MAD state, and Ṁ, ϕ_BH display significant temporal variability, typical for accretion in the MAD state. At certain times in the simulation, ϕ_BH reached a value much higher than its saturation value, with the mass accretion rate plummeting at the same time. During these phases, the absolute magnetic flux, Φ, decreased from a local maximum value. This suggests the removal or expulsion of magnetic flux threading the event horizon. Consequently, we considered time intervals over which Φ initially was at a local maximum and then decreased while ϕ_BH surpassed its saturation value as corresponding to flux eruption events.

We identified multiple such events during our simulation and focused our analysis on an event observed over the time interval Δt = 32 191 − 32 360 t_g, with the onset of the event occurring at t = 32 220 t_g, as marked in the time series presented in Fig. 2. We also present results related to three additional events over the analysis intervals Δt = 14 990–15 290 t_g, Δt = 21 890–22 390 t_g, and Δt = 34 190–34 300 t_g in Appendix A. For these three events, the start of the respective analysis intervals corresponds to the onset of the flux eruption. The onset time reported for each event is defined as the time at which the absolute magnetic flux threading the horizon reached a local maximum within the corresponding analysis interval.

Fig. 2. Mass accretion rate, Ṁ, (top panel), dimensional magnetic flux, Φ (middle panel), and normalized magnetic flux, ϕBH, (bottom panel) during the flux eruption event at t = 32 220 tg. The shaded areas of each plot denote the time span of interest.

During the main event, at t = 32 200 t_g, Ṁ and Φ both briefly increased, reaching their respective local maxima at t = 32 207 t_g and t = 32 220 t_g. Both quantities subsequently decreased, while the normalized horizon magnetic flux, ϕ_BH, obtained values higher than its saturation value over the entirety of the time interval of interest. We note that due to the strong and rapid variability of the accretion flow in the MAD state and especially during flux eruption events, the extent of the radial inflow equilibrium region also varied significantly. After the system achieved its quasi-steady state, the radial inflow equilibrium region extended out to r ≥ 40 r_g on average.

3. Formation and motion of vertical flux tubes

Accreting systems in the MAD state are characterized by magnetic flux eruptions. During these episodic events, magnetic flux in the vicinity of the black hole is expelled outward (Porth et al. 2021) and magnetic reconnection forms hot spots in the equatorial accretion disk, as well as vertical magnetic flux tubes that thread the disk (Ripperda et al. 2020, 2022). A deep understanding of the mechanism that triggers flux eruption events and leads to the formation of their characteristic accompanying structures is necessitated by their connection to the observed flaring behavior of accreting environments (Dexter et al. 2020; Hakobyan et al. 2023; Antonopoulou et al. 2025).

3.1. The detachment instability mechanism

Magnetic reconnection plays a central role in the formation of the vertical flux tubes associated with flux eruption events by essentially reconfiguring the equatorial magnetic field on the accretion disk into bundles of vertical magnetic field lines (Ripperda et al. 2020, 2022). These vertical magnetic flux tubes are large-scale structures consisting of coherent vertical magnetic field lines, filled with low-density plasma (Porth et al. 2021; Salas et al. 2024). The generation and subsequent motion of the vertical flux tubes is a physically rich and complex multistage process that we summarize schematically in Fig. 3.

Fig. 3. Schematic depiction of the physical mechanism through which magnetic flux is expelled from the black hole. Magnetic field lines threading the horizon are shown as thin black lines. The bundle of magnetic field lines of interest is depicted as a thick red line. The temporal evolution of the mechanism is shown from left to right: (a) Magnetic flux is advected by the accreting matter toward the black hole. Close to the black hole the magnetic field is mostly equatorial. (b) The equatorial magnetic field slightly above and below the midplane reconnect. (c) Two structures form during the magnetic reconnection of the equatorial magnetic field, a magnetic field loop close to the black hole’s event horizon, and a tube of vertical magnetic field lines. (d) The vertical flux tube, filled with plasma of lower density than its environment, is buoyantly carried away from the black hole.

At the onset of the flux eruption event, the disk magnetic field close to the black hole is dominated by its horizontal component, that is, its component parallel to the midplane. The magnetic field is advected along with the accreting matter, also in a horizontal configuration at small radii (panel a of Fig. 3). As mass slides through the field lines connected to the black hole, the upper parts of the flux tube become magnetically buoyant and push outward, away from the black hole, via a form of the Parker instability (Parker 1955, 1966). This increases the tension of the magnetic field lines. The tension increases further as matter pushes the field lines closer to the black hole, and the field lines eventually become tightly compressed above and below the equator.

At some point, two horizontal magnetic field lines close to the black hole reconnect, one above and one below the midplane, and form an X-point on the equatorial plane, typically at a distance of a few (3–5) r_g from the black hole (panel b of Fig. 3). The reconnection of the two horizontal magnetic field lines forms two structures, a magnetic field loop close to the black hole horizon, and a vertical magnetic flux tube that threads the equatorial plane (panel c of Fig. 3). Part of the reconnecting field is accreted with the plasma, leading to the detachment of the flux tube (panel d of Fig. 3). Neighboring magnetic field lines immediately follow, and the extent of the reconnecting region on the midplane expands. Thus, magnetic field lines are continuously detached from the black hole and added to the flux tube as part of a cascading process, which we term detachment instability. The vertical flux tube, filled with plasma with a lower density than its surrounding environment, is magnetically buoyant and transported outward, away from the black hole. Through this process, the excess magnetic flux concentrated in the vicinity of the horizon is detached from the black hole and is then buoyantly transported away, resulting in the observed decrease in the horizon magnetic flux during flux eruption events.

The detachment instability described in this subsection is a mechanism wholly distinct from the typical MHD instabilities, such as the magnetic Rayleigh-Taylor or the Kelvin-Helmholtz instability that may emerge in a turbulent plasma flow, such as an accretion disk. It is not driven by a tendency of heavy fluid parcels to decrease their potential energy, as in the case of the magnetic Rayleigh-Taylor instability, or by velocity shear between adjacent layers of different characteristics, as in the case of the Kelvin-Helmholtz instability, and neither does it lead to mixing, like these two characteristic MHD instabilities. Instead, the detachment instability is a reconnection-driven process. It arises in the low-density highly magnetized inner disk where horizontal field lines above and below the midplane reconnect, reconfiguring the topology of horizon-threading flux lines into vertical bundles. Subsequently, neighboring field lines follow, producing a cascade of reconnection-induced reconfiguration in which magnetic flux is continuously removed from the black hole and added to the growing flux tube, after which the tube, being strongly magnetized and less dense than its environment, is buoyantly transported outward.

3.2. Numerical results

The physical scenario described in the previous subsection is supported by the results of our standard 3D GRMHD simulation. As explained in the previous section, we identified multiple flux eruption events during our simulation, focusing our analysis on an event starting at t = 32 220 t_g. In order to study the conditions in the equatorial accretion flow during this flux eruption in detail, we obtained outputs with a frequency of 1 t_g. We repeated our analysis for three additional events, one at t = 14 990 t_g, another at t = 21 890 t_g, and the last at t = 34 190 t_g, for results with a lower output frequency of 50 t_g. These results, as well as a convergence check of the results of the flux eruption event of the main text with the output frequency, are presented in Appendix A.

We first examined the distribution and kinematics of the accreting plasma on the equatorial plane by analyzing equatorial slices of the plasma density, ρ, magnetization, σ = b²/ρ, and radial four-velocity, u^r. These quantities indicate the gradual formation of an expanding region of low density, where the plasma is strongly magnetized and moves outward with a strongly positive radial velocity, as shown in Fig. 4.

Fig. 4. Equatorial plasma density, ρ (left), magnetization, σ (middle), and radial four-velocity ur (right) at four different snapshots during the flux eruption event. The first row corresponds to the maximum of Ṁ at t = 32207 tg. An extended region with a high magnetization and low density forms during the event. The low-density plasma in the high-magnetization region has a positive radial four-velocity. The plasma at the footpoint of the flux tube, which begins to form at the edge of the high σ region at negative y, also has a positive radial four-velocity.

The area of the equatorial plane corresponding to the footpoint of the vertical flux tube that is generated during the flux eruption is better identified by considering quantities related to the magnetic field. The first of these quantities is defined as (B^z)²/P_g, with (B^z)²/P_g ≥ 1 in regions where the vertical magnetic field dynamically dominates the matter (Porth et al. 2021). The second quantity related to the magnetic field is σ(B^z)²/P_g. This quantity traces the most energetic part of the flux tube footpoint (Antonopoulou et al. 2025), that is, the region of the midplane of newly formed vertical magnetic field by magnetic reconnection. The third quantity related to the magnetic field is the angle, ψ, between the magnetic field and the equatorial plane, defined as

$$\begin{array}{c}\hfill \psi =arctan\left(\frac{|B^z|}{|B_{\text{hor}}|}\right),\end{array}$$(8)

where B_hor is the horizontal component of the magnetic field. These three quantities are displayed in Fig. 5.

Fig. 5. (Bz)2/Pg (left), σ(Bz)2/Pg (middle), and angle ψ (right) of the magnetic field with the equatorial plane. At t = 32236 tg, a bundle of ordered vertical magnetic field lines has begun to form, with its most energetic part threading the midplane at a distance of about 6 rg from the black hole. The area of the flux tube footpoint expands over time as more equatorial magnetic field lines reconnect and form additional vertical magnetic lines, shown in the middle column.

Initially, only a small region of dominant magnetic field lies close to the black hole, which has formed before the onset of the event and slowly disappears during the event. This region is visible in the first three snapshots of the left column of Fig. 5. In the second snapshot, shown in the left column of Fig. 5, a region where a strong B^z threads the disk has begun to form. This region is the footpoint of the vertical flux tube that is generated during the event. This same region at the same timestamp is also characterized by a high value of σ(B^z)²/P_g, which traces regions of newly reconnected field lines on the midplane. In the next timestamp, after 29 t_g, this region has expanded, as more horizontal field lines reconnect into a vertical configuration. The black region of high σ(B^z)²/P_g in the middle column of Fig. 5 is smaller than the region of strong B^z depicted in the left column of Fig. 5, but it traces the most energetic part of the flux tube footpoint far more accurately. We also note that this region is located at the edge of the low-density and high-magnetization region shown in Fig. 4.

The angle ψ of the magnetic field with the equatorial plane shows that the magnetic field is mostly horizontal in the high-magnetization region, having a strong vertical component at the edge of this region. The configuration of the magnetic field at the edge of the expanding high-magnetization region is highly transient. However, in the large region of ψ ≥ 70 deg located at negative y at late times of our time interval, the dominant vertical magnetic field has achieved a slowly evolving and long-lasting configuration. This region corresponds to the footpoint of the vertical flux tube that formed during the event. This region is filled with plasma with a lower density than that in the surrounding medium. It therefore slowly moves away from the black hole as a result of magnetic buoyancy, as shown by its positive radial four-velocity in Fig. 4. This buoyant outward motion of the vertical flux tube formed during our analysis time interval is similar to the motion of buoyant low-density flux tubes that transport magnetic flux away from the magnetospheric boundary in instances of stellar accretion (Takasao et al. 2018, 2022). Interestingly, this reconnection-based flaring mechanism for the ejection of magnetic flux from the central black hole as a whole has been observed in 3D simulations of accretion onto magnetized protostars (Takasao et al. 2019).

The similarities between accretion onto protostars and the MAD regime of black hole accretion are remarkable. Takasao et al. (2019) performed simulations of accretion onto protostars lacking a magnetosphere of their own. They found that the accreting material leads to the accumulation of magnetic flux onto the protostar, forming a magnetospheric boundary after accretion initiates. Equatorial field lines of opposite polarity reconnect close to the midplane, leading to the development of areas of low density and high magnetization within the disk that are ejected from the accreting system. This mechanism they invoked for the generation of protostellar flares is similar to the detachment instability mechanism, indicating a common regime of accretion onto objects that lack an ab initio magnetospheric boundary, but instead build one after accretion commences through the accumulation of magnetic flux that is dragged in by the accreting plasma. In this common regime of accretion, reconnection at and around the equatorial plane determines the morphology and evolution of the inner accretion flow. This is different to accretion onto objects with a strong magnetosphere of their own, where the magnetic Rayleigh-Taylor instabilities at the magnetospheric boundary determine the inner accretion flow dynamics and morphology (Kulkarni & Romanova 2008, 2009; Blinova et al. 2016; Parfrey & Tchekhovskoy 2024).

4. Azimuthal structure of the equatorial accretion flow

The results presented in the previous section indicate that the distribution of the accreting matter on the equatorial plane changes significantly over the duration of the flux eruption event. Specifically, the equatorial accretion flow becomes increasingly less axisymmetric over the duration of the flux eruption. The non-axisymmetric features that emerge during the eruption event are of major importance to the dynamical evolution of the accretion flow because they enable the infalling plasma to overcome the magnetic barrier close to the black hole (Narayan et al. 2003; Igumenshchev 2008; McKinney et al. 2012; Porth et al. 2021).

4.1. Measuring the accretion disk’s degree of axisymmetry

Numerical simulations of accretion onto black holes have shown that the emergence of magnetically dominated low-density regions and spiral patterns of infalling matter onto the equatorial plane are a feature that always appears during the evolution of accreting systems in the MAD state (Igumenshchev 2008; Tchekhovskoy et al. 2011; Porth et al. 2021; Ripperda et al. 2022). The emergence of these features disrupts the axisymmetry of the equatorial accretion flow and shapes the azimuthal structure of the equatorial matter distribution. In this subsection, we present a simple quantitative measure of the degree of non-axisymmetry exhibited by the equatorial matter distribution, based on the Fourier expansion of the mass density.

The Fourier expansion of the equatorial mass density constitutes a robust method for deciphering the azimuthal structure of the matter distribution on the midplane. A similar method based on Fourier expansions has been used in previous works to uncover the azimuthal structure of standing shock surfaces in simulations of accretion flows (Nagakura & Yamada 2008; Garain & Kim 2023).

We considered the density ρ(t, ϕ)|_r along the circumference of a circle with a radius r on the equatorial plane. At each snapshot in our analysis interval, this density can be expressed as a Fourier series with respect to ϕ (equivalently, a trigonometric series). The coefficients of this series can be calculated as

$$\begin{array}{c}\hfill D_m(t)=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }}{\rho (t,\varphi )|}_r{\mathrm{e}}^{-im\varphi }\mathrm{d}\varphi .\end{array}$$(9)

For m = 0, the above relation returns the average or axisymmetric term in the density Fourier series, with higher order terms giving the coefficients that express the non-axisymmetric corrections to the axisymmetric term. The Fourier series of the density was calculated up to N = 200. This ensured that the series fully converged to the resolution determined by the simulation grid, which is N_ϕ = 192. We employed the standard criterion for GRMHD simulations, which requires the wavelength of the fastest-growing mode of the MRI to be resolved by a minimum of six grid cells in the azimuthal direction. With an azimuthal resolution of N_ϕ = 192, the largest azimuthal mode number that can be reliably resolved is $m_{\text{max}}=\frac{N_{\varphi }}{6}=\frac{192}{6}=32$. Therefore, structures at m > 32 cannot be reliably resolved.

Additionally, we defined the measure of the degree of non-axisymmetry exhibited by the equatorial matter distribution at each timestamp at a specific radius r as

$$\begin{array}{c}\hfill \mathcal{R}(t)=\frac{2{\displaystyle {\sum }_{m=1}^N}{|D_m(t)|}^2}{|D_0{(t)|}^2}.\end{array}$$(10)

ℛ(t) measures the relevant contribution of the non-axisymmetric (${\displaystyle \sum _{m=1}^N}{|D_m(t)|}^2$) and axisymmetric (|D₀(t)|²) components to the Fourier series and is therefore a suitable metric for the degree of non-axisymmetry displayed by the equatorial matter distribution. In Fig. 6 we present ℛ, calculated at every snapshot of our time interval and at radii r = {r_H, 5 r_H, 10 r_H, 20 r_H}, which shows the transition of the accretion flow morphology from almost axisymmetric to highly non-axisymmetric.

Fig. 6. Measure of the equatorial matter distribution degree of non-axisymmetry ℛ at four different radial positions over the time interval of interest.

At all radii that we considered in our analysis, the equatorial matter distribution displayed a strengthening of its non-axisymmetric features, with this effect being more prominent close to the black hole, at r ≃ r_H and r = 5 r_H. At larger distances from the black hole, the amplification of the non-axisymmetric features of the inner accretion flow was considerably less significant. Moreover, ℛ displayed stronger variability closer to the black hole, and it thus revealed the transient nature of the structures that form in the innermost regions of the accretion disk. The enhancement of large-scale non-axisymmetric features during the flux eruption event is much more significant in the inner parts of the equatorial accretion flow.

We also examined the relative contribution of each of the first four non-axisymmetric terms (m = 1 − 4) in the expansion of ρ through the ratios

$$\begin{array}{c}\hfill {\mathcal{C}}_m(t)=\frac{|D_m{(t)|}^2}{max(|D_m{(t)|}^2)},\end{array}$$(11)

presented in Fig. 7. Along the r ≃ r_H circumference, the initial increase in ℛ is mainly a result of the amplification of the m = 2 terms of the Fourier series, while its further steady increase at later times in our time interval can mainly be attributed to the m = 1 term. At r = 5 r_H, the main non-axisymmetric features in the middle of the time interval appear due to the m = 1, 2 terms of the series. At r = 10 r_H, the dominant non-axisymmetric contribution appears to be due to the m = 3 term, while at r = 20 r_H, the m = 2 azimuthal mode dominates over the first half of the time interval, before being surpassed by the m = 1 term. We note that the strong increase in ℛ at r = r_H and 5 r_H at the end of the time interval is a result of the amplification of the m = 1 term.

Fig. 7. Relative contribution of the first four non-axisymmetric terms at all four radii.

4.2. Dominant azimuthal modes

We subsequently used the method based on the Fourier expansion of the equatorial plasma density introduced in Sect. 2 to determine the dominant azimuthal modes of fundamental quantities related to the equatorial matter distribution. Namely, in addition to the plasma density, we calculated the Fourier expansion coefficients of the mass accretion rate per solid angle, defined as

$$\begin{array}{c}\hfill f_M(t,\varphi ){{|}_r=\rho u^r\sqrt{-g}|}_{r,\theta =\pi /2}.\end{array}$$(12)

Subsequently, we calculated the integrals of the powers of each term in the expansions,

$$\begin{array}{c}\hfill {\mathcal{I}}_m={\displaystyle \int {|D_m(t)|}^2}\mathrm{d}t,\end{array}$$(13)

over the entirety of the time interval. We performed these calculations for ρ and f_M at r = {r_H,  5 r_H, 10 r_H,  20 r_H}. The results are shown in Fig. 8.

Fig. 8. Time-integrated squared Fourier coefficients that express the power in each term of the expansion for the equatorial plasma density, ρ(t, ϕ), and mass accretion rate per solid angle, fM(t, ϕ), integrated over the entire time interval from t = 32 191 tg to t = 32 360 tg. For the main flux-eruption event analyzed in the text, the integral corresponding to the m = 2 azimuthal mode is dominant at all radii except r = 10 rH, where the m = 3 mode dominates.

Fig. 8 shows that the dominant azimuthal mode of the expansions of ρ(t, ϕ)|_rH and f_M(t, ϕ)|_rH for the main event is the m = 2 mode, with the m = 1 being a close second. The same behavior is observed at r = 5 r_H, where the m = 2 mode is again dominant. However, the relative contribution of the m = 1 mode is lower than r = r_H. Farther away from the black hole, at r = 10 r_H, the m = 3 mode dominates the non-axisymmetric spectrum, before the m = 2 again becomes the dominant mode at r = 20 r_H. It is interesting to note that our results for the main event indicate the dominance of low azimuthal number modes on the equatorial plane, mainly of the m = 1 and m = 2 modes close to the black hole. Repetition of this analysis for the three additional events at t = 14 990 t_g, t = 21 890 t_g, and t = 34 190 t_g yielded the same qualitative picture: the equatorial inner accretion flow is dominated by low-order azimuthal modes, although the specific mode with the largest time-integrated Fourier power varies between events and radii (see Appendix A). This indicates that the prominence of low-m modes, especially m = 1 and m = 2, is a robust feature of the equatorial inner accretion flow during flux eruptions in the MAD regime. We note that a highly similar picture is observed in three-dimensional simulations of non-relativistic accretion onto magnetized stars, where accreting matter pierces its way through the preexisting magnetospheric boundary of the central object through the development of low m non-axisymmetric modes (Kulkarni & Romanova 2008, 2009; Blinova et al. 2016; Takasao et al. 2022; Zhu et al. 2024), even though in the MAD scenario, the emergence of low-m non-axisymmetric modes is induced by reconnection, as in protostellar accretion (Takasao et al. 2019), and not driven by MHD instability, as in the aforementioned works.

We must also note that the angular momentum of the initial torus can affect the dynamical evolution and morphology of the accreting system significantly during a flux eruption event. Specifically, in simulations with low angular momentum tori, magnetic flux tubes generated during flux eruption events are able to freely and quickly escape from the inner disk due to magnetic buoyancy. On the other hand, strong azimuthal velocity shear in high angular momentum tori stretches the footpoints of magnetic flux tubes within the disks and leads to a stronger mixing of the highly magnetized plasma within flux tubes with their weakly magnetized environment due to velocity shear, as recently shown by Chan et al. (2025). The angular momentum of the initial torus therefore is a parameter that can affect which azimuthal modes dominate the disk structure.

Moreover, McKinney et al. (2012) reported a dominance of the m = 1 mode close to the black hole in the density Fourier expansion. This difference with respect to our results may arise from a combination of methodological and physical factors. On the methodological side, their Fourier analysis was performed on the density averaged over the full θ extent of the disk, whereas our analysis was based on equatorial slices, which probe a different aspect of the flow morphology. On the physical side, differences in the initial torus setup and in the resulting flow dynamics, including the rotation profile and associated shear, might also affect which azimuthal modes become most prominent. This discrepancy can therefore be attributed to a combination of factors related to the analysis of simulation results and its physical setup.

The inner radii of an accretion flow in the MAD state constitute a region where interchange instabilities can develop and be intensified during energetic phases, such as the flux eruption events studied in this work. First, the interface between the highly magnetized low-density region and its environment, that is, the boundary of the flux tube footpoint, are areas where interchange instabilities develop. The differential rotation and local velocity shear of the inner accretion disk is conducive to the emergence of the Kelvin-Helmholtz instability, which affects the outer surface of the flux tube (Porth et al. 2021; Ripperda et al. 2022). Additionally, the inner parts of the accretion disk in the MAD state are strongly sub-Keplerian (Porth et al. 2021; Dhruv et al. 2025). Consequently, the plasma there experiences a negative effective gravity, which is a favorable condition for the emergence of the magnetic Rayleigh-Taylor instability (Chandrasekhar 1961).

While these types of interchange instabilities develop in MADs, they most probably affect the disk at smaller scales, with the structures that form due to them possessing smaller angular extents and surviving for shorter timescales. Such smaller-angular-scale structures correspond to higher-m modes, so their presence is most naturally interpreted as a signature of local fragmentation, shearing, and interchange activity in the inner flow and not a signature of the large-scale morphology imposed by the reconnection-driven detachment instability.

5. Conclusions

We studied the configuration and dynamical evolution of the equatorial inner accretion flow around a rapidly spinning black hole during magnetic flux eruption events. We used the results of a 3D GRMHD simulation over a time interval that encompassed the duration of a flux eruption, with an output frequency of 1 t_g, in order to obtain a high temporal resolution of the transient dynamics of the inner accretion flow. We repeated our analysis for three additional flux eruption events, using results with a lower output frequency of 50 t_g. Our key findings are summarized in the following paragraphs.

During the flux eruption, a bundle of vertical magnetic flux forms, threading the inner equatorial accretion flow. The formation of this structure is initiated by magnetic reconnection of horizontally oriented magnetic field lines that were advected inward by the accretion flow. As the magnetic field near the black hole becomes increasingly compressed and dynamically dominant, oppositely directed horizontal field lines above and below the equatorial plane reconnect, generating an X-point configuration close to the midplane. This reconnection event detaches a portion of the magnetic field from the horizon, reorganizing it into a vertical configuration threading the disk at a few gravitational radii. The resulting vertical flux tube, filled with plasma with a lower density than its surroundings, becomes buoyantly unstable and rises outward, carrying magnetic flux away from the black hole. This mechanism enables the removal of excess magnetic flux from the horizon region and plays a fundamental role in regulating the saturation of the MAD state.

The flux eruption event is accompanied by the amplification of non-axisymmetric features in the equatorial accretion flow. As the magnetic field accumulates and reconnects, magnetically dominated low-density regions and dense accretion streams emerge on the midplane, breaking the initial axisymmetry of the disk. We quantified the evolution of the azimuthal structure by performing a Fourier decomposition of the equatorial mass density, computing the degree of non-axisymmetry at various radii. Our results showed that the non-axisymmetric features are significantly enhanced closer to the black hole, particularly at r ≃ r_H and r = 5 r_H, while the outer disk remains comparatively more axisymmetric. The evolution of the non-axisymmetric features highlights the highly dynamic and transient nature of the inner accretion flow during magnetic flux eruptions in the MAD state.

We further analyzed the azimuthal structure of the equatorial accretion flow by determining the dominant Fourier modes of key physical quantities. In addition to the plasma density, we computed the Fourier expansion coefficients of the mass accretion rate per solid angle, f_M, and integrated the power of each mode over the duration of the flux eruption event. For the main high-cadence event analyzed in the text, our results showed that the equatorial distributions of ρ and f_M are dominated by low azimuthal number modes. Specifically, the mode with the largest time-integrated Fourier power is the m = 2 mode at r_H,  5 r_H, and 20 r_H, while at 10 r_H the m = 3 mode dominates. These results were obtained by analyzing the structure of the equatorial accretion flow using simulation data that were output with a remarkably high frequency of 1 t_g, with this time interval corresponding to the duration of the flux eruption at 32 220 t_g. The repetition of the same analysis for three additional flux eruptions, at t = 14 990 t_g, t = 21 890 t_g, and t = 34 190 t_g, yielded the same qualitative conclusion: the equatorial inner accretion flow is governed by low-order azimuthal structure, although the specific dominant mode varies between events and radii. Overall, the additional events support the robustness of low-m non-axisymmetric structure during flux eruptions, with m = 1 and m = 2 being the most frequently recurrent dominant modes. The data corresponding to these events were output with a lower frequency of 50 t_g, which was selected after confirming the convergence of results for the event at 32 220 t_g for this lower output frequency, as shown in Fig. A.1 in Appendix A. Overall, we found that an output frequency as low as 50 t_g is adequate for performing a quantitative analysis of an equatorial accretion flow in the MAD state that accurately captures the variable inner disk structure and morphological characteristics of large angular scale such as the one performed in this work.

The dominance of large-scale low-m structures, especially in the vicinity of the black hole (r ≤ 5 r_H), suggests that the overall morphology of the equatorial accretion flow during the flux eruption is controlled by the non-axisymmetric features induced by the reconnection-driven detachment instability mechanism and not by small-scale turbulence. A similar prevalence of low-m modes has been reported in simulations of non-relativistic accretion onto magnetized stars, where accretion streams develop through large-scale instabilities at the magnetospheric boundary.

Although interchange and Kelvin-Helmholtz instabilities are likely active in the highly magnetized inner disk, their effects are confined to smaller spatial and temporal scales. The large-scale azimuthal structure observed in our analysis reflects the reorganization of the inner accretion flow due to the detachment instability mechanism that triggers the reconnection-induced formation and buoyant rise of vertical magnetic flux tubes during flux eruption events. This mechanism is remarkably similar to the mechanism for the production of flares during the Newtonian accretion of matter onto protostars proposed by Takasao et al. (2019). The strong similarities between protostellar accretion and the MAD state of accretion onto black holes suggest the existence of a common regime of accretion onto objects that lack their own magnetospheres, but build up variable magnetospheric boundaries by accumulating magnetic flux by the magnetized matter they accrete.

Acknowledgments

We thank the anonymous referee for their comments and suggestions. AL acknowledges financial support by the State Scholarships Foundation (IKY) scholarship program from the proceeds of the “Nic. D. Chrysovergis” bequest. AN has been supported by the Hellenic Foundation for Research and Innovation (ELIDEK) under Grant No 23698. This work was supported by computational time granted from the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility – ARIS – under project ID 16033. This work was also supported by the Sectoral Development Program (OΠΣ 5223471) of the Ministry of Education, Religious Affairs and Sports, through the National Development Program (NDP) 2021-25.

References

Appendix A: Convergence of results with output frequency

For the analysis of the azimuthal structure of the accretion disk during flux eruption events, we chose an output frequency of 1 t_g for the simulation results inside the time interval corresponding to the analyzed event. This choice was made in order to obtain a high temporal resolution of the equatorial accretion flow that allowed us to examine its morphological evolution as well as the formation of vertical flux bundles in great detail. However, the computational expense of such a high output cadence makes it difficult to reproduce this analysis for many events during the simulation. For this reason, we checked the convergence of the results pertinent to the azimuthal modes of the morphological features of the disk for lower frequency outputs.

Figure A.1 displays the normalized integrals ℐ_m of the Fourier series coefficients of ρ and f_M for an output frequency of 1 t_g as shown in Fig. 8, and for lower output frequencies of 25 t_g and 50 t_g. Results obtained by the lower output frequency data do not present significant differences to the ones obtained by the high 1 t_g output frequency analysis, with the qualitative nature of the results remaining unchanged. We observe that in all three cases, the dominant azimuthal mode at all radii is the same, with the exception of the dominant m of f_M at r = r_H, which can be explained by the very rapid variability of the accretion flow’s morphology very close to the black hole. This convergence check indicates that results regarding the azimuthal structure of the inner accretion flow obtained through data output at a significantly lower frequency are exceedingly similar to the ones obtained by our main analysis of the 1 t_g output frequency results. Therefore, analyses built on low output frequency simulation data, as low as 50 t_g, are equally robust.

Fig. A.1. Time integrated squared Fourier coefficients, ℐm, for the equatorial plasma density, ρ(t, ϕ) and mass accretion rate per solid angle, fM(t, ϕ), integrated over the entire time interval from t = 32191 tg to t = 32360 tg. Left column: 1 tg output frequency, the same as Fig. 8. Middle column: 25 tg output frequency. Right column: 50 tg output cadence. Results do not vary significantly with output frequency.

We repeated our analysis of the accretion flow’s structure during the flux eruption we focused on in the main text for three additional flares, at t = 14990 t_g and t = 21890 t_g, before the event of the main text, and one at t = 34190 t_g, after the event analyzed in the main text. The results are presented in Fig. A.2. The analysis of the three additional flares supports the main conclusion that the azimuthal structure of the equatorial inner accretion flow during flux eruptions is dominated by low-order modes. Although the mode with the largest time-integrated Fourier power is not identical at all radii for all events, the additional events consistently show that m = 1 and m = 2 are the most common dominant modes, while m = 3 becomes dominant only in a more limited subset of cases. Some of the differences between events are likely enhanced by the lower output frequency, which cannot capture the most rapid variability of the inner accretion flow equally well. Nevertheless, the qualitative agreement between the different events supports the robustness of the low-m structure reported in the main text.

Fig. A.2. Time integrated squared Fourier coefficients, which express the power in each term of the expansion, ℐm, for the equatorial plasma density, ρ(t, ϕ) and mass accretion rate per solid angle, fM(t, ϕ), for three additional flux eruption events. Left column: flux eruption event at t = 14990 tg. Middle column: flux eruption event at t = 21890 tg Right column: flux eruption event at 34190 tg.

All Figures

	Fig. 1. Time series of the mass accretion rate (top panel) and normalized magnetic flux (bottom panel) on the black hole outer horizon for the entirety of our simulation. The horizontal dashed black line in the plot of ϕBH corresponds to its time average ⟨ϕBH⟩t = 47.
In the text

	Fig. 2. Mass accretion rate, Ṁ, (top panel), dimensional magnetic flux, Φ (middle panel), and normalized magnetic flux, ϕBH, (bottom panel) during the flux eruption event at t = 32 220 tg. The shaded areas of each plot denote the time span of interest.
In the text

Fig. 3. Schematic depiction of the physical mechanism through which magnetic flux is expelled from the black hole. Magnetic field lines threading the horizon are shown as thin black lines. The bundle of magnetic field lines of interest is depicted as a thick red line. The temporal evolution of the mechanism is shown from left to right: (a) Magnetic flux is advected by the accreting matter toward the black hole. Close to the black hole the magnetic field is mostly equatorial. (b) The equatorial magnetic field slightly above and below the midplane reconnect. (c) Two structures form during the magnetic reconnection of the equatorial magnetic field, a magnetic field loop close to the black hole’s event horizon, and a tube of vertical magnetic field lines. (d) The vertical flux tube, filled with plasma of lower density than its environment, is buoyantly carried away from the black hole.

In the text

Fig. 4. Equatorial plasma density, ρ (left), magnetization, σ (middle), and radial four-velocity ur (right) at four different snapshots during the flux eruption event. The first row corresponds to the maximum of Ṁ at t = 32207 tg. An extended region with a high magnetization and low density forms during the event. The low-density plasma in the high-magnetization region has a positive radial four-velocity. The plasma at the footpoint of the flux tube, which begins to form at the edge of the high σ region at negative y, also has a positive radial four-velocity.

In the text

Fig. 5. (Bz)2/Pg (left), σ(Bz)2/Pg (middle), and angle ψ (right) of the magnetic field with the equatorial plane. At t = 32236 tg, a bundle of ordered vertical magnetic field lines has begun to form, with its most energetic part threading the midplane at a distance of about 6 rg from the black hole. The area of the flux tube footpoint expands over time as more equatorial magnetic field lines reconnect and form additional vertical magnetic lines, shown in the middle column.

In the text

	Fig. 6. Measure of the equatorial matter distribution degree of non-axisymmetry ℛ at four different radial positions over the time interval of interest.
In the text

	Fig. 7. Relative contribution of the first four non-axisymmetric terms at all four radii.
In the text

Fig. 8. Time-integrated squared Fourier coefficients that express the power in each term of the expansion for the equatorial plasma density, ρ(t, ϕ), and mass accretion rate per solid angle, fM(t, ϕ), integrated over the entire time interval from t = 32 191 tg to t = 32 360 tg. For the main flux-eruption event analyzed in the text, the integral corresponding to the m = 2 azimuthal mode is dominant at all radii except r = 10 rH, where the m = 3 mode dominates.

In the text

Fig. A.1. Time integrated squared Fourier coefficients, ℐm, for the equatorial plasma density, ρ(t, ϕ) and mass accretion rate per solid angle, fM(t, ϕ), integrated over the entire time interval from t = 32191 tg to t = 32360 tg. Left column: 1 tg output frequency, the same as Fig. 8. Middle column: 25 tg output frequency. Right column: 50 tg output cadence. Results do not vary significantly with output frequency.

In the text

	Fig. A.2. Time integrated squared Fourier coefficients, which express the power in each term of the expansion, ℐm, for the equatorial plasma density, ρ(t, ϕ) and mass accretion rate per solid angle, fM(t, ϕ), for three additional flux eruption events. Left column: flux eruption event at t = 14990 tg. Middle column: flux eruption event at t = 21890 tg Right column: flux eruption event at 34190 tg.
In the text