Synthetic gamma-ray lightcurves of Kerr black-hole magnetospheric activity from particle-in-cell simulations

Context: The origin of ultra-rapid flares of very high-energy radiation from active galactic nuclei remains elusive. Magnetospheric processes, occurring in the close vicinity of the central black hole, could account for these flares. Aims: We aim to bridge the gap between simulations and observations by synthesizing gamma-ray lightcurves in order to characterize the activity of a black-hole magnetosphere, using kinetic simulations. Methods: We perform global axisymmetric two-dimensional general-relativistic particle-in-cell simulations of a Kerr black-hole magnetosphere. We include a self-consistent treatment of radiative processes and plasma supply, as well as a realistic magnetic configuration, with a large-scale equatorial current sheet. We couple our particle-in-cell code with a ray-tracing algorithm, in order to produce synthetic lightcurves. Results: These simulations show a highly dynamic magnetosphere, as well as very efficient dissipation of the magnetic energy. An external supply of magnetic flux is found to maintain the magnetosphere in a dynamic state, otherwise the magnetosphere settles in a quasi-steady Wald-like configuration. The dissipated energy is mostly converted to gamma-ray photons. The lightcurves at low viewing angle (face-on) mainly trace the spark gap activity and exhibit high variability. On the other hand, no significant variability is found at high viewing angle (edge-on), where the main contribution comes from the reconnecting current sheet. Conclusions: We observe that black-hole magnetospheres with a current sheet are characterized by a very high radiative efficiency. The typical amplitude of the flares in our simulations is lower than what is detected in active galactic nuclei. Such flares could result from the variation of parameters external to the black hole


Introduction
Ground-based Cherenkov telescopes have shown that active galactic nuclei (AGN) can be highly variable sources of very high-energy (VHE) emission (> 100 GeV) (Aharonian et al. 2007;Albert et al. 2007;Aleksić et al. 2014). Variability timescales can be shorter than the horizon light-crossing time t g = r g /c, where r g is the gravitational radius of the central supermassive black hole. This constrains emission models, as the size of the emitting region must be on the order of r g by virtue of causality.
This variability was primarily observed in blazars, AGN with jets lying close to our line of sight. However, TeV γ-ray flares were also detected from the nuclei of the radio galaxies M87 (Aharonian et al. 2006;Aliu et al. 2012) and Centaurus A (Aharonian et al. 2009) for instance, these galaxies having jets misaligned with our line of sight by more than 15 • . This suggests that VHE flares are a widespread feature in AGN. The variability timescale ∆t = 2 days of M87* VHE flares is comparable to the horizon light-crossing time t g 0.4 days. M87 is of paramount importance, and has attracted considerable attention because it is close enough that its nucleus can be resolved by radio interferometry (Event Horizon Telescope Collaboration et al. 2019). This has allowed observers to establish a connection between VHE flares and a brightening of the radio core (Acciari et al. 2009), which occurred simultaneously on several occasions. These observations also ruled out other potential compact VHE emission sites, such as knots in the relativistic jet. Thus, we might be able to link the formation of such a jet with processes at play in the close vicinity of the central black hole.
The extreme variability of these flares challenges conventional models of AGN (Rieger & Aharonian 2012). Because this emission seems to originate from the vicinity of the black hole, we are motivated to study nonthermal magnetospheric processes as a possible source (Katsoulakos & Rieger 2018;Levinson & Rieger 2011;Hirotani & Pu 2016;Levinson 2000). This model is applicable to low-luminosity AGN, such as M87* or Sgr A*. If the luminosity of the accretion flow is low enough, the plasma density can drop below the Goldreich-Julian value, which is required to screen the electric fields generated by the dragging of the magnetic field lines by the black hole (Wald 1974). Consequently, spark gaps arise, accelerating particles to very high energies; these energetic particles then scatter off soft photons from the accretion flow to produce VHE emission. Subsequent pair production is triggered by the annihilation of TeV photons with soft photons. This fresh plasma screens the electric field, quenching particle acceleration and nonthermal radiation. As the plasma inevitably flows away from the gap, the electric field is restored and VHE emission resumes. Hence, depending on the gap size, this model may account for the variability of the VHE emission observed from AGN. An electromagnetic cascade takes hold, feeding on the black hole to produce VHE emission and pair plasma. An important takeaway from this model is the possibility to activate the Blandford-Znajek (BZ) process by providing the plasma necessary to establish a quasi-force-free magnetosphere (Blandford & Znajek 1977). Since the presence of a gap has an impact on the global structure of the magnetosphere, and because the electromagnetic cascade is a highly nonlinear phenomenon, numerical simulations are well suited to gain insight into this problem. Parfrey et al. (2019) demonstrated the feasibility of performing global general relativistic particle-in-cell (GRPIC) simulations of a black hole magnetosphere. They modeled a black hole immersed in an initially vertical magnetic field, but used a simplified treatment for plasma supply that could only mimic how an electromagnetic cascade develops. This was improved upon in Levinson & Cerutti 2018, Chen & Yuan 2020, and Crinquand et al. 2020, where a self-consistent treatment of inverse Compton (IC) scattering and pair production was implemented. In Crinquand et al. (2020, hereafter C20), we simulated a monopole magnetosphere to capture the intrinsic activity of spark gaps, and showed that the BZ process could be successfully activated, as the magnetosphere was filled with pair plasma produced in the ergosphere. The size of the gap was consistent with sub-horizon variability.
However, in the case of isolated magnetospheres, more realistic configurations with a large-scale poloidal magnetic field should display an equatorial current sheet (Komissarov 2004;Komissarov & McKinney 2007). This current sheet would originate from the need to close the electric current system, since negative currents flow from both poles (if the spin axis is aligned with the magnetic field). Such a situation can come up if the accretion flow is truncated at large radius, causing the accretion to pause for a while, as can happen in magnetically arrested disk simulations (Narayan et al. 2003). Still, large-scale and intermittent ergospheric current sheets are expected to develop naturally in accreting black hole magnetospheres as well (e.g., Ripperda et al. 2020), highlighting the need to understand their importance. Magnetic reconnection, an intermittent phenomenon that is known to accelerate particles very efficiently, is ubiquitous in such current sheets (Kagan et al. 2015). It could also be responsible for variable VHE emission (Cerutti et al. 2012;Christie et al. 2019;Mehlhaff et al. 2020). It is unclear how such a current sheet can affect the pair discharge mechanism, and what the relative contributions of the polar cap and the current sheet emissions are. In this paper we extend the work carried out in C20, no longer neglecting the equatorial reconnection activity. Now that the polar cap activity has been characterized, we can aim toward more realistic magnetic configurations. In addition, contrary to previous kinetic studies, here we make a point of extracting physical observables from our simulations by implementing a more efficient treatment of ray tracing. This allows us to synthesize gamma-ray light curves, assimilating HE and VHE emission to IC processes. This paper is divided into two main sections. In Sect. 2 we describe our numerical scheme and setup, and present our new simulations of magnetospheric activity. In Sect. 3 we focus on producing observables from PIC simulations. We present synthetic light curves for the simulations carried out in Sect. 2 and in C20. The post-process treatment is described in Appendix A.

Numerical techniques
In this section only we set c = 1.
Metric In this work we perform 2D global GRPIC simulations, using a general relativistic version of the PIC code Zeltron (Cerutti et al. 2013, first introduced in Parfrey et al. (2019). The background spacetime is described by the Kerr metric, with a spin parameter a ∈ [0, 1[. The code uses spherical Kerr-Schild coordinates (t, r, θ, ϕ), which are not singular at the event horizon (see Komissarov 2004 for an expression of the coefficients of the metric). We use the 3+1 formulation of general relativity (MacDonald & Thorne 1982) in order to evolve the particles and fields with respect to a universal coordinate time t. This formulation naturally introduces fiducial observers (FIDOs), whose wordlines are orthogonal to the spatial hypersurfaces of constant t as privileged observers. In an axisymmetric and stationary spacetime such as the one described by the Kerr metric, they are also zero angular momentum observers (ZAMOs). In this formulation, the metric can be rewritten such that the line element ds 2 reads (1) In this equation α is the lapse function (the redshift of a FIDO with respect to the coordinate time t), β is the shift vector (the 3-velocity of a FIDO with respect to the coordinate grid), whereas h i j denotes the spatial 3-metric associated with the spatial hypersurfaces of constant t. The gravitational radius of the black hole is denoted r g .

Electromagnetic fields
We solve the electromagnetic field equations derived by Komissarov (2004) where B and D are respectively the magnetic and electric fields locally measured by FIDOs (they are physical observables), and H and E are auxiliary fields defined by The auxiliary current density J is related to the electric current density measured by FIDOs j via J = α j − ρβ, ρ being the electric charge density measured by FIDOs. These fields are defined on a spatial Yee grid (Yee 1966). Equations (2) and (3) resemble classical Maxwell-Ampère and Maxwell-Faraday equations, so they can be solved by classic finite-difference time-domain schemes, with additional steps due to the introduction of the auxiliary fields. Since the Maxwell-Gauss equation ∇ · D = 4πρ is not enforced by the code, we have to regularly perform divergence cleaning .
Particles We simulate pair plasma. In the 3 + 1 formalism, the Hamiltonian of a positron (or electron) of charge q = ±e, mass m e and 4-velocity u µ , moving in an electromagnetic field with 4-potential A µ , reads The particle's equations of motion are deduced from Eq. (6) and Hamilton's equations (Hughes et al. 1994;Dodin & Fisch 2010;Bacchini et al. 2018Bacchini et al. , 2019Parfrey et al. 2019): where F = q( D + (u/Γ) × B) is the Lorentz force, Γ = 1 + h jk u j u k is the FIDO-measured Lorentz factor of the particle, v is its 3-velocity with respect to the grid, and u/Γ its FIDO-measured 3-velocity. The electric current density J source term involved in Eq. (3) is determined by the contributions qv of each particle.
Photons In order to treat plasma supply self-consistently in our simulations, we use the radiative transfer algorithm introduced in Levinson & Cerutti (2018) and C20 (see the Supplemental Material therein for more details) in order to include IC scattering and γγ pair production. We include high-energy photons as a neutral third species that propagates along null geodesics. All particles evolve in a background soft radiation field, putatively emitted by the accretion flow, which makes the propagating medium opaque. We assume that this radiation field is static, homogeneous (with uniform density n 0 in any FIDO frame), isotropic, and mono-energetic (with energy ε 0 ). The opacity of the medium for all particles is parameterized by the fiducial optical depth τ 0 = n 0 σ T r g , where σ T is the Thomson cross-section. At every time step a lepton can scatter off a background soft photon to produce a high-energy photon by IC scattering, whereas a high-energy photon can annihilate with a soft photon to produce an e ± pair.
Two photons of energies ε and ε , colliding with an angle θ 0 , can only produce a pair provided ε ε (1 − cos θ 0 ) /2 ≥ (m e c 2 ) 2 (Gould & Schréder 1967). One of our main goals is to characterize high-energy emission by synthesizing light curves from our simulations, that can directly be compared to observations. To do so we need to keep track of the IC photons emitted below the pair-creation threshold (m e c 2 ) 2 /ε 0 . These photons are able to escape the soft photon field: they are responsible for the magnetospheric component of the high-energy emission received from Earth. However, these photons do not participate in the plasma dynamics. It is therefore critical to save the information they carry, as detailed in Sect. 3, before discarding them from the simulation. We do not include synchrotron and curvature radiation, which will be considered in a future work.
Parameters The simulations are axisymmetric: the particles move in 3D space, with all quantities being invariant by rotation around the spin axis of the black hole. The simulation domain is r ∈ [r min = 0.9 r h , r max = 10 r g ], θ ∈ [0, π], with r h = r g (1 + √ 1 − a 2 ) the radius of the event horizon. Fields are defined on a grid evenly spaced in log 10 (r) and θ. The spin parameter a is fixed at 0.99. We focus on the high optical depth regime, and run simulations with τ 0 = 30, 50, and 80.
The normalized magnetic field is defined byB 0 = r g (eB 0 /m e c 2 ), with B 0 the magnetic field strength at the horizon, whereas the normalized radiation field energy isε 0 = ε 0 /m e c 2 . We keptB 0 = 5 × 10 5 andε 0 = 5 × 10 −3 fixed in the simulation, in accordance with the criteria described in C20. This ensures that secondary pair-produced particles have high Lorentz factors, and that primary particles can be accelerated to energies above the pair-creation threshold. The magnetosphere is initially filled with high-energy photons, distributed uniformly and isotropically from r = r h up to r = 4 r g , with the energy ε 1 = 200 m e c 2 . From there the system takes about 100 r g /c to reach a steady state.
We use a damping layer at the outer boundary to absorb all outgoing electromagnetic waves in order to mimic an open boundary. Particles are removed if r ≤ r h or r ≥ r max . We performed our runs with a grid resolution 1536 (r) × 1024 (θ), with the requirement that we resolve the plasma skin depth everywhere, which was checked a posteriori. The fiducial density of the simulations is the Goldreich-Julian density n GJ = B 0 ω BH /(4πce), which is needed to screen the electric field induced by the rotation of the black hole. In this equation, ω BH = ca/(2r h ) is the black hole angular velocity.

Magnetic configuration
We simulate a generic magnetic configuration with large-scale poloidal field. This choice is the natural setup of the BZ mechanism (but see Parfrey et al. 2015 andMahlmann et al. 2020), and it is suggested by GRAVITY observations of the Galactic center (Gravity Collaboration et al. 2018). General relativistic magnetohydronamics (GRMHD) simulations of accretion flows also hint toward a paraboloidal geometry of the magnetic field lines (Komissarov & McKinney 2007;McKinney et al. 2012). The initial poloidal magnetic field in the magnetosphere is defined for θ ≤ π/2 by the following flux function (Tchekhovskoy et al. 2010) where r 0 and ν are free parameters. We chose r 0 = 10 r g and ν = 3 in our runs. The geometry of the initial magnetic field lines is shown in the upper panel of Fig. 1. In addition, in opposition to the work conducted in C20, this magnetic configuration naturally produces anti-parallel field lines at the equator because they are dragged toward the black hole during the simulation. This allows an equatorial current sheet to develop due to the discontinuity of the magnetic field. However, it is not known whether a configuration with an equatorial current sheet extending beyond the ergosphere is stable. In the Wald configuration studied by Parfrey et al. (2019) the current sheet did not extend beyond the ergosphere.
A magnetosphere with such an initial magnetic field quickly dies out after a few tens of r g /c (see Sect. 2.6). This occurs because the current sheet extends to the outer boundary of the box, which is endowed with open boundary conditions. Magnetic reconnection at the current sheet ejects plasmoids and magnetic flux from the simulation box. Too much energy and flux are lost by the simulation box, and the black hole almost completely expels the magnetic field lines threading it. Unlike pulsars, black holes do not have a conducting surface and cannot sustain a magnetic field. Therefore, we implemented a static and perfectly conducting disk as a boundary condition for the electromagnetic fields in the equatorial plane. This disk extends from its inner radius r = r in to the outer boundary of the box, for θ ∈ [π/2 − θ 0 , π/2 + θ 0 ], and we fixed r in = 6 r g and θ 0 = 0.02 in all simulations. The resulting setup is represented the lower panel of Fig. 1. Magnetic field lines crossing this disk are frozen in. The magnetic flux cannot escape the simulation box, which prevents our simulations from decaying entirely (see Fig. 5). We exclude the study of the magnetic linkage that can exist between Article number, page 3 of 12  Fig. 1: Schematic magnetic configuration. (a) Initial poloidal magnetic field lines, according to Eq. (9). (b) Magnetic configuration in steady state. Poloidal magnetic field lines are shown as red solid lines, except the last closed magnetic field line, which is in black. The equatorial current sheet (blue shaded area) is prone to the plasmoid instability. The conducting disk, represented by the gray shaded area, extends from r in to the outer edge of the simulation box. Two emitting zones are highlighted: the polar cap (low inclination with respect to the spin axis) and the current sheet.
the black hole and the disk, which is deferred to future work. We do not claim to simulate a realistic accretion disk, but rather to provide the physical conditions suitable to the study of the intrinsic behavior of the magnetosphere. The disk is merely included as a boundary condition for the fields. We are not interested in the zone surrounding the disk, and focus on the magnetosphere itself, that is, the zone enclosed by the field lines crossing the ergosphere. For this reason, in all subsequent figures we choose to leave the inner radius of the disk out of the represented domain. We also checked that there was no significant numerical diffusion, and hence no unphysical slippage of the field lines.

General features
We first describe the general features of our simulations before addressing the influence of magnetic field transport and the long-term evolution of the magnetosphere.
Structure The structure of the magnetosphere is shown in Fig. 2. The right panel shows H ϕ , which quantifies the poloidal current through a loop at constant (r, θ) according to Ampère's law ∇ × H = 4πJ/c. This toroidal field is nonzero on the field lines connected to the black hole, penetrating the ergosphere; therefore, a nonvanishing flux of energy and angular momentum can flow along those lines. The left panel shows the radial component of the current density. The electric current system is consistent with what is expected for a black hole magnetosphere in the force-free regime with a > 0 (Komissarov 2004). Our simulations have Ω · B > 0 in both hemispheres. In the upper hemisphere an electric field pointing toward the black hole is gravitationally induced by the frame-dragging of magnetic field lines. Negative poloidal currents are generated, which help screen the initial nonzero D · B, thus giving rise to a negative H ϕ . The situation is opposite in the lower hemisphere. By symmetry, H ϕ must vanish in the equatorial plane. The subsequent current sheet carries a positive electric current, closing the electric current system. This positive current flows along the separatrix, i.e., the last magnetic field line connected to the black hole, which defines the magnetospheric boundary.
Pair creation The equatorial current sheet is prone to plasmoid instability (Uzdensky et al. 2010), which mediates fast magnetic reconnection. Magnetic energy is dissipated and deposited into particles, leading to intense pair creation. Figure 3 shows a snapshot of the photon density above the pair-creation threshold and particle density, both in logarithmic scale.
We confirm that the mechanism described in C20 is still operating in this new configuration. Bursts of pair creation occur in an intermittent manner near the inner light surface 1 at intermediate latitudes. This fresh plasma mostly follows the magnetic field lines; therefore, it mainly flows close to the magnetospheric boundary. Inside the bursts the plasma density is marginally denser than the Goldreich-Julian density, and the 1 Light surfaces are the two surfaces separating subluminal and superluminal rotation for a point orbiting a Kerr black hole with angular velocity Ω (Komissarov 2004). The relevant inner light surface is defined by taking Ω as the velocity of the magnetic field lines. It is located within the ergosphere. outflowing plasma is highly magnetized. Pair creation is almost quenched near the rotation axis. We checked that in this zone the 4-current is null, although it is spacelike near the horizon at intermediate latitudes. In addition, the acceleration of particles in the X-points of the current sheet triggers pair creation and high-energy photon emission. The plasma density can reach 10 3 n GJ in the current sheet plasmoids.
Dynamics The magnetosphere displays an interesting dynamical phenomenon that is responsible for the replenishment of the magnetic field threading the black hole (see Fig. 4). Starting from an initial state similar to that shown in Fig. 2, plasma accumulates near the Y-point of the magnetosphere 2 . This plasma is supported by the magnetic pressure inside the magnetosphere. When the magnetosphere can no longer sustain the plasma, which roughly occurs when the particle energy density exceeds the magnetic energy density, a giant plasmoid forms and suddenly plunges into the black hole. This corresponds to the breakdown of the force-free approximation. The weakly magnetized plasma plunges due to the gravitational pull of the black hole, and works against the magnetic tension of field lines. As this giant plasmoid rushes inward, it pulls inward vertical magnetic field lines that were not crossing the event horizon initially. This replenishes the magnetic flux of the black hole. After the black hole swallows the giant plasmoid, the magnetosphere goes back to its initial state, until a new giant plasmoid is formed.

Long-term evolution
The outcome of the simulation is shown by the black and blue curves in Fig. 5, which represent the evolution of the magnetic flux Φ through the upper hemisphere of the event horizon with time. The magnetosphere experiences the dynamic cycles described in the previous section for about 300 r g /c, but the magnetic flux Φ decays secularly. It settles at a steady value after a time 500 r g /c. The steady state of the simulation resembles the Wald setup (Parfrey et al. 2019), as can be seen in Fig. 6. The field lines are much more vertical and, more importantly, the Y-point is located very close to the boundary of the ergosphere. In this steady state there are no more giant plasmoid accretion cycles. The current sheet is still disrupted by the tearing instability, so that small plasmoids fall toward the black hole. The magnetic flux escape by magnetic reconnection is exactly balanced by the supply of magnetic flux caused by the inflowing plasmoids pulling vertical field lines. Therefore, without external forcing, the only stable configuration for the magnetosphere is Wald-like. This is reminiscent of pulsar magnetospheres where the Y-point naturally migrates toward an equilibrium position at the light cylinder (Spitkovsky 2006). This configuration is close to force-free, except in the current sheet.
It should be noted that because the magnetic field strength has dropped significantly, the maximum Lorentz factor that particles can achieve is no longer much larger than the pair-creation threshold. Being slightly starved, large gaps can momentarily open up. This is merely an effect of our limited scale separation and does not affect the general conclusion. We also note that the final value for the magnetic flux does not depend on τ 0 in the range of parameters we have tested. If the opacity of the medium is high enough, the magnetosphere can reach a state close to force-free, irrespective of the plasma supply details.

Magnetic field transport
We are interested in maintaining the dynamic state and impeding magnetic field decay, since this variable state is promising for the prospect of high-energy flares. Therefore, we added the possibility of supplying magnetic flux to the central black hole in order to study the response of the magnetosphere either free or forced. To this end, we did not inject magnetic flux in the whole simulation box, but rather advected the frozen-in field lines that are initially crossing the perfectly conducting disk. We added a small toroidal electric field E acc = − (V 0 /c) × B only in the conducting disk. We ran another set of simulations of varying τ 0 , but this time with V 0 /c = 0.05. This setup can mimic inward magnetic flux transport in accretion flows (Lubow et al. 1994). The latter value of V 0 that we used is consistent with ideal MHD simulations of accretion disks (Jacquemin-Ide et al. 2020). By virtue of Faraday's law applied to a loop of radius r = r in in the equatorial plane, the magnetic flux through a surface enclosed by this loop must increase steadily for V 0 0 and remain constant for V 0 = 0. In other words, magnetic field lines that have been transported below r = r in at θ = π/2 must remain below r in from then on. This is why we place the inner boundary of the conducting disk at sufficient distance r in = 6 r g from the black hole. We choose not to run simulations with V 0 0 for as long as simulations with V 0 = 0 because the magnetic flux within r = r in would accumulate near the ergosphere.
In the presence of magnetic field line transport (V 0 0) the magnetosphere is able to remain in a dynamic state of periodic giant plasmoid accretion events (Fig. 5), and the inflow of magnetic flux compresses the magnetosphere and compensates the secular decay. The evolution of Φ in this dynamic state is represented in the upper panel in Fig. 8 for different optical depths, with the four blue dots representing successive snapshots relative to Fig. 4. As a magnetized giant plasmoid is swallowed by the black hole, the magnetic flux experiences a sharp rise. Between two successive giant plasmoid accretion events the magnetic flux decays almost exponentially with time due to magnetic reconnection. We observed that the characteristic decay time of Φ barely depends on τ 0 . On the other hand, the frequency of these accretion cycles is controlled by the fiducial optical depth τ 0 : it increases with increasing τ 0 . This occurs because mass loading at the Y-point is more efficient at high optical depth, which results in more frequent cycles of accretion. These cycles are illustrated in Fig. 7, which represents a spacetime diagram of the flux function A ϕ in the equatorial plane. They occur with a period of around 15 r g /c. The slow transport of magnetic field lines from the conducting disk to the black hole is also visible between 3 r g and 4 r g .
The lower panel in Fig. 8 shows the time evolution of the Poynting flux through the event horizon for the simulations with magnetic field transport. It is defined as where h denotes the determinant of the spatial 3-metric, the integration is performed over the event horizon B at r = r h , and is the radial component of the Poynting vector, i.e., the flux of electromagnetic energy-at-infinity (Komissarov 2004). We also observe sharp rises in the Poynting flux, synchronized with those in Φ. This comes as no surprise since the output power is expected to scale as Φ 2 if the Blandford-Znajek process is activated (Blandford & Znajek 1977;Tchekhovskoy et al. 2010). In the case of a pure split-monopole magnetosphere the total Poynting luminosity is L 0 = B 2 0 ω 2 BH /6. The measured luminosity is lower than this estimate because some flux is removed from the event horizon during an initial transient, due to the initial conditions not being an equilibrium state. This also explains why Φ is consistently below 2πr 2 g B 0 . The time-averaged luminosity corresponds to the total energy extracted from the black hole by the Blandford-Znajek process. In Sect. 3 all energy fluxes and light curves are normalized by this average luminosity L BZ .

Toy model for magnetic flux decay
The decay of the magnetic flux Φ through the event horizon, in the absence of any source term due to inflowing plasmoids, is a consequence of dissipation of magnetic energy by magnetic reconnection. We provide a toy model to account for the order of magnitude of the characteristic time T , assuming axisymmetry. The magnetic flux Φ can be expressed as where B + is the upper hemisphere of the event horizon. Faraday's law allows us to express the time derivative of Φ as the circulation of E along a loop of radius r h in the equatorial plane: dΦ dt = 2π∂ t A ϕ (r h , π/2) = −2πcE ϕ (r h , π/2).
In a purely axisymmetric and stationary magnetosphere, we have V in = V in ∂ θ toward the reconnection region. Just above and below the current sheet, the electric field reads E = − (V in /c) × B. At such high magnetizations, the outflow velocity is very close to c. We then define a global dimensionless reconnection rate R as where h is evaluated at r = r h and θ = π/2. Here V in is not measured by the FIDO, but with respect to the grid. We also assume that the configuration of field lines at the event horizon is close to a split monopole, so that B r (r = r h ) does not depend much on θ. The magnetic flux can be written as Φ = B r (r h )S, where S = 2π r 2 h + (ar g ) 2 is half the area of the event horizon (Bicak & Janis 1985). Ultimately, the evolution of the magnetic flux is governed by If the global reconnection rate is time-independent, the magnetic flux decreases exponentially with a characteristic decay time From Figs. 8 and 5 we measure the slope of the exponential decay, and obtain a reconnection rate R = 0.02 ± 0.002, corresponding to a decay time T 50 r g /c. We also measured the local reconnection rate using Eq. (14) and found values ranging from 0.02 to 0.04, which is consistent with the flux decay time. This model naturally explains why the characteristic decay time is the same in all simulations.
The local reconnection rate in collisionless relativistic reconnection has been determined by numerous numerical studies (e.g., Werner et al. 2018), with typical values between 10 −2 and 10 −1 . To compare with the decay time T , as measured by an observer at infinity, one needs to take into account gravitational dilation. Assuming the current sheet is roughly comoving (radially inward) with the Kerr-Schild FIDO at r = r h and θ = π/2, by definition of the lapse function α, the local reconnection rate can be estimated as R/α 1.66 R. Our value of R is consistent with measurements of the local collisionless reconnection rate, although slightly lower.
As mentioned in Sect. 2.2, we also ran simulations with no conducting disk. In these simulations the equatorial current sheet quickly extends across the whole box. The initial magnetic energy is quickly dissipated, whereas the mechanism previously analyzed cannot take place; there is no available vertical magnetic flux that could compensate for this decay. The simulation is exhausted of particles and dies out after a time on the order of T , which is solely determined by the reconnection rate.

Numerical method
As discussed in Sect. 2.1, we must save the information carried by photons that are below the pair-creation threshold. Given their initial positions, directions, and emission times, our goal is to reconstruct light curves for different viewing angles (with respect to the spin axis). This task has already been performed in flat spacetime (Cerutti et al. 2016) in the context of pulsar magnetospheres. Here this approach must be generalized to photons propagating in a curved spacetime. We neglect the gravitational influence of the black hole beyond a given r out , which we fix at 200 r g : for r ≥ r out , photons are considered to propagate in straight lines. We need to integrate the null geodesics of the photons from their emission points up to r = r out in order to compute their final directions and times of flight. Then, just as in flat spacetime, the light curve can be reconstructed if the directions of the outgoing photons are known at r = r out .
Keeping these photons in the simulation box and integrating their equations of motion with the PIC algorithm, even with a looser constraint on the time step, would be too demanding computationally. As it happens there is no need to solve the entire geodesic since the only relevant information is the initial and final coordinates (t, r, θ, ϕ) of the photons. Instead, we use the public ray-tracing code geokerr (Dexter & Agol 2009) . This code is optimized to integrate numerically null geodesics in the Kerr metric, and allows us to directly compute the final coordinates. If a photon produced by IC scattering is measured to be under the pair-creation threshold, its relevant information is dumped to a file that will be processed by geokerr before being discarded. All diagnostics can then be performed in post-processing. The light curve reconstruction procedure and the coupling between the two codes are detailed in Appendix A.

Results
We applied the previously outlined procedure to our three simulations with V 0 /c = 0.05, along with the highest opacity monopole simulation presented in C20 (which hadB 0 = 5 × 10 5 , ε 0 = 5 × 10 −3 and τ 0 = 30). The time resolution is ∆t = 0.0098r g /c, and the angular resolution is ∆α obs = π/24. We compute the energy flux per unit of solid angle. Some resulting light curves, normalized by the time-averaged BZ power of each simulation L BZ and by the solid angle ∆Ω obs = 2π sin α obs ∆α obs , are shown in Fig. 9. The monopole light curves do not depend much on the viewing angle, especially at intermediate latitudes. One example is shown in Fig. 9. This is expected since the monopole simulations show little structure in the orthoradial direction, and photons are mainly emitted radially by particles flowing along the magnetic field lines. We find a bolometric luminosity of L γ = dL γ / dΩ dΩ 0.04 L 0 0.04 L BZ . This is consistent with the dissipation rate of electromagnetic energy measured in C20, and confirms the hypothesis that the dissipated Poynting flux is mainly transferred to photons below the pair-creation threshold. Although the light curve shows signs of rapid variability, which is consistent with the small size of the gap, no flare can be detected. The incoherent process of pair creation along various magnetic field lines hinders the occurrence of large amplitude flares.
The situation is rather different for the disk simulations with external forcing (Fig. 9). These light curves show pronounced differences if viewed face-on (line of sight close to the spin axis, low α obs ) or edge-on (line of sight close to the equato-rial plane, α obs π/2). At low α obs they exhibit strong variability. During a flaring event the flux doubles within a rising time 2 r g /c. The periodicity of these flares is around 10 r g /c, in agreement with the periodicity of the giant plasmoid accretion cycles. Conversely, light curves observed at α obs π/2 are remarkably smooth, with no sign of variability at all. In order to understand these qualitative differences, we constructed two light curves, associated with the sites of emission of the photons. We distinguished the polar cap as the zone defined by θ ∈ [0 • , 60 • ] ∪ [120 • , 180 • ], and the current sheet, defined by θ ∈ [60 • , 120 • ] and r ∈ [r h , r in ]. Unsurprisingly, photons emitted in the polar cap mainly contribute to the emission at low viewing angles, whereas the emission at α obs π/2 is mainly due to the current sheet.
In the simulations with external forcing the average bolometric luminosity L γ ranges between 0.3 L BZ and 0.5 L BZ . Therefore, a very significant fraction of the BZ power is converted to IC luminosity. The radiative efficiency is much higher than in the monopole simulations. The total luminosity of photons emitted in the polar cap is around 5% of the BZ luminosity, a fraction that is similar to the monopole high-energy bolometric luminosity, although the polar cap emission is much more variable in these simulations. The current sheet and equatorial plasmoids emit 30% to 40% of the luminosity.
High variability should not be expected from magnetospheres observed edge-on. This stems from the fact that the formation of plasmoids in the current sheet, and the subsequent emission of high-energy photons, is inherently incoherent. Furthermore, photons emitted in this region travel along complex null geodesics that can have several turning points in θ. These geodesics are likely to differ significantly from simple radial rays; this adds to the decoherence that impacts current sheet emission, and erases any strong variability. On the other hand, photons emitted from the polar cap follow more direct geodesics toward the observer at infinity, such that the variability of the primary process is imprinted in the light curve. Therefore, the polar cap shows more pronounced variability in these simulations than in the monopole case. This indicates that the gap dynamics cannot be studied with no consideration for the global magnetospheric structure: the magnetospheric dynamics enhance the activity of the gap.
The power spectral densities of the light curves at α obs = 11.2 • for the simulations with V 0 /c = 0.05, computed using the stingray package (Huppenkothen et al. 2019), are shown in Fig. 10 (after logarithmic rebinning). At frequencies 0.1 c/r g f 2 c/r g , it is nicely described by a red-noise power law ∝ f −p , with p = 2.00 ± 0.13. A spectral break is visible around 0.1 c/r g . The spectral break frequency is consistent with the characteristic timescale associated with giant plasmoid accretion events. Resolving this spectral break observationally would require data acquisition for much longer than 10 r g /c (more than 4 days in the case of M87*). Beyond 2 c/r g , the power spectra are similar to white noise. Most of the power is distributed at lower frequencies: flux variations on long timescales dominate those on short timescales. The value of the index p is in agreement with that measured by Aharonian et al. (2007) from the AGN PKS 2155-304, although this measure may depend on whether the AGN is in a flaring state or not (H. E. S. S. Collaboration et al. 2017). We found that the value of p does not depend much on τ 0 . The characteristic value of the plasma frequency ν p = 8πe 2 n/Γ /m e /2π is 50 c/r g , and lies beyond the frequency range presented on the spectrum.

Discussion
We acknowledge that the results described in Sect. 2 would slightly differ in 3D since nonaxisymmetric modes would also allow the interchange of tenuous magnetospheric plasma with dense unmagnetized plasma at the Y-point. In particular, it is possible that Φ would not experience sharp and periodic peaks. Nonetheless, we believe this mechanism should still hold and allow the black hole to retain a significant magnetic flux and luminosity on a timescale longer than the characteristic reconnection decay time T .
In this paper we were interested in the pure magnetospheric response of a low-luminosity AGN. We find that a substantial fraction of the BZ electromagnetic luminosity is channeled into IC photons, especially if magnetic flux is supplied to the black hole by the accretion flow. This high radiative efficiency is the most salient feature of our simulations, especially in the presence of external forcing. We note that emission from equatorial latitudes is smoothed out, such that high-energy variability should primarily be expected from the polar caps. Variability from the polar caps is enhanced with respect to the monopole simulations. However, even with constant external forcing, the intrinsic activity of a steady-state black hole magnetosphere does not reproduce the most dramatic features of AGN flares: a fluxdoubling time below r g /c and an increase in the flux by a factor of at least 5, rather than 2 in our modeling. The variability seems to be well characterized by a red-noise power law down to 2 c/r g . Because of the nonaxisymmetric modes mentioned earlier, 3D simulations are likely to show even less variability. It should also be noted that in our simulations, the background radiation field is mono-energetic. Using a more realistic power law would have reduced the variability in the gap screening because the pair-creation threshold would not have been defined at a single photon energy. This makes it even less likely that realistic gamma-ray flares can be accommodated with our numerical setup.
We are able to quantify how fast the magnetic flux through the upper hemisphere of the event horizon decays, and find that external forcing is necessary in order to sustain a dynamic magnetospheric state. A free magnetosphere naturally tends toward a steady state similar to the Wald configuration. Recently, Ripperda et al. (2020) considered the possibility of magnetic reconnection powering infrared and X-ray flares in the Galactic center, and studied this scenario by resistive GRMHD simulations. They found that only in the magnetically arrested disk setup could there be a flaring state, during which plasmoids formed in an equatorial current sheet are heated to relativistic temperatures. This is consistent with our finding that in the absence of magnetic flux supply, the magnetosphere reaches a somewhat quiescent state and cannot produce flares.
Although the numerical setup used here can seem idealized, we note that the magnetic configuration that we have studied is actually the natural outcome, at the ergospheric scale, of any larger-scale configuration of an isolated magnetosphere (Komissarov & McKinney 2007). Therefore, we think that our findings concerning the generic magnetospheric response have broad applicability. If of magnetospheric origin, rather being a manifestation of the intrinsic variability due to the pair production mechanism, flares could be interpreted as the fast response of a black hole magnetosphere to a sudden change in the external parameters. This conclusion was also reached by Levinson & Cerutti (2018) and Kisaka et al. (2020) through radiative 1D GRPIC simulations. For example, a variation in the accretion rate would cause the density of soft photons to increase, leading to an increase in τ 0 . The velocity of magnetic field transport could also change, if a very magnetized plasma blob should accrete toward the magnetosphere. In that sense, black hole magnetospheres differ fundamentally from pulsar magnetospheres: pulsar activity is determined by parameters that are characteristic to the pulsar itself, and therefore remains quite stable. Another possibility is that flares result from a magnetosphere-disk interaction. This is suggested by the GRAVITY observation of a hot spot orbiting near the innermost circular orbit around Sgr A* (Gravity Collaboration et al. 2018).