© The Authors 2025

1. Introduction

Soon after the advent of the standard theory of accretion disks (Shakura & Sunyaev 1973; Novikov & Thorne 1973), it was recognized (Bisnovatyi-Kogan & Blinnikov 1976; Liang & Price 1977) that a classic geometrically thin, optically thick accretion disk cannot supply the hard X-ray radiation (up to 100 keV or so) observed in black hole (BH) X-ray binaries (XRBs). This led to the postulation that such disks may exist in close proximity to a much hotter and more tenuous plasma. The putative plasma was dubbed, by analogy with the hot and tenuous X-ray emitting plasma enshrouding the Sun, an accretion disk corona (ADC; Bisnovatyi-Kogan & Blinnikov 1976; Liang & Price 1977).

Pursuing the solar analogy farther, Galeev et al. (1979) suggested that convective instabilities in the disk generate magnetic fields that subsequently buoyantly rise up out of the disk, leading to an arcade of ultramagnetized magnetic loops tethered to, and manipulated by, the heavy accretion disk below (see also Tout & Pringle 1992). In this scenario, magnetic fields are not only the hallmark of the coronal morphology, but also constitute the main vehicle through which free energy from the underlying disk rotation can be transformed into hard X-ray emission. For example, relativistic magnetic reconnection (Blackman & Field 1994; Lyutikov & Uzdensky 2003; Lyubarsky 2005) can efficiently transfer magnetic energy to the plasma, offsetting the extremely rapid radiative losses (Galeev et al. 1979; Di Matteo 1998; Liu et al. 2002; de Gouveia dal Pino & Lazarian 2005; Uzdensky & Goodman 2008; Goodman & Uzdensky 2008; Singh et al. 2015; Kadowaki et al. 2015; Khiali et al. 2015; Uzdensky 2016; Beloborodov 2017; though magnetized turbulence may also drive particle energization; Singh et al. 2015; Kadowaki et al. 2015; Grošelj et al. 2024; Nättilä 2024).

Besides merely dissipating magnetic energy, reconnection also reconfigures the magnetic field lines, thereby modifying how they transport energy and angular momentum through the corona. Powering high-energy radiation is thus not the only important role reconnection may play in ADCs. It may also profoundly impact global aspects of accretion, including concomitant ejection phenomena.

For example, reconnection may compete with Keplerian shear and random magnetic footpoint motion to set the coronal scale height and regulate the balance between closed and open field-line bundles (Tout & Pringle 1996; Uzdensky & Goodman 2008). By shaping the closed portion of the corona, reconnection sets the efficiency with which angular momentum is transported from inner to outer closed-loop footpoints, directly impacting accretion (Heyvaerts & Priest 1989; Goodman 2003). At the same time, the influence of reconnection on open field lines is equally important, since these can channel a magneto-centrifugal wind (Blandford & Payne 1982) and thus directly shed angular momentum off the disk (Konigl 1989).

Furthermore, in the event that a coronal magnetic link is established between the accretion disk and a spinning central BH, magnetic reconnection can throttle the magnetic flux threading the event horizon, thereby controlling the extraction of BH spin energy in the form of a Blandford & Znajek (1977) jet. If such a coupling is repeated, with multiple opposing-polarity coronal loops successively fed to the BH and ejected via the Blandford-Znajek (BZ) mechanism (Parfrey et al. 2015), the nascent jet is embedded with a striped magnetic configuration that may dissipate much farther downstream (Giannios & Uzdensky 2019). Clearly, then, by shaping ADCs, reconnection is intricately tied to the more general phenomena of BH accretion, outflow, energy extraction, and jet dissipation.

These mainly theoretical remarks are underscored by observations, especially those of outbursting BH XRB state transitions. In a summary of the main observed features of such outbursts, Fender et al. (2004) pointed out that the binary’s transition to a softer X-ray spectrum and reduced X-ray variability is generally hailed by one or more discrete radio-jet ejections (see also Fender et al. 2009; Miller-Jones et al. 2012; Russell et al. 2019 as well as the reviews by Remillard & McClintock 2006; Done et al. 2007; Ingram & Motta 2019). This bolsters the expected picture in which reconnection-mediated coronal structure, reconnection-powered X-ray flaring, and jet launch are all connected, as indeed has already been appreciated in the literature. Tagger et al. (2004) and de Gouveia dal Pino & Lazarian (2005), for example, posited that the radio-jet ejections are created by reconnection in the inner corona along field lines coupling the accretion disk to the central BH.

Modeling ADCs requires, as a first step, choosing a suitable mathematical description of the coronal plasma. Here, the framework of force-free electrodynamics is often employed, since, barring thin reconnecting current sheets, the corona is thought to be well magnetized so that the plasma inertia and the bulk electromagnetic force both approximately vanish. As a result, scenarios in which coronal activity above razor-thin accretion disks links to jet launch have chiefly been studied in the force-free framework (Parfrey et al. 2015; Yuan et al. 2019; Mahlmann et al. 2020). These works have shown that magnetically coupling a spinning BH to its accretion disk can, for the right coronal loop structure, permit strong BZ jet production. The basic constraint is that the loops must be much larger than the BH. Otherwise, rotational shear is insufficient to open the field lines threading the BH horizon (Uzdensky 2005; Parfrey et al. 2015). At the same time, if the coronal loops encircling the inner one – i.e., the one that couples the disk to the BH – are too large compared to the inner loop, they can squeeze the otherwise jet-like outflow, disrupting it through kink instabilities (Yuan et al. 2019; Mahlmann et al. 2020).

While the force-free approximation works well to describe the corona, at least outside of reconnecting current sheets, it breaks down inside the heavy, geometrically thin accretion disk, where gravitational and pressure forces become important. This is why, in the force-free studies of Parfrey et al. (2015), Yuan et al. (2019), and Mahlmann et al. (2020), accretion is prescribed as a simplified boundary condition rather than as a dynamical disk. The coronal magnetic loops are simply attached to the equator by hand, and inward footpoint dragging is imposed as a proxy for accretion.

To model the disk properly, one must use a more general plasma framework, such as magnetohydrodynamics (MHD), that can capture the non-negligible disk-plasma inertia. However, even for MHD simulations, realizing the razor-thin disk limit is numerically challenging since it imposes prohibitive resolution requirements. Nevertheless, recent works have made great strides toward the thin-disk regime. One recent result is that thin disks can efficiently transport magnetic flux inward in radius (Liska et al. 2019; Jacquemin-Ide et al. 2021; Scepi et al. 2024), which defies preceding theoretical expectations (Lubow et al. 1994). Moreover, even though an efficient disk dynamo has not been observed in the thin-disk limit (Musoke et al. 2023), state-of-the-art simulations of thicker disks (Liska et al. 2020) exhibit the generation of large-scale magnetic loops that are accreted inward and fed to the BH (Jacquemin-Ide et al. 2024). Even in the thick-disk regime, when such large loops accrete onto the BH, they result in efficient BZ jet activation (Chashkina et al. 2021). Taken together, MHD studies appear to validate, at least in part, the simple prescription of inward magnetic loop dragging adopted by force-free models (Parfrey et al. 2015; Yuan et al. 2019; Mahlmann et al. 2020). They further suggest that, even if the disk is, in reality, thicker, and the magnetization contrast less pronounced between disk and corona, the basic looped magnetic structure falling onto the BH may still be preserved (Chashkina et al. 2021; Jacquemin-Ide et al. 2024).

We therefore assume, in this work, that the basic coronal magnetic configuration and dynamics can be studied simply by imposing magnetic footpoint motion as an equatorial boundary condition. We attempt to model this situation in a way that retains one major element that, in the context of a collisionless coronal plasma, neither force-free electrodynamics nor MHD can capture: self-consistent energetics. In particular, we aim for the electromagnetic energy dissipated through collisionless reconnection, as well as the reconnection rate itself, to be governed by first-principles plasma physics as opposed, as in MHD or force-free models, to numerical prescriptions¹. Such self-consistency is absolutely necessary to capture reconnection-powered particle acceleration and radiation, and also to ascertain the role of reconnection in regulating the global coronal structure and BZ jet formation. To achieve our targeted self-consistency, we perform simulations based on the particle-in-cell (PIC) technique (Birdsall & Langdon 1991).

Though needed for our purposes, the PIC method is not perfect. It generally incurs a high computational cost, which limits the separation between plasma microscales and global scales as well as the maximum spatiotemporal extent of simulations. We comment on these limitations as they arise in the body of this work.

Nevertheless, because of its attractive self-consistency property, the PIC method has already been used to simulate local models of reconnection in regimes relevant to ADCs (Werner et al. 2019; Sironi & Beloborodov 2020; Mehlhaff et al. 2021, 2024; Sridhar et al. 2021, 2023; Gupta et al. 2024). However, we are here keenly interested in the impact of reconnection on the global magnetic configuration, which requires us to perform global simulations. Global PIC models featuring a magnetic coupling between a central spinning BH and its accretion disk have recently been presented by El Mellah et al. (2022, 2023). However, their work did not include advection of multiple loops onto the BH. We aim to fold in this fundamental aspect, and so to unveil how reconnection regulates particle acceleration, coronal structure, and jet formation throughout the loop-accretion cycle. In this first paper, we focus on the corona-jet interaction in the regime (uncovered by previous force-free works; Uzdensky 2005; Parfrey et al. 2015; Yuan et al. 2019; Mahlmann et al. 2020) where the loop size is much larger than the BH, leading to efficient jet launch; in a companion publication (Crinquand et al., in prep.), we tackle the complementary regime where the loops are smaller, and the main reconnection dynamics come from the corona interacting with itself.

This paper is structured as follows. In Section 2, we describe our simulation setup, which is heavily inspired by the aforementioned force-free models (Parfrey et al. 2015; Yuan et al. 2019; Mahlmann et al. 2020). In Section 3, we present the results from our simulation, focusing on those aspects for which the PIC approach is uniquely adapted: the global energetics and (in consequence) the radiative signatures. We discuss implications of our findings for ADCs in the context of XRBs, changing-look active galactic nuclei (AGNs), Sagittarius A*, and the striped jet model of Giannios & Uzdensky (2019) in Section 4. We summarize our main findings in Section 5.

2. Setup

We present axisymmetric pair-plasma simulation using the general relativistic PIC (GRPIC) code GRZELTRON (Parfrey et al. 2019). Our simulation setup, sketched in Fig. 1, closely resembles that used by Parfrey et al. (2015) for their force-free modeling. Here, we summarize the main qualitative features of our setup; we detail technical aspects in the subsections that follow.

Fig. 1. Initial conditions (IC) and boundary conditions (BC) of our simulation setup. Sketch is not to scale.

We simulated the upper-half space near a Kerr BH. To model the plasma corona, we attached a tightly packed train of ultramagnetized, purely poloidal (i.e., with no initial azimuthal component), equally sized magnetic loops to the simulation equatorial boundary. At the simulation outset, we began dragging the loop footpoints azimuthally, at the local Keplerian velocity, as well as slowly radially inward. The Keplerian motion is prograde: aligned with the BH spin axis. Rotational shear across loop footpoints builds up an azimuthal component in the magnetic field, magnetically pressurizing the loops and, hence, causing them to inflate toward infinity. Simultaneously, inward field-line dragging feeds magnetic flux to the BH, enforcing a coupling with the corona and allowing a BZ jet to be launched.

We chose a rather extreme dimensionless Kerr spin value of a = 0.99. While we suspect that our results would hold for a broad range of (prograde) spins, the near-extreme case provides a number of numerical conveniences. For example, it extremizes the rotation of field lines threading the BH horizon and thus renders the closed field-line region smaller. This makes it easier to fit all of the regions of distinct magnetic topology into our simulation box (El Mellah et al. 2022). In addition, a near-extreme spin tightens the radius, r_ISCO, of the innermost stable circular orbit. This limits the zone between the event horizon and r_ISCO where our simplified accretion prescription is least realistic.

Our equatorial disk boundary resembles a classic, optically thick, geometrically thin accretion disk (Shakura & Sunyaev 1973). By imposing the disk motion by hand, we ignore both the internal disk physics and the potential backreaction from the corona. However, we expect a wide separation of timescales between the slow thin-disk evolution and the comparatively very rapid coronal dynamics. In this sense, our setup can be viewed as a means of probing the coronal impulse response to a single loop size. Our radial footpoint accretion simply scans across the possible magnetic footpoint configurations for this one loop size.

We employed loops much larger than the central BH. While one does not typically expect large loops to be produced in situ by a thin disk, recent GRMHD results suggest that such loops could be dragged inward from larger radii (Jacquemin-Ide et al. 2024). Alternatively, they could be built up by repeated mergers of smaller loops (Tout & Pringle 1996; Uzdensky & Goodman 2008).

2.1. The 3 + 1 formalism

Our simulations adopt a 3 + 1 formalism, in which the general relativistic line element is expressed as

$$\begin{array}{cc}\hfill \mathrm{d}{s}^2& =g_{\mu \nu }\mathrm{d}{x}^{\mu }\mathrm{d}{x}^{\nu }\hfill \\ \hfill & =({\beta }_l{\beta }^l-{\alpha }^2)c^2\mathrm{d}{t}^2+2{\beta }_i\mathrm{d}{x}^{i}c\mathrm{d}t+h_{\mathit{ij}}\mathrm{d}{x}^i\mathrm{d}{x}^j.\hfill \end{array}$$(1)

Here, α is the lapse function, β_i are the shift vectors, h_ij = g_ij is the spatial part of the metric, Greek indices span space and time (taking values 0−3), Latin indices span only space (values 1−3), and the Einstein summation convention applies. The lapse function and shift vectors define fiducial observers (FIDOs), with four-velocities

$$\begin{array}{c}\hfill k^{\mu }=g^{\mu \nu }{k}_{\nu }=g^{\mu \nu }{(-\alpha ,0,0,0)}_{\nu }=(1,-{\beta }^1,-{\beta }^2,-{\beta }^3)/\alpha \end{array}$$(2)

normalized such that k^μk_μ = −1. Unless otherwise stated, we report all simulated quantities as measured by FIDOs.

Throughout this work, we adopt a Kerr spacetime expressed via horizon-penetrating spherical Kerr-Schild coordinates, (x⁰, x¹, x², x³) = (t, r, θ, ϕ). In these coordinates, the metric reads

$$\begin{array}{cccccc}\hfill \alpha & =\frac{1}{\sqrt{1+z}}\hfill & \hfill {\beta }^r& =\frac{z}{1+z}\hfill & {\beta }^{\theta }& ={\beta }^{\varphi }=0\\ \hfill h_{\mathit{rr}}& =1+z\hfill & \hfill h_{\theta \theta }& =\mathrm{\Sigma }\hfill & h_{\varphi \varphi }& =\frac{A{sin}^2\theta }{\mathrm{\Sigma }}\\ \hfill h_{r\varphi }& =-a{r}_g{sin}^2\theta (1+z)\hfill & \hfill h_{r\theta }& =h_{\theta \varphi }=0,\hfill \end{array}$$(3)

where

$$\begin{array}{cccc}\hfill \mathrm{\Sigma }& \equiv r^2+a^2{r}_g^2{cos}^2\theta \hfill & \hfill z& \equiv 2r_{g}r/\mathrm{\Sigma }\hfill \\ \hfill A& \equiv {(r^2+a^2{r}_g^2)}^2-a^2{r}_g^2\mathrm{\Delta }{sin}^2\theta \hfill & \hfill \mathrm{\Delta }& \equiv r^2-2r{r}_g+a^2{r}_g^2,\hfill \end{array}$$(4)

r_g = GM/c² is the gravitational radius, and M is the BH mass. The determinant of the spatial metric is h = Σ²(1 + z)sin²θ.

Maxwell’s equations can be stated covariantly as (Jackson 1975)

$$\begin{array}{cc}\hfill {\mathrm{\nabla }}_{\mu }{F}^{\nu \mu }& =\frac{4\pi }{c}{I}^{\nu }\mathrm{and}\hfill \\ \hfill {\mathrm{\nabla }}_{\mu }\ast F^{\nu \mu }& =0,\hfill \end{array}$$(5)

where ∇_μ is the covariant derivative, I^μ is the four-current density, and F^μν and *F^μν are, respectively, the electromagnetic field tensor and its Hodge dual. In the 3 + 1 formalism adopted in this paper, equations (5) are recast to closely resemble their flat-space form as (Komissarov 2004)

$$\begin{array}{cc}\hfill {\partial }_t\mathbf{B}& =-c\mathrm{\nabla }\times \mathit{E},\hfill \\ \hfill {\partial }_t\mathbf{D}& =c\mathrm{\nabla }\times \mathit{H}-4\pi \mathit{J},\hfill \\ \hfill \mathrm{\nabla }·\mathbf{D}& =4\pi \rho ,\mathrm{and}\hfill \\ \hfill \mathrm{\nabla }·\mathbf{B}& =0,\hfill \end{array}$$(6)

where we have assumed a stationary metric. Here, Dⁱ = F^iνk_ν = −αFⁱ⁰ and Bⁱ = − * F^iνk_ν = α * Fⁱ⁰ are, respectively, the electric and magnetic fields observed by FIDOs; ρ = −I^μk_μ/c is the FIDO-observed charge density, Jⁱ is a current density² given in terms of the timelike Killing vector ζ^μ = (1, 0, 0, 0) of the metric as Jⁱ = (I^μζⁱ − Iⁱζ^μ)k_μ = αIⁱ; and E and H are the auxiliary fields

$$\begin{array}{cc}\hfill \mathbf{E}& =\alpha \mathbf{D}+\beta \times \mathbf{B}\mathrm{and}\hfill \\ \hfill \mathbf{H}& =\alpha \mathbf{B}-\beta \times \mathbf{D}.\hfill \end{array}$$(7)

2.2. Initial conditions

Our simulations start in vacuum with D^r = D^θ = D^ϕ = B^ϕ = 0. We set the initial poloidal components of the magnetic field, B^r and B^θ, according to a prescribed magnetic flux function A_ϕ(r, θ, t = 0) such that

$$\begin{array}{cc}\hfill B^r(r,\theta ,0)& =\frac{1}{\sqrt{h}}{\partial }_{\theta }{A}_{\varphi }(r,\theta ,0)\mathrm{and}\hfill \\ \hfill B^{\theta }(r,\theta ,0)& =-\frac{1}{\sqrt{h}}{\partial }_r{A}_{\varphi }(r,\theta ,0).\hfill \end{array}$$(8)

We used this prescription to thread the equatorial boundary with a sequence of magnetic loops, as shown in Fig. 1.

We chose the functional form of A_ϕ(r, θ, 0) so that the loops all: have approximately the same diameter, λ = 10r_g; carry the same amount of magnetic flux; and alternate in magnetic polarity (the direction that the magnetic field circulates around each loop). These properties can all be achieved by selecting A_ϕ on the equator as

$$\begin{array}{c}\hfill A_{\varphi }(r,\pi /2,0)=\{\begin{array}{cc}{B}_0{r}_g^{2}sin(\pi r/\lambda ),\hfill & r<20r_g\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array},\end{array}$$(9)

and extending this into the rest of the simulation domain by requiring the resulting magnetic field to be potential. As a simplification, instead of using the exact potential-field that results from condition (9), we adopted the solution that would follow if r, θ, and ϕ were the standard flat-space spherical coordinates:

$$\begin{array}{c}\hfill A_{\varphi }(r,\theta ,0)=rsin\theta {\displaystyle {\int }_0^{\infty }}B(\kappa )J_1(\kappa rsin\theta )e^{-\kappa rcos\theta }\mathrm{d}\kappa ,\end{array}$$(10)

where

$$\begin{array}{c}\hfill B(\kappa )=\kappa {\displaystyle {\int }_0^{\infty }}{A}_{\varphi }(R,\pi /2,0)J_1(\kappa R)\mathrm{d}R\end{array}$$(11)

and J₁ is a Bessel function of the first kind. Because expression (10) ignores spacetime curvature, our initial magnetic field is not formally current free. This results in a mild transient that quickly evacuates our simulation in the first several r_g/c – well in advance of the times we analyze in this work. Our main objective in using expression (10) is not to start from a formal equilibrium, but rather to embed our simulation with an initially looped magnetic topology (Fig. 1).

2.3. Boundary conditions

At the θ = 0 and r = r_max boundaries, we employed axisymmetric and open boundary conditions, respectively. For the open boundary at r = r_max, this involves setting up a perfectly matched layer (PML; Birdsall & Langdon 1991; Berenger 1994; Cerutti et al. 2015) where the fields are smoothly damped to zero over a finite shell extending from r = r_pml = 27r_g to r = r_max = 30r_g. Particles that reach r_max are deleted. Our choice of r_max = 30r_g, limited by computational cost in PIC simulations, restricts our computational domain to the jet-launch region. To help circumvent this limitation, in Section 3.2.2 we extrapolate our results to larger radii via a phenomenological model.

We placed the inner radius, r = r_min, of our simulation underneath the BH event horizon, $r_H=r_g(1+\sqrt{1-a^2})\simeq 1.14r_g$, setting r_min = 0.9r_H. This conveniently renders the numerics insensitive to how the r = r_min surface is handled. For simplicity, we enforced a zero radial derivative here on D, E, B, and H, and we deleted all particles that touched r_min.

The main nontrivial boundary condition in our simulations is applied across a thin, one-cell thick disk lining the equator. This is the main control surface that we used to drive the simulation dynamics. Beginning at time t = 0, we threw the footpoints of magnetic field lines threading the disk into rotation at the local Keplerian angular velocity. We did this by applying a poloidal electric field, $E_r=\sqrt{h}{\mathrm{\Omega }}_{\mathrm{K}}{B}^{\theta }$ and $E_{\theta }=-\sqrt{h}{\mathrm{\Omega }}_{\mathrm{K}}{B}^r$, where

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{K}}=\sqrt{\frac{\mathit{GM}}{r_g^3}}\frac{1}{{(rsin\theta /r_g)}^{3/2}+a},\end{array}$$(12)

and r sin θ ≃ r is the distance to the symmetry axis.

Besides Keplerian rotation, we also prescribed a time-varying A_ϕ inside the equatorial disk, continuous with our initial condition (9), of the form

$$\begin{array}{c}\hfill A_{\varphi }(r,\pi /2,t)=\{\begin{array}{cc}{B}_0{r}_g^{2}sin\left(\frac{\pi }{\lambda }(r+v_{0}t)\right),\hfill & r<20r_g\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}.\end{array}$$(13)

This pulls the footpoints of the magnetic field lines threading the equator inward at speed v₀ = 0.01c, mimicking the dragging of magnetic loops by an underlying accretion disk. Our intent in choosing v₀ = 0.01c is not to model a real accretion speed. Instead, we merely aim to have the accretion time for a single loop, t_acc = λ/v₀, much longer than all the timescales governing the coronal dynamics, the longest of which is the time it takes to reconnect the flux in one loop, t_rec ∼ λ/β_recc, where β_rec ≃ 0.1 is a typical relativistic magnetic reconnection rate. As long as the hierarchy t_acc ≫ t_rec (i.e., v₀ ≪ β_recc ≃ 0.1c) is respected, the imposed inward field-line dragging just causes the simulation to pass through a sequence of quasisteady states, each one corresponding to a different set of magnetic footpoint locations.

2.4. Plasma supply

For simplicity, we adopted an ad hoc volumetric plasma injection scheme. At every timestep, for any cell in which the overall (electron+positron) plasma number density, n, fell below a prescribed floor of

$$\begin{array}{c}\hfill n_{\mathrm{fl}}=n_0{\left(\frac{r_g}{r}\right)}^2,\end{array}$$(14)

we injected electron-positron pairs to bring the density back up to n_fl. Plasma was injected with zero mean FIDO-measured velocity and with a FIDO-frame temperature kT = m_ec². The power law n_fl ∝ r⁻² maintains a roughly constant (modulo sinusoidal modulations) equatorial plasma magnetization, σ ≡ B_iBⁱ/4πnm_ec², since, from equations (8) and (13), the square of the vertical field threading the disk, ∼h_θθB^θB^θ, also decays with a 1/r² profile. We weighted the injected particles to achieve a simulation average of about 20 macroparticles per grid cell.

2.5. Synchrotron radiative cooling

We applied synchrotron radiative cooling to the plasma. We relegate a study involving inverse Compton (IC) cooling, which is also potentially important – and, indeed, perhaps even dominant – in this context to a future work. Our goal here is not to capture the radiative physics exactly, but merely to ensure that reconnection-accelerated particles promptly cool, radiating away their energy in the form of high-energy photons while remaining close to their acceleration sites. This is the expected situation in highly magnetized BH ADCs (Beloborodov 2017; Mehlhaff et al. 2021). We calibrated synchrotron cooling so that all particles with relativistic (FIDO-measured) Lorentz factors γ ≫ 1 cooled down to transrelativistic energies, γ ≳ 1, within a characteristic lightcrossing time, r_g/c, of the BH.

Adopting synchrotron cooling only as a proxy for the expected prompt radiative losses (possibly generated, in reality, via the IC mechanism) prevents us from predicting the observable spectrum of radiation. However, in the prompt-cooling regime, the bolometric high-energy luminosity simply tracks the instantaneous magnetic dissipation (e.g., Mehlhaff et al. 2024). Thus, we expect the spectrally integrated light curves presented in Section 3.3 to be somewhat insensitive to the choice (synchrotron or IC) of radiation mechanism. This should, however, be checked by a future work, since, even if different radiative processes preserve the bolometric luminosity, they can lead to other differences, such as in the angular distribution of emitted power away from the source.

2.6. Parameter values and grid selection

We chose the parameters B₀ and n₀ in conjunction with our simulation grid so that the plasma in our simulation (at least outside of reconnection sites) is highly magnetized and force free, as expected in BH ADCs. In addition, in order to avoid spurious numerical effects, we needed to ensure that important plasma microscales were resolved by our simulation grid. These objectives required us to meet three simultaneous conditions:

1. A healthy plasma supply, with enough particles to carry the current necessary for the force-free electromagnetic fields.

2. A high magnetization σ ≫ 1 (except potentially inside reconnecting current sheets).

3. A fine enough grid to resolve the plasma skin depth everywhere.

Conditions 1 and 3 are most stringent near the BH horizon, where the plasma is densest, the magnetic field strongest, and the rotation most rapid (i.e., shortest skin depth and highest demand for electrical charge). Therefore, meeting these two conditions closest to the BH suffices to meet them throughout our simulation. Condition 2, on the other hand, tends to become more strained at large r as σ decays. To ensure σ ≫ 1 at all radii, we set the fiducial magnetization near the BH horizon, σ₀ = B₀²/4πn₀m_ec², as high as we can afford while respecting conditions 1 and 3. We then checked a posteriori that σ ≫ 1 throughout our simulation box.

To ensure an adequate plasma supply, we chose n₀ as a multiple, η, of the fiducial Goldreich-Julian density, n_GJ = B₀Ω_BH/4πce, where Ω_BH = ac/2r_H is the BH angular velocity. We set η ≡ n₀/n_GJ = 3. This takes care of condition 1.

To handle conditions 2 and 3, we recast the magnetic field in terms of the nominal gyroradius ρ₀ ≡ m_ec²/eB₀, which permits us to rewrite the horizon-scale skin depth, d_e0 = (m_ec²/4πn₀e²)^1/2, and magnetization, σ₀ = B₀²/4πn₀m_ec², as

$$\begin{array}{c}\hfill d_{e0}^2=\frac{1}{\eta }\frac{m_e{c}^2}{e{B}_0}\frac{c}{{\mathrm{\Omega }}_{\mathrm{BH}}}=\frac{2{\rho }_0{r}_g}{\eta }\left(\frac{1+\sqrt{1-a^2}}{a}\right)\simeq 0.8{\rho }_0{r}_g\end{array}$$(15)

and

$$\begin{array}{c}\hfill {\sigma }_0=\frac{1}{\eta }\frac{e{B}_0}{m_e{c}^2}\frac{c}{{\mathrm{\Omega }}_{\mathrm{BH}}}=\frac{2}{\eta }\frac{r_g}{{\rho }_0}\left(\frac{1+\sqrt{1-a^2}}{a}\right)\simeq 0.8\frac{r_g}{{\rho }_0}·\end{array}$$(16)

The values a = 0.99 and η = 3 are assumed at the ends of both lines. These expressions show that the need for a well-resolved skin depth is in direct tension with that for a high magnetization, linking high σ directly to computational cost. To ease the numerical burden, we employed a logarithmically stretched grid, keeping Δ(lnr) constant, which concentrates resolution toward the BH. We then set ρ₀ = r_g/4000, yielding d_e0 ≃ 0.01r_g and σ₀ ≃ 3000 – sufficient to maintain σ ≫ 1 throughout the simulation domain – and we laid out N_r = 1024 cells in r and N_θ = 512 cells in θ (with Δθ constant). This grid results in a skin-depth resolution of d_e0/Δr ≃ 4 at r = r_H.

While the parameters above respect the expected high-σ, charge-satiated conditions in ADCs, their actual values are probably much lower than in real systems. For example, for an XRB containing a 10 M_⊙ BH with inner coronal magnetic field strength ∼10⁷ G (Mehlhaff et al. 2021), we would have r_g/ρ₀ ∼ 10¹⁰. Future, larger simulations could improve the separation between microphysical and global scales, but PIC simulations of this problem will probably always require a reduced scale separation.

3. Results

3.1. Overall dynamics

A snapshot showcasing the main dynamical aspects of our simulation is presented in Fig. 2. Rotational shear imposed by the equatorial disk and BH builds up an azimuthal component in the magnetic field of closed loops, pressurizing and inflating them toward infinity. As the loops inflate, sharp, current-sheet discontinuities develop in the magnetic field and begin reconnecting. Reconnecting current sheets separate differentially rotating open field lines, appearing as discontinuities in the map of field-line angular velocity, $\mathrm{\Omega }=-E_{\theta }/\sqrt{h}{B}^r$, in Fig. 2.

Fig. 2. Snapshot of the simulation analyzed in this work. Left: Local average particle Lorentz factor ⟨γ⟩. Right: Angular velocity, $\mathrm{\Omega }=-E_{\theta }/\sqrt{h}{B}^r$, of magnetic field lines. Rotational shear along magnetic field lines causes the initially closed magnetic topology to partially open. A fully open topology is prevented by reconnection, which allows some closed loops to persist. Reconnecting current sheets are visible as sites of particle acceleration (left) and discontinuities in field-line angular velocity (right). A movie of this simulation is available online.

The competition between rotational shear and reconnection sets the size of the open field-line regions (Uzdensky & Goodman 2008). On the one hand, rotation inflates magnetic loops, pumping energy into the field and promoting a progressively higher-energy open configuration (with the theoretical maximum-energy limit being a fully open field; Aly 1984, 1991; Sturrock 1991). On the other hand, reconnection allows the field to snap back, permitting closed field lines to persist and preventing a fully open (maximum-energy) configuration from being reached. An advantage of our fully kinetic model is that the reconnection rate is set self-consistently, and thus so too are the relative sizes of the open and closed field-line regions. Our finite simulation domain can potentially artificially accentuate open field-line regions, however, by cutting off reconnection at the outer boundary.

Despite rapid reconnection, prominent bundles of open magnetic flux persist, transmitting Poynting flux from their footpoints toward infinity, either in the form of a BZ jet funnel (for field lines attached to the BH; Blandford & Znajek 1977) or in the form of a force-free wind (for field lines attached to the disk; Blandford 1976). This transmission is imperfect, however, because some of the injected Poynting flux gets consumed by reconnection. The dissipated energy is used to accelerate particles and, thus, current sheets appear as thin, red-hot regions in Fig. 2. The energy budget of this system, especially the efficiency with which Poynting flux is transmitted to infinity versus consumed by reconnection, is the subject of Section 3.2. The observable synchrotron radiation produced by the accelerated particles is the focus of Section 3.3.

The open flux bundles in Fig. 2 have a nearly radial, monopolar shape at large r. This is similar to the final field structure predicted by earlier studies of axisymmetric force-free equilibria (Barnes & Sturrock 1972; Lynden-Bell & Boily 1994), wherein an initially closed magnetic field, line-tied to a differentially rotating disk, is quasistatically twisted open. Far away from the disk region where the flux is tied, the open field becomes nearly radial. Modulo reversals across current sheets, the strength of this field is approximately that of a monopole with flux equal to the total unsigned open flux threading the disk.

The sites of most intense particle acceleration in Fig. 2 coincide with reconnection along field lines coupling the disk to the BH. This indicates the BH-disk interaction as an important driver of the dynamical activity. Reconnection still occurs on disk-disk field lines, but it is less intense, owing to the lower field strength and rotational shear.

Over long timescales, visible in the animated movie to accompany Fig. 2, one can see the slow advection, at speed v₀ = 0.01c, of the magnetic footpoints toward the BH. Once the footpoints of the innermost loop get pushed onto the horizon, vigorous reconnection begins to annihilate the loop, ejecting its flux away in the form of a rapid-fire barrage of plasmoids launched toward infinity.

Figure 2 and its accompanying film also show, at the outer disk radii (close to r = 20r_g) the formation of new magnetic loops. This is a self-consistent byproduct of our dynamical driving (discussed in Section 2.3) of A_ϕ on the equatorial boundary of our simulation. Because new loops are formed in situ, we are able to simulate as many loop advection and ejection cycles as we wish without needing to fit all of the loops onto the equatorial disk initially. We are therefore able to present three such cycles, though we thread the disk with only two loops at t = 0.

3.2. Energy budget

Here, we analyze the energy budget in our system. First, in Section 3.2.1, we empirically track the fate of the Poynting flux injected by the equatorial disk and BH. We quantify how much of this Poynting flux is transmitted through our simulation versus converted to particle kinetic energy, and we characterize the efficiency with which a BZ jet is launched. Subsequently, in Section 3.2.2, we present a toy phenomenological model to explain the simulated dissipation in terms of the rate, β_rec, of relativistic magnetic reconnection. We summarize the key takeaways from Sections 3.2.1 and 3.2.2 in Section 3.2.3.

3.2.1. Measurements from simulation

Our main analysis technique to quantify the energy budget in our system is to evaluate Poynting’s theorem, which can be phrased in the 3 + 1 formalism as (Komissarov 2004)

$$\begin{array}{c}\hfill {\partial }_t\left[\frac{1}{8\pi }(\mathit{E}·\mathit{D}+\mathit{H}·\mathit{B})\right]+\mathbf{\nabla }·\left[\frac{c}{4\pi }(\mathit{E}\times \mathit{H})\right]=-\mathit{J}·\mathit{E}.\end{array}$$(17)

We integrated Eq. (17) over the surface shown in Fig. 3, decomposing into contributions from the terms defined in the same figure. These terms are: the Poynting flux injected from the disk and BH, ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}$; the Poynting flux escaping the simulation, ${\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{Poynt}}$; the Poynting flux consumed to energize the plasma, ${\dot{\mathcal{E}}}_{\mathrm{J}.\mathrm{E}}$; the rate of change of electromagnetic energy, ${\dot{\mathcal{E}}}_{\mathrm{EMfields}}$; and the Poynting-theorem residual (which should equal zero). We note that here, as throughout the remainder of the text, we refer to the quantities in Poynting’s theorem according to their familiar flat-space names, such as “Poynting flux” or “electromagnetic energy density”, even though equation (17) is written in terms of Poynting flux and electromagnetic energy density measured at infinity (as opposed to locally by FIDOs).

Fig. 3. Integration surface (pink) over which we evaluate Poynting’s theorem. It extends from r1 = rH up to r2 = 25rg < rpml and from θ = 0 to θ2 = 89° (just outside the equatorial disk). We decompose contributions to Poynting’s theorem into the terms shown. In the expressions, S ≡ cE × H/4π. We define ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}\equiv {\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{BH}}+{\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Disk}}$.

We present time series of the energy fluxes defined above in Fig. 4. We also present the time-averaged values of these fluxes in Table 1. In both the figure and table, we define L_BZ ≡ r_g⁴B₀²Ω_BH²/12c as the BZ power extracted from the upper hemisphere of a BH threaded by a monopolar magnetic field of horizon strength B₀ (Blandford & Znajek 1977; Tchekhovskoy et al. 2011; Crinquand et al. 2020). In our simulations, a monopole roughly describes the field lines at large r as well as the open field lines crossing the horizon.

Fig. 4. Top: Time series of the terms, defined in Fig. 3, contributing to Poynting’s theorem. Bottom: Total magnetic flux (Eq. (18)) piercing the BH horizon, including decomposed contributions (Eq. (19)) from open (BH-Inf) and disk-connected (BH-Disk) field lines.

In Fig. 4, the input and output Poynting fluxes evolve roughly in phase, with approximately two thirds of the injected Poynting flux being transmitted out of the simulation box. This agrees with the time-average result from Table 1 that $⟨{\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{Poynt}}⟩\simeq 2⟨{\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}⟩/3$. The remaining one third of the injected electromagnetic energy is used to energize the plasma through magnetic reconnection, with $⟨{\dot{\mathcal{E}}}_{\mathrm{J}.\mathrm{E}}⟩\simeq ⟨{\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}⟩/3$.

The contribution from ${\dot{\mathcal{E}}}_{\mathrm{EMfields}}$ to Poynting’s theorem goes to zero on time average, since we average, by design, over an integer number of advection cycles. Theoretically, we expect ${\dot{\mathcal{E}}}_{\mathrm{EMfields}}$ to approach zero at all times, not just on the time average, in the v₀ → 0 limit. Otherwise, the simulation would not attain a quasisteady state for each magnetic footpoint configuration. We have checked this explicitly by running simulations (not presented) with varying v₀, which show that a smaller value results in a smaller contribution from ${\dot{\mathcal{E}}}_{\mathrm{EMfields}}$ at any given time to the balance of Poynting’s theorem. The same exercise, on the other hand, does not change the time averages of the various $\dot{\mathcal{E}}$ terms presented in Table 1. We therefore conclude that the time-averaged values presented in Table 1 represent the true slow-advection limit.

The bottom panel of Fig. 4 shows the evolution of the instantaneous magnetic flux,

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

piercing the BH upper hemisphere. The Φ time series follows a clean sinusoid, reflecting the driving boundary condition (13). The energy budget is highly sensitive to Φ: more horizon flux generally corresponds to higher energy injection (${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}$), dissipation (${\dot{\mathcal{E}}}_{\mathrm{J}.\mathrm{E}}$), and transmission (${\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{Poynt}}$). This emphasizes the result already observed in the discussion of Fig. 2 (Section 3.1) that the dynamics are strongly driven by the magnetic coupling between the disk and the BH. Maxima in |Φ| represent moments where this coupling is quite efficient, powering the most vigorous reconnection and producing coincident maxima in ${\dot{\mathcal{E}}}_{\mathrm{J}.\mathrm{E}}$.

Notably, however, the BH flux, Φ, is not exactly in phase with the total injected power, ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}$. Moreover, the latter occasionally exceeds L_BZ. This indicates that the energetics of this system are not purely due to the BZ process; an appreciable amount of the injected energy also comes from the disk. To isolate the individual contributions from the BH and disk requires a more detailed energy budget assessment.

We therefore complement the Poynting theorem analysis presented above with a more fine-grained view. Here, we decompose the Poynting flux flowing into and out of our integration surface of Fig. 3 into contributions from different field-line bundles. To enable this analysis, we first classified all the points in our simulation domain based on their magnetic connectivity. Points lying on field lines connecting the disk to the BH are labeled “BH-Disk”; those on field lines coupling the disk to itself are labeled “Disk-Disk”; points on black-hole ingrown field lines as “BH-BH”; those inside open magnetic flux tubes anchored to the disk or BH as “Disk-Inf” or “BH-Inf”, respectively; and those on field lines tethered neither to the BH nor to the disk as “Detached”. This labeling scheme is illustrated in Fig. 5.

Fig. 5. Classification of spatial regions in our simulation domain based on magnetic connectivity. A movie of this figure is available online.

We note that ingrown BH-BH field lines are highly transient structures in our simulation. They only briefly occur at the moment when the final shred of magnetic flux in the innermost loop is ejected, by reconnection, to infinity: the corresponding morsel of reconnected flux still attached to the BH (the ingrown magnetic hair) quickly sinks through the horizon. In this way, the BH serves as a magnetic flux sink – a role it plays in addition to injecting electromagnetic energy via the BZ mechanism.

Our magnetic connectivity classification allows us to decompose the magnetic flux threading the BH into its contributions from BH-Inf (open) and BH-Disk (disk-connected) field lines. Formally, we define

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_b=2\pi {\displaystyle {\int }_0^{\pi /2}}{\mathrm{\Theta }}_b(r_H,\theta )B^r\sqrt{h}\mathrm{d}\theta ,\end{array}$$(19)

where Θ_b(r, θ) is equal to one if the point (r, θ) lies on flux bundle b (e.g., b = BH–Inf) and zero otherwise. The total horizon flux can be written Φ = ∑_bΦ_b = Φ_{BH − Inf} + Φ_{BH − Disk}. Here, the sum over field-line bundles recovers the only two sets of field lines, BH-Inf and BH-Disk, that, by definition, yield nonzero contributions to Φ. The lower panel of Fig. 4 shows time series of Φ_{BH − Inf} and Φ_{BH − Disk} alongside that of Φ. Loud periods in the overall dissipation, ${\dot{\mathcal{E}}}_{\mathrm{J}.\mathrm{E}}$, follow peaks in |Φ_{BH − Inf}|. This shows that the most intense magnetic reconnection in our system preys on the jet funnel, eating away at the BZ jet to power particle acceleration and (as we shall see in Section 3.3) bright high-energy emission.

In Fig. 6, we use our magnetic connectivity classification to analyze the energy flux flowing into and out of our Poynting integration surface (Fig. 3) on the Disk-Inf and BH-Inf field-line bundles. From here onward, we refer to these two bundles, respectively, as the “disk wind” and the “jet”. Although we focus on just the jet and disk wind here, for completeness, we present in Appendix A the energy budget including contributions from all the magnetic connectivity regions.

Figure 6 shows that the transmission of Poynting flux through the system, in both the jet and wind, is extremely well correlated with injection from the corresponding energy source: either the BH or the disk. The transmission efficiency for each source is nearly time-independent, with the jet and wind field lines both relaying about two thirds of their initial Poynting flux to the outer edge of the domain. This is, of course, also true of the time-averaged transmission efficiencies reported in Table 2. We further note from Fig. 6 that, while the disk wind is roughly stationary, the jet turns on and off with a duty cycle of roughly one half. This permits the jet to temporarily become more powerful than the wind even though the wind carries more Poynting flux on average.

Fig. 6. Top: Time series of Poynting flux entering (escaping) our integration surface – depicted in Fig. 3 – on open field lines threading the BH, ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{BH}-\mathrm{Inf}}$ (${\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{BH}-\mathrm{Inf}}$), and disk, ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Disk}-\mathrm{Inf}}$ (${\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{Disk}-\mathrm{Inf}}$). Bottom: Time series of open (BH-Inf), closed (BH-Disk), and total horizon-piercing magnetic flux as in Fig. 4.

The lower panel of Fig. 6 shows, like Fig. 4, the time series of Φ, Φ_{BH − Inf} and Φ_{BH − Disk}. The injected jet power, ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{BH}-\mathrm{Inf}}$, is very well correlated with the open magnetic flux, Φ_{BH − Inf}, through the horizon. This is precisely what one expects from the BZ mechanism, here clearly exhibited thanks to our field-line classification scheme. We explicitly demonstrate that the BZ correlation, ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{BH}-\mathrm{Inf}}\propto {\mathrm{\Phi }}_{\mathrm{BH}-\mathrm{Inf}}^2$, is obeyed in Fig. 7. The proportionality constant can be explained by noting that ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{BH}-\mathrm{Inf}}$ should approach f(a)L_BZ, where f(a) = 1 + 1.38(Ω_BHr_g/c)² − 9.2(Ω_BHr_g/c)⁴ is the high-spin correction factor identified by Tchekhovskoy et al. (2010), as Φ_{BH − Inf} tends toward the full fiducial flux 2πB₀r_g². This gives ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{BH}-\mathrm{Inf}}/f(a)L_{\mathrm{BZ}}={({\mathrm{\Phi }}_{\mathrm{BH}-\mathrm{Inf}}/2\pi B_0{r}_g^2)}^2$, as in Fig. 7.

Fig. 7. Injected Poynting power versus magnetic flux carried by open field lines threading the BH. Each data point represents one snapshot of the simulation. Color indicates time phase-folded on the loop advection and ejection period, 1000rg/c. The expected correlation, explained in the text, is drawn in black with f(a) evaluated as f(0.99) = 0.93.

3.2.2. A toy dissipation model

To take a closer look at the physics underlying our dissipation measurements, we present radial dissipation profiles from our simulation in Fig. 8. Each profile is defined with respect to an integration surface similar to that depicted in Fig. 3 as

$$\begin{array}{c}\hfill \dot{ϵ}(r)\equiv {\dot{ϵ}}_{\mathrm{inj}}(r)-{\displaystyle {\int }_0^{2\pi }}{\displaystyle {\int }_0^{\pi /2}}{S}^r(r,\theta \prime )\sqrt{h}\mathrm{d}\theta \prime \mathrm{d}\varphi \prime ,\end{array}$$(20)

where

$$\begin{array}{cc}\hfill {\dot{ϵ}}_{\mathrm{inj}}(r)\equiv & {\displaystyle {\int }_0^{2\pi }}{\displaystyle {\int }_0^{\pi /2}}{S}^r(r_H,\theta \prime )\sqrt{h}\mathrm{d}\theta \prime \mathrm{d}\varphi \prime \hfill \\ \hfill & -{\displaystyle {\int }_{r_H}^r}{\displaystyle {\int }_0^{2\pi }}{S}^{\theta }(r\prime ,{\theta }_2)\sqrt{h}\mathrm{d}r\prime \mathrm{d}\varphi \prime \hfill \end{array}$$(21)

is the Poynting flux injected by the BH plus the disk up through radius r. In the above expressions, as in Fig. 3, S = cE × H/4π and θ₂ = 89°.

Computing dissipation profiles using Eq. (20) instead of directly via the J ⋅ E term of Eq. (17) misses the contribution from the time-derivative term in Poynting’s theorem. However, this term disappears when we average over an integer number of loop advection cycles. Thus, after averaging, the profiles defined by Eq. (20) reconstruct the exact J ⋅ E dissipation.

Understanding the profiles $\dot{ϵ}(r)$ presented in Fig. 8 is the main goal of this section. To do so, we formulate below a toy dissipation model heavily inspired by that presented by Cerutti et al. (2020) for the case of reconnection in pulsar winds. Within this model, the logarithmic dependence, $\dot{ϵ}(r)\propto ln(r/r_0)$, of the profile follows naturally from magnetic reconnection. The proportionality constant is, moreover, the collisionless relativistic reconnection rate, β_rec, well known to be of order 0.1.

Fig. 8. Dissipation profiles computed according to Eq. (20). Instantaneous profiles are faint black lines; the time-averaged (over 200 ≤ ct/rg ≤ 3200) profile is a thick black line. We fit Eq. (28), with βrec a free parameter, to the time-averaged profile between r = r0 = 4rg and r = 16rg. The fit is drawn as a dashed red line between r = 4rg and r = 16rg with a constant vertical offset for clarity.

To begin with, we concentrate on dissipation within the innermost current sheet along the jet funnel. Dissipation farther out – in current sheets along disk-disk field lines – also occurs, but it does not consume Poynting flux from the jet. Moreover, outer current sheets power very little high-energy radiation (Section 3.3). This is because the local magnetic field strength and accelerated particle energies (Fig. 2 and its accompanying film) are lower, diminishing the energies of the emerging photons.

As can be seen from Figs. 2 and 5 and their accompanying films, reconnection kicks in at a finite distance, r₀, away from the BH, typically around 4r_g or so. Since this is fairly far from the event horizon, we approximate space as flat throughout the dissipation zone. We therefore use, for the remainder of this section, only the fields B and E and work exclusively in flat orthonormal spherical coordinates (all lower indices with, e.g., E² = E_r² + E_θ² + E_ϕ²).

In addition, as previously observed in Section 3.1, the poloidal magnetic field far away from the BH is nearly monopolar. We therefore approximate the radial magnetic field at r > r₀ = 4r_g as that of a magnetic monopole:

$$\begin{array}{c}\hfill B_r=\pm B_{r_0}\frac{r_0^2}{r^2},\end{array}$$(22)

where B_r0 is the radial field strength at r = r₀ and the overall sign changes across reconnection current sheets. Within the jet funnel, this field rotates at the constant angular velocity Ω = Ω_BH/2 (Fig. 2), which can be associated with a light cylinder at r sin θ = R_LC = 2c/Ω_BH ≃ 4r_H/a ≃ 4r_g. The azimuthal field produced by this rotation is

$$\begin{array}{c}\hfill B_{\varphi }=\mp \frac{rsin\theta }{R_{\mathrm{LC}}}{B}_r=\mp B_{r_0}\frac{r_0^2}{R_{\mathrm{LC}}}\frac{sin\theta }{r}·\end{array}$$(23)

For r sin θ > R_LC, the field is mostly azimuthal (|B_ϕ|> |B_r|).

We assume that the jet-funnel current sheet is a thin cone-like slab with opening angle θ_CS, where θ_CS, based on Figs. 2 and 5 and their accompanying films, is close to 45°. The slab extends outward starting from r = r₀ and has a full angular thickness Δθ = δ/r (constant physical thickness δ). Dissipation causes the profile $\dot{ϵ}(r)$ to grow from its initial value, $\dot{ϵ}(r_0)$, at the beginning of the current sheet by the amount

$$\begin{array}{cc}\hfill \dot{ϵ}(r)-\dot{ϵ}(r_0)& ={\displaystyle {\int }_0^{2\pi }}{\displaystyle {\int }_0^{\pi /2}}{\displaystyle {\int }_{r_0}^r}\mathit{E}·\mathit{J}{r}^{\prime 2}sin\theta \mathrm{d}r\prime \mathrm{d}\theta \mathrm{d}\varphi \hfill \\ \hfill & \simeq 2\pi \delta {\displaystyle {\int }_{r_0}^r}\mathit{E}·\mathit{J}r\prime sin{\theta }_{\mathrm{CS}}\mathrm{d}r\prime .\hfill \end{array}$$(24)

In the second step, the integral over θ only activates across a thin region of thickness Δθ = δ/r′ centered on θ = θ_CS.

To evaluate the remaining radial integral in Eq. (24), we need to determine the electric field and current density inside the current sheet. Since the condition r sin θ > R_LC holds throughout most of the current sheet, where θ = θ_CS ≃ 45° and r ≥ r₀ = 4r_g ≃ R_LC, the reconnecting magnetic field is mostly azimuthal. This field reverses by an amount ΔB_ϕ ∼ 2B_ϕ across the reconnection layer, requiring a radial current density of

$$\begin{array}{c}\hfill J_r\simeq \frac{c\mathrm{\Delta }{B}_{\varphi }}{4\pi \delta }\sim \frac{c{B}_{\varphi }}{2\pi \delta }·\end{array}$$(25)

Meanwhile, the reconnection electric field in the heart of the layer is E_r = β_recB_ϕ, where β_rec ∼ 0.1 is the collisionless relativistic reconnection rate.

Plugging these estimates for E_r and J_r into equation (24) gives the dissipation profile

$$\begin{array}{cc}\hfill \dot{ϵ}(r)-\dot{ϵ}(r_0)& \simeq c{\beta }_{\mathrm{rec}}{B}_{r_0}^2\frac{r_0^4}{R_{\mathrm{LC}}^2}{sin}^3{\theta }_{\mathrm{CS}}{\displaystyle {\int }_{r_0}^r}\frac{\mathrm{d}r\prime }{r\prime }\hfill \\ \hfill & =c{\beta }_{\mathrm{rec}}{B}_{r_0}^2\frac{r_0^4}{R_{\mathrm{LC}}^2}{sin}^3{\theta }_{\mathrm{CS}}ln\left(\frac{r}{r_0}\right).\hfill \end{array}$$(26)

The thickness, δ, of the current sheet cancels out. To simplify expression (26), we normalize it by the Poynting power, ${\dot{ϵ}}_{\mathrm{inj}}(r_0)$, injected at r = r₀. Since this radius is before substantial dissipation occurs, we can set $\dot{ϵ}(r_0)$ to zero in Eq. (20), which yields

$$\begin{array}{cc}\hfill {\dot{ϵ}}_{\mathrm{inj}}(r_0)& =\frac{c}{4\pi }{\displaystyle {\int }_0^{2\pi }}{\displaystyle {\int }_0^{\pi /2}}{(\mathit{E}\times \mathit{B})}_r{r}_0^{2}sin\theta \mathrm{d}\theta \mathrm{d}\varphi \hfill \\ \hfill & =\frac{c}{2}{\displaystyle {\int }_0^{\pi /2}}{B}_{\varphi }^2{r}_0^{2}sin\theta \mathrm{d}\theta =\frac{c}{3}{B}_{r_0}^2\frac{r_0^4}{R_{\mathrm{LC}}^2}·\hfill \end{array}$$(27)

Using this result to normalize Eq. (26), we can write

$$\begin{array}{c}\hfill \frac{\dot{ϵ}(r)-\dot{ϵ}(r_0)}{{\dot{ϵ}}_{\mathrm{inj}}(r_0)}\simeq 3{\beta }_{\mathrm{rec}}{sin}^3{\theta }_{\mathrm{CS}}ln\left(\frac{r}{r_0}\right)\simeq {\beta }_{\mathrm{rec}}ln\left(\frac{r}{r_0}\right).\end{array}$$(28)

In the final step, we noted, using θ_CS ≃ 45°, that 3sin³θ_CS ≃ 1.

In Fig. 8, we fit Eq. (28) to the time-averaged dissipation profiles measured from our simulation. Unlike in our toy model, the measured profiles exhibit some dissipation at radii smaller than r₀ = 4r_g. This near-horizon dissipation is probably spurious, related to our loop advection prescription that starts to lose realism there. Thus, when fitting to the measured profiles, we consider only radii r > 4r_g. In addition, the dissipation profiles all steepen beyond r = 16r_g or so. Such steepening is expected due to dissipation in the outer current sheet, which is barely contained in our simulation. To avoid this steepening influencing the fit, we consider only r < 16r_g.

The measured profiles of Fig. 8 are normalized by the injected energy ${\dot{ϵ}}_{\mathrm{inj}}(r_2)$ up to outer radius r₂ = 25r_g. Poynting’s theorem dictates that ${\dot{ϵ}}_{\mathrm{inj}}(r_2)=\dot{ϵ}(r_2)-\dot{ϵ}(r_0)+{\dot{ϵ}}_{\mathrm{inj}}(r_0)$. Plugging this into Eq. (28) yields

$$\begin{array}{c}\hfill \frac{\dot{ϵ}(r)-\dot{ϵ}(r_0)}{{\dot{ϵ}}_{\mathrm{inj}}(r_2)}=\frac{{\beta }_{\mathrm{rec}}ln(r/r_0)}{1+{\beta }_{\mathrm{rec}}ln(r_2/r_0)}\simeq {\beta }_{\mathrm{rec}}ln\left(\frac{r}{r_0}\right),\end{array}$$(29)

where, in the last step, we used β_rec ∼ 0.1 to estimate 1 + β_recln(r₂/r₀)≃1.2 ≃ 1. Thus, normalization by ${\dot{ϵ}}_{\mathrm{inj}}(r_2)$ as in Fig. 8 does not substantially change the fitted β_rec values.

Extrapolating Eq. (28) beyond our numerical domain implies that our simulation underestimates the total dissipation. At extremely large distances, our toy model even implies that the injected Poynting flux becomes fully dissipated at r = r_diss = r₀exp(1/β_rec) ∼ 10⁵r_g (between 10³r_g and 10⁹r_g assuming 0.05 < β_rec < 0.2). Taking into account reconnection in the outer current sheet would pull in the radius of complete dissipation even closer. However, since the outer current sheet cannot siphon Poynting flux off of the jet, we think that this estimate of r_diss is appropriate for determining the scale where the jet becomes fully dissipated.

However, whether complete dissipation truly occurs at such enormous radii is quite uncertain. At these scales (10⁵r_g = 5 pc for a 10⁹ M_⊙ BH), the jet may have already started interacting with the ambient medium, modifying our predicted profile (28), if not disrupting the jet entirely. Even without an ambient medium, the jet could still become intrinsically unstable to disruption, for example by kink modes.

We conclude, therefore, that our simulation strictly underestimates total dissipation, much of which occurs at larger, unsimulated radii. Indeed, the toy model of this section implies that reconnection along the jet wall provides an efficient mechanism for converting an initially Poynting-flux-dominated flow into particle kinetic energy and radiation. However, whether Poynting flux dissipation continues according to profile (28) all the way to total depletion depends on additional factors (such as the ambient medium and current-driven instabilities) that are beyond the scope of our work.

3.2.3. Energy budget summary

We wish to emphasize the following key takeaways from our measurements (Section 3.2.1) and model (Section 3.2.2) of the energy budget in our simulation:

1. Our simulation dissipates about one third of the Poynting flux injected by the equatorial disk and the BH; it lets the remaining two thirds pass through the box.

2. Loop inflation leads to open magnetic flux bundles attached to both the BH and the disk. We refer to these bundles, respectively, as the jet and disk wind.

3. The jet and wind individually transmit roughly two thirds – equivalent to about 0.2L_BZ and 0.3L_BZ, respectively – of their initial Poynting flux out of the simulation.

4. The wind is quasistationary, but the jet flickers.

5. Jet flickering results from alternating periods of open magnetic flux accumulating on the BH and reconnection eating that flux away.

6. We present a model where reconnection on the jet wall yields the radial Poynting flux dissipation profile, β_recln(r/r₀), with β_rec ∼ 0.1 the relativistic reconnection rate.

7. Our model suggests that dissipation would continue at radii beyond our simulated domain. Thus, points 1 and 3 above must be interpreted as lower bounds on energy dissipation and upper bounds on Poynting flux transmission.

We caution that all of these results apply to the regime of our simulation, involving loops much larger than the BH: λ = 10r_g ≫ r_H. In a companion publication (Crinquand et al., in prep.), we illustrate the effects of a smaller loop size, λ ≲ r_H, a regime in which many of the points above no longer hold.

3.3. Synthetic observables

We calculated synchrotron emission from the particles in our simulation and ray-traced (Crinquand et al. 2021, 2022) the resulting photons to observers at infinity. Here, we focus on the high-energy radiation, which is produced by particles accelerated via reconnection. To avoid contamination from unaccelerated particles (with order-unity Lorentz factors), we only counted photons with energy-at-infinity exceeding 10ℏeB₀/m_ec. In addition, we ignored photons produced at altitudes lower than r cos θ = r_H above the disk. This excludes spurious near-horizon emission resulting from our ad hoc accretion prescription at the equatorial boundary. Figure 9 presents the light curve of bolometric luminosity, L_bol, of all escaping photons (i.e., neither absorbed by the BH nor by the equatorial disk) as well as the simultaneous horizon magnetic flux calculated per Eqs. (18) and (19).

Fig. 9. Top: Bolometric luminosity calculated for photons escaping to infinity with energy exceeding 10ℏeB0/mec. Here, ct/rg corresponds to time of photon emission; not time of photon reception. Light-travel-time effects are thus absent in this light curve. Bottom: Instantaneous open (BH-Inf), closed (BH-Disk), and total magnetic flux threading the BH horizon, as in Fig. 4.

The loop advection and ejection cycles, with period λ/v₀ = 1000r_g/c, have a dramatic effect on L_bol. Bolometric output peaks near maximum |Φ|. This corresponds to the moment when reconnection is actively shredding away the open field lines attached to the BH, leading to vigorous particle acceleration. As a result, even though peak brightness periods align with maxima in the total flux, |Φ|, they follow slightly behind peaks in the open flux, |Φ_{BH − Inf}|. Magnetic reconnection stochastically ejects the flux in each loop in the form of plasmoids, creating jagged short-timescale features on the falling side of each luminosity peak. Once the innermost loop has been completely ejected (Φ ≃ 0), the luminosity goes through a minimum. From there, it slowly grows again as the next loop is brought into place. This gradual brightening results from reconnection intensifying again as the differential rotation and field strength in the next loop become more pronounced at progressively smaller radii. Loud periods (|Φ| maxima) and quiet periods (Φ ≃ 0) alternate on the loop-advection period, λ/v₀, and exhibit a remarkable brightness contrast of 10³.

Besides the bolometric light curve, we also present, in Fig. 10, light curves measured by distant observers at different inclination angles, i, measured from the BH spin axis. The main features present in the L_bol time series are also present in the light curves measured by individual observables: on secular timescales (of order the loop advection period, 1000r_g/c), a brightness contrast of 10³ or more between maximum and minimum activity levels; on short timescales and particularly on the falling segment of the luminosity envelope, abrupt reconnection-driven subflaring.

Fig. 10. Light curves measured by distant observers at different inclinations i. The average light curve over all observers between i = 0° and i = 90° is also shown for reference (i-avg.). The time-axis is reported in retarded time, correcting for the light-travel-time to the observer (placed an arbitrary distance robserver away from r = 0).

Despite overall qualitative agreement with the main features of Fig. 9, the observer-specific light curves of Fig. 10 are quantitatively more nuanced depending on inclination i. Observers (e.g., i = 40.5°) looking along the brightest inner current sheet on the jet wall witness faster and more intense variability. For example, in the inset of Fig. 10, the i = 40.5° light curve fluctuates by up to a factor of 10 on timescales as short as ∼r_g/c.

This is a light-travel-time effect, similar to what occurs in blazars: because the emitting plasma travels nearly at the speed of light toward the observer, the radiation is beamed and the arrival times of the photons are compressed, rendering taller and sharper the associated subflares in the light curve. The relativistic motion here probably has both fluid and kinetic origins. On the fluid level, the motion of fast plasmoids, especially small ones, carries the plasma radially outward at relativistic speeds (Sironi et al. 2016; Petropoulou et al. 2016). There is also likely a contribution from the kinetic beaming effect, wherein the highest energy particles (and, hence, their radiation) are selectively beamed along the reconnection outflow (Cerutti et al. 2012, 2013, 2014; Mehlhaff et al. 2020).

We complement the light curves in Fig. 10 with corresponding power spectral densities (PSDs) in Fig. 11. Each PSD is defined as dP/df  ≡  (|X(f)| + |X(−f)|)/2, where $X(f)\equiv {\displaystyle {\int }_{-\infty }^{\infty }}L(t)exp(-2\pi ift)\mathrm{d}t$ is the Fourier Transform of the corresponding observer’s light curve, L(t). Our secular loop-advection timescale 1000r_g/c can clearly be seen in Fig. 11: the spectral power for each PSD is concentrated at frequencies f ∼ 10⁻³c/r_g. In contrast to this long-timescale narrow-frequency driving, reconnection produces much faster and broader-band variability, resulting in extended tails in the PSD for each observer. The privileged i = 40.5° observer that looks directly along the inner reconnection current sheet witnesses the most rapid variability timescales. Hence, the i = 40.5° PSD contains an excess at high frequencies that extends all the way up to the Nyquist limit of c/2r_g.

Fig. 11. Power spectral densities (PSDs), dP/df, for observers at different inclinations i. These PSDs are calculated from the light curves in Fig. 10, including the one averaged over all inclination angles. Each PSD is compensated by f, which shows the spectral power per logarithmic frequency interval, and normalized by P ≡ ∫0∞dP/df df.

We cannot fully determine the longest variability timescales that reconnection can produce. This is because, at low frequencies in Fig. 11, the narrow-frequency feature associated with the loop-advection timescale, λ/v₀ = 1000r_g/c, overlaps the broadband reconnection-powered variability. A low-frequency cutoff on the power-law PSD produced by reconnection is a natural outcome of the self-similar plasmoid hierarchy (Shibata & Tanuma 2001; Uzdensky et al. 2010), where the slowest variability timescales are governed by the growth and ejection of the largest plasmoids the current sheet can produce (Petropoulou et al. 2016). Definitively measuring such a cutoff would require pushing the loop-driving feature in our PSDs to lower frequencies, demanding either a slower advection speed, v₀, or a larger loop size, λ, and, hence, more expensive simulations. In the absence of such simulations, we infer a conservative bound on the slowest variability frequencies yielded by reconnection. This bound can be deduced from the excess in the i = 40.5° PSD of Fig. 11. Because the excess persists down to frequencies f ∼ 10⁻²c/r_g and is linked exclusively to anisotropy (beaming and variability compression) associated with reconnection, we infer that reconnection can drive variability on timescales at least as long as ∼100r_g/c.

In summary, the high-energy synchrotron observables from our simulations feature:

1. A large, secularly driven brightness contrast of at least 10³ (possibly more depending on observer inclination) between bright and quiet phases of the loop advection and ejection cycle.

2. Peak intensity periods coincident with peak magnetic flux threading the BH.

3. Rapid subflaring due to reconnection between the BH and equatorial disk. Subflaring begins in the bright period, close to peak magnetic flux on the horizon, and continues during the transition to the quiet phase, as the horizon flux wanes.

4. Enhanced variability at inclinations matching the opening angle of the innermost current sheet. Subflares at these inclinations can rise by up to an order of magnitude in as little as ∼r_g/c.

5. Reconnection-driven broadband variability extending across at least two decades in timescales, from ∼r_g/c up through ∼100r_g/c. Costlier simulations could determine if reconnection can yield even slower variability, expanding this range.

4. Discussion

4.1. Hard-to-soft X-ray binary state transitions

Outbursting BH XRBs pass through a sequence of empirically defined states: first the hard state, where the X-ray spectrum is dominated by nonthermal (power-law) emission peaking in the hard X-rays; next an intermediate state where the overall luminosity is high and the spectrum begins to soften; subsequently the soft state, where a soft, ∼1 keV quasithermal spectrum dominates the X-ray emission; finally, an eventual return to the hard state (Remillard & McClintock 2006; Done et al. 2007). These X-ray states also have radio counterparts. During the hard state, a compact, steady jet is detected in the radio; in the soft state, the radio emission is consistent with no jet; the hard-to-soft transition is heralded by one or more discrete radio ejections before the jet eventually extinguishes as the binary settles into the soft state (Fender et al. 2004, 2009; Miller-Jones et al. 2012; Russell et al. 2019; Wood et al. 2023).

We would like to suggest a picture where radio ejections observed during hard-to-soft XRB transitions stem from the launching of magnetic loops as in our simulation: one radio ejection per loop. In this picture, magnetic reconnection links X-ray flaring to radio ejecta by playing dual roles. First, it powers high-energy emission (Section 3.3). Second, it also eats away at the magnetic field lines forming the jet funnel (Section 3.2.1), eventually completely untying their footpoints from the BH, and, as a result, ejecting the magnetic flux originally contained in each coronal loop toward infinity.

In such a scenario, radio emission is ultimately powered by the Poynting flux injected into the jet by the BH. Observations suggest that this Poynting flux must silently propagate for some distance before dissipating. This is because discrete radio ejections generally light up after the main X-ray transition: already at a considerable distance from the BH. It is only by backtracing the motion of the ejecta across the sky that a launch time nearly coincident with the X-ray transition can be inferred (Miller-Jones et al. 2012; Russell et al. 2019; Wood et al. 2021, 2023).

In the context of our model, the movement of the radio ejecta seems to preclude their production by continuous reconnection along the jet wall (Section 3.2.2), which would yield a standing structure. Moreover, such reconnection must not entirely consume the available Poynting energy in the jet. Some of it must remain to allow for more sudden and explosive dissipation by different mechanisms farther downstream.

Among the possible downstream dissipation mechanisms are shocks that form as the ejected material interacts with the ambient medium or catches up to slower previously launched ejecta (Jamil et al. 2010; Malzac 2014). Alternatively, given that the polarity of each ejected loop alternates, a sequence of ejected loops would, at large scales, produce a magnetically striped jet. The Poynting flux in such a jet could be dissipated by direct reconnection between adjacent stripes (Drenkhahn 2002; Drenkhahn & Spruit 2002; Lyubarsky 2010; Giannios & Uzdensky 2019; Zhang & Giannios 2021). We return to elaborate this point in Section 4.3.

Besides the change to their X-ray spectrum, BH XRBs also frequently exhibit a few characteristic X-ray timing properties during hard-to-soft state transitions. Often their X-ray variability diminishes and quasi-periodic oscillations (QPOs) at low frequencies (≲10 Hz) in their X-ray power spectral densities change form (i.e., type B QPOs may appear; Fender et al. 2009; Belloni 2010; Ingram & Motta 2019). Of these features, the reduction in X-ray variability would be expected from our model. As seen in the light curves of Figs. 9 and 10, the high-energy radiation powered by magnetic reconnection dims as the flux associated with the innermost loop is ejected away. In such a case, one expects the steadier radiation produced by the underlying accretion disk (not modeled in our setup) to stabilize the overall X-ray output.

All of these remarks seem to produce a rather consistent picture of X-rays and radio ejecta in hard-to-soft XRB transitions being coproduced by the reconnection-mediated ejection of coronal loops. However, there are a few important caveats to keep in mind, particularly with respect to X-ray timing. First, X-ray variability, including low-frequency QPOs, associated with hard-to-soft transitions tends to be characterized at frequencies below 10 Hz or so (Remillard & McClintock 2006; Belloni 2010). This is much slower than the reconnection-driven variability we measure in Section 3.3, which occurs from frequencies as high as c/r_g ∼ 10⁴ Hz down to 10⁻²c/r_g ∼ 10² Hz (assuming a 10 M_⊙ BH mass). As was mentioned in Section 3.3, probing whether reconnection can yield variability on even slower timescales requires more expensive simulations.

Additionally, given the absence of discrete features in our synthetic PSDs (Fig. 11), QPOs are probably not driven by reconnection. Instead, they likely stem from physics outside our model, in particular from effects linked to the accretion disk such as Lense-Thirring precession (Stella & Vietri 1998; Fragile et al. 2007; Ingram & Motta 2019; Liska et al. 2021). Such processes, though not originating in the corona, could perturb it, potentially leading to QPO signatures in the hard X-rays.

4.2. The peculiar changing-look AGN 1ES 1927+654

Some AGNs undergo rapid changing-look events, where, broadly speaking, their spectral types suddenly change from Type I – defined by clear Doppler-broadened (∼10³ km s⁻¹) emission lines – to Type II – characterized by the absence of broad lines – or vice versa (e.g., Tohline & Osterbrock 1976; LaMassa et al. 2015; MacLeod et al. 2016; Yang et al. 2018, and references therein). In December 2017, 1ES 1927+654 became the first AGN to be observed changing look in real time (Trakhtenbrot et al. 2019; Ricci et al. 2020). Among the several peculiar aspects of its behavior, perhaps the most striking is that, t ≃ 160 days after its initial optical/UV brightening by a factor of 10^1 − 2, the X-ray luminosity dropped by nearly four orders of magnitude, bottoming out at t ≃ 200 days before rebounding, by t ≃ 300 days, to just beyond its pre-flare level. Throughout this dramatic drop-off and rebound, the intraday X-rays were themselves variable by nearly two orders of magnitude (Trakhtenbrot et al. 2019; Ricci et al. 2020).

Scepi et al. (2021) and Laha et al. (2022) speculated that the peculiar X-ray activity of 1ES 1927+654 observed in 2017−2018 could be the result of the destruction and reformation of an X-ray emitting magnetized ADC. The optical/UV brightening, in this scenario, signals an increase in the accretion rate in the underlying disk, which inwardly advects magnetic flux of opposing polarity to that of the field that threads the disk initially. The fresh magnetic flux first reconnects with the originally present flux, destroying the magnetized corona, and then accumulates in the inner region, reconstructing a corona of opposite magnetic polarity to the first one.

The scenario sketched by Scepi et al. (2021) bears a striking resemblance to our simulation. One might imagine, for example, that the coronal configuration at the beginning of the event is similar to that at the beginning of our simulation: the BH is initially threaded by a large magnetic flux, with reconnection on BH-disk field lines powering the variable X-ray emission. The optical/UV brightening signals an increase in the accretion speed, v₀. The inner loop then gets pushed onto the BH and ejected, causing a dip in the X-ray emission. Finally, as the next loop is brought into place, the X-rays rebound. Essentially, the observations probe two peak activity states in our simulation: for example, from 200 to 1500 in the abscissa of Fig. 10.

We wish to stress that the absolute time elapsed, t_acc = λ/v₀ = 10³r_g/c, during this sequence of events in our simulation is not very important; it just reflects our artificially fast accretion speed, v₀ = 0.01c. Indeed, for a central BH mass of 2 × 10⁷ M_⊙ (Trakhtenbrot et al. 2019), the observed duration of the X-ray drop-off and rebound, Δt_obs ∼ 100 d ∼ 10⁵r_g/c, is much longer than t_acc = 10³r_g/c. More relevant is the agreement in amplitude of the long-timescale X-ray modulation between our simulation and the changing-look event: in both cases the variation is by a factor of 10^3 − 4. In our model, this modulation is dictated by the coronal response to the slow magnetic footpoint motion in the disk. Insofar as the response is quasistatic, the precise duration of the loop-ejection cycle is irrelevant, and so a comparison focusing on the amplitude of the X-ray modulation is justified.

Unlike the long timescale, t_acc, imposed artificially by our loop-driving prescription, the faster variability generated by reconnection in our model (e.g., in the inset of Fig. 10) is completely physical. We would like to suggest that such reconnection-powered subflaring is responsible for the fastest variability reported by Ricci et al. (2020): the intra-day X-ray fluctuations. The problem here is that the observing cadence, roughly 1 d⁻¹ ∼ 10⁻³c/r_g (Ricci et al. 2021), was too slow to probe the rapid timescales expected from reconnection, which can be as short as a few r_g/c (Section 3.3). In our proposed scenario, then, the observations represent essentially uncorrelated reconnection-powered fluctuations³. Focusing again, therefore, on a comparison of variability amplitudes instead of direct timescales, the intra-day X-ray variability by a factor of 10^1 − 2 agrees well with that in our simulated light curves by up to a factor of 10 (e.g., the i = 40.5° light curve in the inset of Fig. 10). All-in-all, from the standpoint of variability amplitude on both long and short timescales, our simulations seem to corroborate the scenario of Scepi et al. (2021).

Caveats here include that no discrete radio ejections were observed to accompany the optical/UV/X-ray activity (Laha et al. 2022; though we note that such ejections, along with QPOs, did recently accompany an increase in X-ray activity in the same object; Meyer et al. 2025; Laha et al. 2025). Furthermore, our simulations resemble the case of a classic, optically thick, geometrically thin accretion disk (Shakura & Sunyaev 1973) sandwiched by a magnetized corona. The accretion timescales for such a disk around a supermassive BH are known to be far too long to explain a changing-look timescale of several months (Dexter & Begelman 2019). Nevertheless, recent MHD simulations (Jacquemin-Ide et al. 2024) show evidence that thicker disks can efficiently produce and advect inward their own magnetic field loops, producing a circumdisk magnetic configuration that is actually not so different from our model.

4.3. Striped jets

In our model, reconnection allows a BH to untie magnetic field lines from its accretion disk and twist them up into a BZ-like jet. In each loop accretion cycle, the resultant jet funnel comprises a magnetic helix of opposing polarity to the preceding one. A series of loop advection and ejection cycles thus produces a corresponding train of alternating-polarity magnetic slabs in the jet. These slabs are initially very elongated along the jet axis. However, as they expand along the collimation profile of the jet, they spread out transversely while maintaining an approximately fixed spacing in the longitudinal direction. The far downstream region of the jet thus becomes magnetically “striped”, with axially thin but laterally wide stripes of reversing magnetic polarity (Giannios & Uzdensky 2019): much like a striped pulsar wind (Coroniti 1990).

As the stripe aspect ratio thins, the stripes become increasingly prone to reconnect with one another. If strong inter-stripe reconnection can be triggered, it drives particle acceleration (with commensurate radiation) and bulk jet acceleration in tandem. In this way, magnetic striping could explain a number of key aspects of AGN and gamma-ray-burst jets, including the need for a high jet Lorentz factor, in situ particle acceleration far downstream of the central engine, and a broad range of variability timescales (Drenkhahn 2002; Drenkhahn & Spruit 2002; Giannios 2006; Lyubarsky 2010; Giannios & Uzdensky 2019; Zhang & Giannios 2021).

However, due to the tremendous range of scales covered by BH jets (from the event horizon to beyond the host galaxy for the most extreme AGNs), global numerical striped jet models have so far proven elusive. Computational work has instead focused on key subcomponents of the problem, especially whether accreting magnetic loops can indeed translate into powerful stripes in the jet (Parfrey et al. 2015; Yuan et al. 2019; Mahlmann et al. 2020; Nathanail et al. 2020, 2022; Chashkina et al. 2021). Our work complements these preceding efforts by being the first to analyze striped-jet launching using a first-principles GRPIC framework (though see also El Mellah et al. 2022, 2023, similar studies without jet striping). Our simulations thus incorporate completely self-consistent energetics. These enable us to confirm a key finding of previous work: the advection of large (λ ≫ r_g) loops leads to efficient jet activation (Parfrey et al. 2015; Mahlmann et al. 2020; Chashkina et al. 2021).

In addition, the self-consistent particle acceleration and radiation in our simulations enable us to report a clear observable counterpart to the striped-jet launching mechanism. Magnetic reconnection, as it launches each loop (i.e., eventual stripe) into the jet funnel, consumes a good fraction of the overall energy budget: even slightly more than the Poynting flux that is transmitted into the jet (since $⟨{\dot{\mathcal{E}}}_{\mathrm{E}.\mathrm{J}}⟩>⟨{\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{BH}\to \infty }⟩$ in Tables 1 and 2). This powers efficient nonthermal particle acceleration, leading to rapidly variable high-energy emission accompanying each stripe ejection.

While reconnection creates a potentially observable high-energy counterpart to jet striping, it also dissipates the remaining available Poynting flux in the stripes rather efficiently (Section 3.2.2), rendering them progressively more kinetic-energy-dominated. This opens up the possibility that the stripes – if produced as in our simulation – may become depleted of their magnetic energy before reaching distances where they could otherwise reconnect efficiently. However, we cannot constrain this possibility with high confidence, since our toy model (Section 3.2.2) that we use to extrapolate reconnection-powered dissipation to larger scales ignores other factors that become important there.

Other caveats on the implications of our results in the context of striped jets include our 2D setup, which cannot probe important effects such as the potential disruption of the nascent jet due to non-axisymmetric current-driven (i.e., kink) instabilities (Yuan et al. 2019; Mahlmann et al. 2020). Furthermore, in this preliminary study, we have neither examined retrograde disks, which could lead to even more efficient jet launching (Parfrey et al. 2015), nor the impact of varying loop sizes, especially λ ≲ r_g. We shall report the effects of alternate loop sizes in a forthcoming companion publication (Crinquand et al., in prep.).

4.4. Sgr A* flares

The supermassive BH at the center of our Galaxy, Sagittarius A* (Sgr A*), exhibits dramatic variability in the infrared (IR) and X-ray bands (Witzel et al. 2018). The GRAVITY Collaboration (2018) recently associated IR flaring with the orbital motion of hot spots close to the BH. El Mellah et al. (2022, 2023) proposed a model for the GRAVITY observations based on GRPIC simulations very similar to ours (see also Aimar et al. 2023). Within their model, reconnection occurs on field lines coupling the BH to its accretion disk, just like the innermost current sheet in our simulations. As a result, particles are accelerated and ejected away from the BH in the form of hot plasmoids (flux ropes in 3D). El Mellah et al. (2023) explored the hypothesis that such flux ropes are the structures observed by GRAVITY as hot spots.

As discussed by El Mellah et al. (2023), one difficulty with this scenario is that the synchrotron power from a single flux rope decays too quickly to explain the observed flare duration of 30 − 60 min, corresponding to 90 − 180r_g/c for a BH mass of 4.1 × 10⁶ M_⊙. They therefore speculated that the GRAVITY hot spots are conglomerates of multiple flux ropes. Individual flux ropes could then explain the subflaring observed on timescales of ∼10 min (Genzel et al. 2003; Dodds-Eden et al. 2011), or about 30r_g/c.

Our results reveal a mechanism through which a multiple-flux-rope hot spot could produce the right range of variability timescales, as speculated by El Mellah et al. (2023). In Section 3.3, we show that the ensemble of particle acceleration sites and plasmoids (2D flux ropes) in the reconnection layer generates a range of variability timescales from r_g/c up to 100r_g/c (though we cannot rule out still longer timescales). This roughly coincides with the observed range: from 30r_g/c to 90 − 180r_g/c. The absence of observed variability as fast as r_g/c could stem from the viewing-angle effect identified in Section 3.3: only observers looking along the current sheet witness pronounced variability at the fastest r_g/c timescales.

While our results might ease the tension in variability timescales between observations and the model of El Mellah et al. (2022, 2023), there are a number of remaining issues. For example, as El Mellah et al. (2023) also point out, the open spiral trajectories traced by plasmoids away from the BH make it difficult to reproduce the observed closed loops in the Stokes (Q, U)-plane (GRAVITY Collaboration 2018; Wielgus et al. 2022). Moreover, the need for an ensemble of flux ropes to form an orbiting hot spot requires strong non-axisymmetry in the reconnection current sheet. This might arise from a non-axisymmetric embedding of magnetic field lines in the disk, but it is a possibility that neither we nor El Mellah et al. (2022, 2023) explore. Thus, whether the scenario discussed here underlies Sgr A* flares remains unclear.

5. Conclusions

In this paper, we use the GRZELTRON code (Parfrey et al. 2019) to present an axisymmetric GRPIC model of a BH fed by its ADC. This is the first PIC model to include the coupling between a magnetized plasma corona and a BH as well as an accretion prescription. Our model confirms and extends previous similar numerical studies based on force-free electrodynamics (Parfrey et al. 2015; Yuan et al. 2019; Mahlmann et al. 2020). In particular, we confirm the activation of the BZ mechanism for the case when large – but nevertheless zero-net-flux – coronal magnetic loops (with diameters λ ≫ r_g) are fed to a central rapidly spinning BH. We extend preceding works by employing a first-principles plasma description, which allows us to compute fully self-consistent energy budget diagnostics (Section 3.2) and radiative signatures (Section 3.3).

In Section 3.2, we trace the flow of Poynting flux through our simulation, evaluating Poynting’s theorem both across our whole simulation domain and along specific field-line bundles. We also present a phenomenological model for the Poynting flux dissipation through magnetic reconnection. Globally, we find that about one third of the Poynting flux injected into our simulation (by the BH plus the disk) is consumed through coronal reconnection. The other two thirds are transmitted through our numerical domain via a BZ jet (for open field lines connected to the BH) and a force-free wind (for open field lines attached to the disk). This 2:1 breakdown between dissipation and transmission also turns out to characterize the fate of the Poynting flux injected separately along jet field lines and disk-wind field lines.

However, we wish to discourage the interpretation of this 2:1 transmission-dissipation partition as a quantitative prediction. Instead, this ratio must be interpreted in light of our analytic model of Section 3.2.2. That model reproduces the radial profile of the fraction of dissipated Poynting flux, β_recln(r/r₀), from our simulation, where β_rec ≃ 0.1 is the relativistic magnetic reconnection rate. Extrapolating this formula implies that dissipation should continue at radii beyond our numerical domain, albeit more slowly (scaling logarithmically in r). At the same time, it is not clear that the profile, β_recln(r/r₀), remains accurate all the way out to the critical radius, r_diss = r₀exp(1/β_rec), where it predicts complete dissipation. A number of factors could intervene before then: the ambient medium could disrupt the force-free jet propagation; the jet could become intrinsically unstable (e.g., to kink modes); or the jet magnetization could decline as the result of bulk acceleration, altering the reconnection rate in situ. In view of these uncertain factors, we wish to cast the 2:1 transmission-dissipation ratio witnessed in our simulation simply as an upper bound. Dissipation continues beyond the simulation, scaling initially as β_recln(r/r₀) but eventually changing shape due to physics relevant at larger scales.

We wish to underscore that, while the final stages of jet dissipation remain uncertain, our work reveals a mechanism through which a nascent jet can begin efficiently converting Poynting flux to particle kinetic energy and radiation straight from launch. The triggering of some additional dissipation mechanism at larger scales, though certainly possible, is, in principle, not needed. Left to their own devices, jets launched via the inflation of a zero-net-flux magnetic loop, as in our simulation, will completely dissipate their energy by the time they reach r_diss.

One practical consequence of these remarks is that we cannot guarantee – though we also cannot exclude – that the magnetically striped BZ jet produced by our simulation retains a large reservoir of magnetic energy to be consumed later on (Giannios & Uzdensky 2019, see discussion in Section 4.3). However, we do witness a very robust electromagnetic counterpart to striped-jet launch. As magnetic reconnection injects successive stripes into the jet, it also powers particle acceleration and bright high-energy radiation at the jet base. High-energy flaring may thus herald jet striping. We characterize such flaring in the synchrotron channel in Section 3.3.

In particular, in Section 3.3, we reconstruct and analyze high-energy synchrotron light curves from our simulation. These are computed using the general relativistic ray-tracing module developed for GRZELTRON by Crinquand et al. (2021, 2022). Independently of inclination angle, the light curves show a dramatic brightness contrast of at least 10³ on slow timescales linked to the loop advection-ejection period. The closer a given loop is accreted in radius, the stronger the magnetic field and the rotational shear across its footpoints become. For our large loop sizes (λ = 10r_g), this means that the high-energy radiation from reconnection on the innermost loop, which couples the BH to the accretion disk, dominates that from all the other loops. This radiation peaks just after the maximum in open magnetic flux piercing the horizon, when vigorous reconnection preys on and whittles down the jet funnel. Reconnection along the jet wall is the main source of high-energy coronal output (a picture also recently advocated by Sridhar et al. 2025).

The secular brightness modulations on the loop advection and ejection period provide an envelope on top of which is superimposed much faster stochastic variability associated with reconnection. These rapid stochastic fluctuations are more sensitive than the secular envelope to observer inclination. For observers looking along the bright jet-wall reconnection layer, radiation is Doppler-compressed, akin to what occurs in blazars. This intensifies the stochastic variability component, leading to subflares occurring on timescales as fast as r_g/c and with brightness contrasts of up to a factor of 10.

The superimposed secular and stochastic variations in our simulations of, respectively, a factor of 10³ and 10 are on approximately the right levels to explain the secular and stochastic X-ray variability observed in the peculiar changing-look event of the AGN, 1ES 1927+654 (Trakhtenbrot et al. 2019). Thus, our model fits well with the picture advanced by Scepi et al. (2021) and Laha et al. (2022) in which this changing-look event corresponds to a change in polarity of the ADC associated with the AGN. This scenario is discussed in greater detail in Section 4.2.

Additionally, we measure (Section 3.3) a range of reconnection-powered variability timescales: from r_g/c at the shortest up to (at least) 100r_g/c at the longest. As discussed in Section 4.4, this range of timescales is roughly commensurate with that observed in the context of Sgr A* IR flares (Genzel et al. 2003; Dodds-Eden et al. 2011; GRAVITY Collaboration 2018). Our work thus corroborates the hypothesis of El Mellah et al. (2023), where a bundle of reconnecting flux powers the observed range of variability timescales. However, the difficulty of a reconnection layer – as studied by El Mellah et al. (2022, 2023) and by us here – to reproduce the orbiting hot spots measured by GRAVITY Collaboration (2018) and their concomitant polarization signatures (Wielgus et al. 2022) remains unresolved.

Our work reveals an important connection between coronal activity and jet launch: BZ jet power (Section 3.2.1) and coronal emission (Section 3.3) are coupled through a single reconnection layer. Thus, reconnection-mediated ejection of a coronal loop by a BH might explain radio-jet ejections accompanying X-ray outbursts during XRB hard-to-soft state transitions (Fender et al. 2004, see Section 4.1). Such a situation also naturally gives rise to the reduction in broadband X-ray variability that typically coincides with these transitions (Ingram & Motta 2019). As the stochastically fluctuating (Section 3.3) hard X-rays associated with coronal reconnection diminish, the presumably steadier soft X-ray radiation coming from the disk will dominate the light curve. However, discrete features (QPOs) seen in X-ray power spectral densities during XRB state transitions are likely produced by some other mechanism, perhaps linked to the accretion disk, that we do not model.

Our work focuses on the regime expected to launch powerful jets: that is, where the accreting magnetic loops are much larger than the central BH. In a companion publication (Crinquand et al., in prep.), we present the complementary and highly qualitatively different limit where the loops are smaller with respect to the BH. Besides the large loop size, a strong assumption of our work is that of axisymmetry. Relaxing this assumption in the future will permit an assessment of the impact of non-axisymmetric modes on jet launch (Yuan et al. 2019; Mahlmann et al. 2020) and on the reconnection layer itself. Another extension of our work would be to add polarization to the ray-tracing analysis. This would help constrain our model in light of recent X-ray polarimetry measurements by the IXPE instrument (e.g., Krawczynski et al. 2022; Chattopadhyay et al. 2024). Finally, and perhaps most excitingly, our model represents a fundamental link between BH ADCs on the one hand and the jet properties and dissipation that manifest at large scales on the other. Elucidating this link even further via first-principles modeling – and with the continued support of the complementary force-free and MHD approaches – is paramount in the quest to reveal the mysteries of BH jets more broadly.

Data availability

Movies associated with Figs. 2 and 5 are available at https://www.aanda.org

Acknowledgments

The authors would like to thank the anonymous reviewer for providing a thorough report that improved the quality of the manuscript. The authors further gratefully acknowledge stimulating discussions with Jonathan Ferreira, Jeremy Goodman, Claire Guépin, Gilles Henri, Geoffroy Lesur, Krzysztof Nalewajko, Pierre-Olivier Petrucci, Nicolas Scepi, Adrien Soudais, Dmitri Uzdensky, and Greg Werner. In addition, the authors thank Adrien Soudais for generously providing the colormap for the left-hand panel of Fig. 2. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 863412). Computing resources were provided by the TGCC under the allocation A0150407669 made by GENCI.

References

Appendix A: Complete Poynting theorem analysis

Here, we complement our analysis of Section 3.2.1, providing the complete decomposition of the energy budget in our simulation based on the identification of field line bundles illustrated in Fig. 5. We present, in Fig. A.1, the time-averaged values of every term in Poynting’s theorem at hierarchical granularity levels.

As shown by Fig. A.1, a good fraction of the Poynting flux injected by the BH gets directly channeled along BH-disk coupling field lines to the disk. Though the BH (BH box in the figure) injects 64 percent of the overall input Poynting flux into the simulation, 24 percent of this is immediately deposited, via BH-disk coupling field lines (BH-disk boxes in the figure) to the disk. Thus, while naively one might conclude that the BH is responsible for the majority (64 percent) of the simulation energy budget, in reality, of the energy that gets processed through the corona – either reaching the edge of the box on open field lines or being dissipated as J ⋅ E heating – only 64 − 24 = 40 percent comes from the BH, while 36 + 24 = 60 percent comes from the disk. In our companion publication (Crinquand et al. in prep.), we show that for smaller loop sizes λ ≲ r_H, this effect becomes more extreme: a larger fraction of the magnetic flux on the BH remains closed, and so most of the energy injected by the BH is immediately absorbed by the disk. In this way, for smaller loop sizes, closed field lines connecting the BH to the disk choke jet field lines, and with them the BH’s ability to transmit energy to infinity.

Fig. A.1. In each box, averages – taken over 200 ≤ ct/rg ≤ 3200 – of contributions to Poynting’s theorem expressed as a fraction of the total input Poynting flux, $⟨{\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}⟩=0.82L_{\mathrm{BZ}}$ (cf. Table 1). Energy sources are on the left of the central hatched bar; energy sinks are on its right. The decomposition becomes more fine-grained toward the middle. Percentages indicate branching ratios between granularity levels. For example, from the BH-Inf box and its connection to the Poynting box on the energy sink side, one reads $⟨{\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{BH}-\mathrm{Inf}}⟩=0.21⟨{\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}⟩=0.33⟨{\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{Poynt}}$). The sum of values in every column equals unity to within percent-level error, as required for energy conservation. We exclude contributions from BH-BH field lines and the rate-of-change of the electromagnetic field energy (${\dot{\mathcal{E}}}_{\mathrm{EMfields}}$) because these vanish at the percent level. Except in the cases of colored arrows crossing the hatched bar, it is not possible to determine branching ratios from an arbitrary source to an arbitrary sink. To represent this, black arrow-connections to the hatched bar indicate a sink-ambiguous deposit to (for the energy sources), or a source-ambiguous withdrawal from (for the energy sinks), the available energy budget.

All Tables

All Figures

	Fig. 1. Initial conditions (IC) and boundary conditions (BC) of our simulation setup. Sketch is not to scale.
In the text

Fig. 2. Snapshot of the simulation analyzed in this work. Left: Local average particle Lorentz factor ⟨γ⟩. Right: Angular velocity, $\mathrm{\Omega }=-E_{\theta }/\sqrt{h}{B}^r$, of magnetic field lines. Rotational shear along magnetic field lines causes the initially closed magnetic topology to partially open. A fully open topology is prevented by reconnection, which allows some closed loops to persist. Reconnecting current sheets are visible as sites of particle acceleration (left) and discontinuities in field-line angular velocity (right). A movie of this simulation is available online.

In the text

Fig. 3. Integration surface (pink) over which we evaluate Poynting’s theorem. It extends from r1 = rH up to r2 = 25rg < rpml and from θ = 0 to θ2 = 89° (just outside the equatorial disk). We decompose contributions to Poynting’s theorem into the terms shown. In the expressions, S ≡ cE × H/4π. We define ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}\equiv {\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{BH}}+{\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Disk}}$.

In the text

	Fig. 4. Top: Time series of the terms, defined in Fig. 3, contributing to Poynting’s theorem. Bottom: Total magnetic flux (Eq. (18)) piercing the BH horizon, including decomposed contributions (Eq. (19)) from open (BH-Inf) and disk-connected (BH-Disk) field lines.
In the text

	Fig. 5. Classification of spatial regions in our simulation domain based on magnetic connectivity. A movie of this figure is available online.
In the text

Fig. 6. Top: Time series of Poynting flux entering (escaping) our integration surface – depicted in Fig. 3 – on open field lines threading the BH, ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{BH}-\mathrm{Inf}}$ (${\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{BH}-\mathrm{Inf}}$), and disk, ${\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Disk}-\mathrm{Inf}}$ (${\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{Disk}-\mathrm{Inf}}$). Bottom: Time series of open (BH-Inf), closed (BH-Disk), and total horizon-piercing magnetic flux as in Fig. 4.

In the text

	Fig. 7. Injected Poynting power versus magnetic flux carried by open field lines threading the BH. Each data point represents one snapshot of the simulation. Color indicates time phase-folded on the loop advection and ejection period, 1000rg/c. The expected correlation, explained in the text, is drawn in black with f(a) evaluated as f(0.99) = 0.93.
In the text

Fig. 8. Dissipation profiles computed according to Eq. (20). Instantaneous profiles are faint black lines; the time-averaged (over 200 ≤ ct/rg ≤ 3200) profile is a thick black line. We fit Eq. (28), with βrec a free parameter, to the time-averaged profile between r = r0 = 4rg and r = 16rg. The fit is drawn as a dashed red line between r = 4rg and r = 16rg with a constant vertical offset for clarity.

In the text

	Fig. 9. Top: Bolometric luminosity calculated for photons escaping to infinity with energy exceeding 10ℏeB0/mec. Here, ct/rg corresponds to time of photon emission; not time of photon reception. Light-travel-time effects are thus absent in this light curve. Bottom: Instantaneous open (BH-Inf), closed (BH-Disk), and total magnetic flux threading the BH horizon, as in Fig. 4.
In the text

	Fig. 10. Light curves measured by distant observers at different inclinations i. The average light curve over all observers between i = 0° and i = 90° is also shown for reference (i-avg.). The time-axis is reported in retarded time, correcting for the light-travel-time to the observer (placed an arbitrary distance robserver away from r = 0).
In the text

	Fig. 11. Power spectral densities (PSDs), dP/df, for observers at different inclinations i. These PSDs are calculated from the light curves in Fig. 10, including the one averaged over all inclination angles. Each PSD is compensated by f, which shows the spectral power per logarithmic frequency interval, and normalized by P ≡ ∫0∞dP/df df.
In the text

Fig. A.1. In each box, averages – taken over 200 ≤ ct/rg ≤ 3200 – of contributions to Poynting’s theorem expressed as a fraction of the total input Poynting flux, $⟨{\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}⟩=0.82L_{\mathrm{BZ}}$ (cf. Table 1). Energy sources are on the left of the central hatched bar; energy sinks are on its right. The decomposition becomes more fine-grained toward the middle. Percentages indicate branching ratios between granularity levels. For example, from the BH-Inf box and its connection to the Poynting box on the energy sink side, one reads $⟨{\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{BH}-\mathrm{Inf}}⟩=0.21⟨{\dot{\mathcal{E}}}_{\mathrm{inj},\mathrm{Poynt}}⟩=0.33⟨{\dot{\mathcal{E}}}_{\mathrm{out},\mathrm{Poynt}}$). The sum of values in every column equals unity to within percent-level error, as required for energy conservation. We exclude contributions from BH-BH field lines and the rate-of-change of the electromagnetic field energy (${\dot{\mathcal{E}}}_{\mathrm{EMfields}}$) because these vanish at the percent level. Except in the cases of colored arrows crossing the hatched bar, it is not possible to determine branching ratios from an arbitrary source to an arbitrary sink. To represent this, black arrow-connections to the hatched bar indicate a sink-ambiguous deposit to (for the energy sources), or a source-ambiguous withdrawal from (for the energy sinks), the available energy budget.

In the text