© B. Crinquand et al. 2021

1. 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 M 87 (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 M 87* VHE flares is comparable to the horizon light-crossing time t_g ≃ 0.4 days. M 87 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 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 M 87* 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. 1 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.

2. Magnetospheric activity

2.1. 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, 2015), 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² 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_ij 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)

(2)

(3)

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

(4)

(5)

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 (Cerutti et al. 2015).

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

(6)

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. 2018, 2019; Parfrey et al. 2019):

(7)

(8)

where F = q(D + (u/Γ)×B) is the Lorentz force, 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₀ in any FIDO frame), isotropic, and mono-energetic (with energy ε₀). The opacity of the medium for all particles is parameterized by the fiducial optical depth τ₀ = n₀σ_Tr_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 θ₀, can only produce a pair provided ε ε′(1−cos θ₀)/2 ≥ (m_ec²)² (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_ec²)²/ε₀. 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 the radius of the event horizon. Fields are defined on a grid evenly spaced in log₁₀(r) and θ. The spin parameter a is fixed at 0.99. We focus on the high optical depth regime, and run simulations with τ₀ = 30, 50, and 80.

The normalized magnetic field is defined by , with B₀ the magnetic field strength at the horizon, whereas the normalized radiation field energy is . We kept and 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 ε₁ = 200 m_ec². 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₀ω_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.

2.2. 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; Mahlmann et al. 2020), and it is suggested by GRAVITY observations of the Galactic center (GRAVITY Collaboration 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)

(9)

where r₀ and ν are free parameters. We chose r₀ = 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.

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 rin 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.

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 − θ₀, π/2 + θ₀], and we fixed r_in = 6 r_g and θ₀ = 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 Sect. 2.4). We exclude the study of the magnetic linkage that can exist between 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.

2.3. 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.

Fig. 2. Maps of the steady-state quantities. Left panel: normalized electric current density Jr/ecnGJ. Right panel: normalized toroidal component Hφ/B0rg. Poloidal magnetic field lines are represented by black solid lines. The dashed red line indicates the surface of the ergosphere. The quantities are averaged between t = 100 rg/c and 110 rg/c.

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.

Fig. 3. Snapshot of the simulation with τ0 = 50 and V0/c = 0.05 at t = 126 rg/c. Left panel: logarithm of the normalized density of photons above the threshold, normalized by nGJ. Poloidal magnetic field lines are represented by white solid lines. Right panel: logarithm of the normalized plasma density, normalized by nGJ. The dashed black line indicates the surface of the ergosphere, and the yellow solid line on the right panel shows the inner light surface. A movie is available online.

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¹ 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 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³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². 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.

Fig. 4. Snapshots of the logarithm of σ = B2/8πnmec2Γ from a simulation with τ0 = 50. The zones with σ ≃ 1 are in white. Poloidal magnetic field lines are represented by solid black lines. The dashed white line indicates the surface of the ergosphere. The full temporal evolution is available as an online movie.

2.4. 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.

Fig. 5. Evolution of the magnetic flux through the upper hemisphere of the event horizon for different V0 and τ0.

Fig. 6. Snapshot of the simulation with τ0 = 50 and V0 = 0 at t = 611 rg/c. Left panel: logarithm of the normalized photon density, normalized by nGJ. Poloidal magnetic field lines are represented by white solid lines. Right panel: logarithm of the normalized plasma density, normalized by nGJ. The dashed black line indicates the surface of the ergosphere.

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 τ₀ 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.

2.5. 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₀/c) × B only in the conducting disk. We ran another set of simulations of varying τ₀, but this time with V₀/c = 0.05. This setup can mimic inward magnetic flux transport in accretion flows (Lubow et al. 1994). The latter value of V₀ that we used is consistent with ideal MHD simulations of accretion disks (Jacquemin-Ide et al. 2021).

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 and remain constant for V₀ = 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 for as long as simulations with V₀ = 0 because the magnetic flux within r = r_in would accumulate near the ergosphere.

In the presence of magnetic field line transport (V₀ ≠ 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 τ₀. On the other hand, the frequency of these accretion cycles is controlled by the fiducial optical depth τ₀: it increases with increasing τ₀. 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.

Fig. 7. Spacetime diagram of Aφ in the equatorial plane for the simulation with τ0 = 50 and V0/c = 0.05. The black solid lines are the contours of Aφ, which represent poloidal magnetic field lines.

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

(10)

Fig. 8. Time evolution of the normalized magnetic flux through the upper hemisphere of the event horizon (upper panel) and of the normalized Poynting flux , integrated over the event horizon (lower panel), for the three simulations at V0/c = 0.05. The black disks, relative to the τ0 = 50 simulation, indicate the times of the snapshots relative to Fig. 4.

where h denotes the determinant of the spatial 3-metric, the integration is performed over the event horizon ℬ at r = r_h, and

(11)

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 Φ² 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 . 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 .

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⟩.

2.6. 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

(12)

where ℬ₊ 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:

(13)

In a purely axisymmetric and stationary magnetosphere, we have E_φ = 0 everywhere. However, this component arises in and near the current sheet because there is an inflow of plasma at velocity 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

(14)

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)𝒮, where is half the area of the event horizon (Bicak & Janis 1985). Ultimately, the evolution of the magnetic flux is governed by

(15)

If the global reconnection rate is time-independent, the magnetic flux decreases exponentially with a characteristic decay time

(16)

From Figs. 5 and 8 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⁻² and 10⁻¹. 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.

3. Synthetic light curves

3.1. 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.

3.2. Results

We applied the previously outlined procedure to our three simulations with V₀/c = 0.05, along with the highest opacity monopole simulation presented in C20 (which had , and τ₀ = 30). The time resolution is Δt = 0.0098 r_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.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.

Fig. 9. Light curve for the monopole simulation at τ0 = 30 (described in C20), normalized by , at αobs = 33.5° (black solid line), and light curves for the disk simulation with τ0 = 50 and V0/c = 0.05 at two different viewing angles: αobs = 11.2° (red dashed line) and αobs = 71.2° (blue solid line). These two light curves are normalized by the time-averaged value of LBZ shown in Fig. 8.

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 equatorial 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₀/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 M 87*). 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 2017). We found that the value of p does not depend much on τ₀. The characteristic value of the plasma frequency is 50 c/r_g, and lies beyond the frequency range presented on the spectrum.

Fig. 10. Power spectral density of the light curves at viewing angle αobs = 11.2°, for the simulations with V0/c = 0.05, as a function of the frequency f. The black dashed line shows the red-noise scaling f−2, down to a spectral break located close to 0.1 c/rg. The shaded areas indicate the error bars on the power spectra.

4. 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 flux-doubling 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 τ₀. 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 2018).

Movies

Movie 1 associated with Fig. 3 (density_fig_3) Access here

Movie 2 associated with Fig. 4 (sigma_fig_4) Access here

Acknowledgments

We thank Geoffroy Lesur and Amir Levinson for fruitful discussions that helped improve the article, and the anonymous referee for constructive and insightful suggestions. 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 TGCC and CINES under the allocation A0070407669 made by GENCI. A. P. acknowledges support by the National Science Foundation under Grant No. 2010145.

References

Appendix A: Light curve reconstruction

Let us denote by p^μ the 4-momentum of a photon with initial coordinates (t₀, r₀, θ₀, φ₀). The motion of the photon is specified by several constants of motion: the energy at infinity ℰ = −p_t, the component of the angular momentum relative to the spin axis of the black hole L = −p_φ, and the Carter constant (Carter 1968). We also define the re-scaled quantities l = L/ℰ and q² = Q/ℰ². The photon trajectory is entirely determined by its initial position, the parameters l and q², and two choices of sign in its initial radial and polar directions. It is independent of the photon energy. The PIC code Zeltron, at every time step, saves the quantities q, l, r₀, θ₀, t₀, φ₀, and E, as well as the initial sgn(ṙ) and sgn. The code geokerr is then used to integrate every saved geodesics up to r = r_out, and outputs the final photon coordinates (t_f, r_f, θ_f, φ_f). If the photon ends up at infinity, the spherical components of the normalized 3-velocity can be reconstructed (Cunningham & Bardeen 1973; Gralla et al. 2018):

(A.1)

(A.2)

(A.3)

The photons are then collected at infinity. Let us define a Cartesian basis (e_x, e_y, e_z), where e_z is aligned with the spin axis of the black hole and where φ is measured in the (e_x, e_y) plane with respect to e_x. On this basis, the final direction of a photon is defined by its co-latitude α and its azimuth ψ. Our simulations being axisymmetric, we assume that observers are solely characterized by a viewing angle α_obs, and are able to collect all photons whose final direction is α = α_obs, irrespective of ψ. The angle α is determined by

(A.4)

whereas ψ is determined and (with v^x = v ⋅ e_x and ). The eventual time delay for a photon at position (r_r, θ_f, φ_f), propagating toward an observer at (α_obs, ψ), is computed as (Cerutti et al. 2016)

(A.5)

The total time delay, from emission to collection by an observer, is t_delay = t_f − t₀ + δt. All in all, a photon arriving at r_out with an angle θ_f at a time t_f contributes to the luminous intensity per unit solid angle, at the time t₀ + t_delay and in the direction θ_f:

(A.6)

We do not keep any spectral information because of the low separation of scales: the fluxes we compute are integrated over all photon energies.

Zeltron uses spherical Kerr-Schild coordinates (t, r, θ, φ), whereas geokerr uses spherical Boyer-Lindquist coordinates . The coordinate transformation from Boyer-Lindquist to Kerr-Schild reads (Komissarov 2004; Font et al. 1999)

(A.7)

(A.8)

along with and . This yields, in integrated form,

(A.9)

(A.10)

The t and φ coordinates of emission are transformed from Kerr-Schild to Boyer-Lindquist before being dumped by Zeltron. The covariant components x_t, x_θ, and x_φ of any one-form x_μ are left unchanged by this transformation, so all constants of the motion have the same expression in both coordinate systems.

From the values of q, l, and r₀ for a given photon, it is possible to determine whether this photon will end up in the black hole or at infinity, by solving for the roots of the radial potential R(r) = (r² + a² − al)² − (r² + a² − 2r)(q² + (a − l)²) (Gralla et al. 2018). We perform this test so that only the photons that make it to infinity and actually contribute to the light curve are saved. For computational reasons, we only keep a fixed fraction of these photons (between 5% and 10%), which is sufficient to reconstruct the angularly resolved light curve.

We ran several sanity checks after coupling geokerr and Zeltron. While geokerr conserves the constants of the motion by construction, Zeltron does not. First we checked that for a large sample of photons whose trajectories are integrated by Zeltron, ℰ and Q are conserved to a relative accuracy of 10⁻¹⁰. Knowing that Zeltron’s integrator is reliable, we compared trajectories integrated by the two codes, and checked that they matched to a similar accuracy, regardless of the complexity of the geodesic (number of turning points in r or θ) or the number of points computed on the geodesic by geokerr.

All Figures

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 rin 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.

In the text

	Fig. 2. Maps of the steady-state quantities. Left panel: normalized electric current density Jr/ecnGJ. Right panel: normalized toroidal component Hφ/B0rg. Poloidal magnetic field lines are represented by black solid lines. The dashed red line indicates the surface of the ergosphere. The quantities are averaged between t = 100 rg/c and 110 rg/c.
In the text

Fig. 3. Snapshot of the simulation with τ0 = 50 and V0/c = 0.05 at t = 126 rg/c. Left panel: logarithm of the normalized density of photons above the threshold, normalized by nGJ. Poloidal magnetic field lines are represented by white solid lines. Right panel: logarithm of the normalized plasma density, normalized by nGJ. The dashed black line indicates the surface of the ergosphere, and the yellow solid line on the right panel shows the inner light surface. A movie is available online.

In the text

	Fig. 4. Snapshots of the logarithm of σ = B2/8πnmec2Γ from a simulation with τ0 = 50. The zones with σ ≃ 1 are in white. Poloidal magnetic field lines are represented by solid black lines. The dashed white line indicates the surface of the ergosphere. The full temporal evolution is available as an online movie.
In the text

	Fig. 5. Evolution of the magnetic flux through the upper hemisphere of the event horizon for different V0 and τ0.
In the text

	Fig. 6. Snapshot of the simulation with τ0 = 50 and V0 = 0 at t = 611 rg/c. Left panel: logarithm of the normalized photon density, normalized by nGJ. Poloidal magnetic field lines are represented by white solid lines. Right panel: logarithm of the normalized plasma density, normalized by nGJ. The dashed black line indicates the surface of the ergosphere.
In the text

	Fig. 7. Spacetime diagram of Aφ in the equatorial plane for the simulation with τ0 = 50 and V0/c = 0.05. The black solid lines are the contours of Aφ, which represent poloidal magnetic field lines.
In the text

	Fig. 8. Time evolution of the normalized magnetic flux through the upper hemisphere of the event horizon (upper panel) and of the normalized Poynting flux , integrated over the event horizon (lower panel), for the three simulations at V0/c = 0.05. The black disks, relative to the τ0 = 50 simulation, indicate the times of the snapshots relative to Fig. 4.
In the text

	Fig. 9. Light curve for the monopole simulation at τ0 = 30 (described in C20), normalized by , at αobs = 33.5° (black solid line), and light curves for the disk simulation with τ0 = 50 and V0/c = 0.05 at two different viewing angles: αobs = 11.2° (red dashed line) and αobs = 71.2° (blue solid line). These two light curves are normalized by the time-averaged value of LBZ shown in Fig. 8.
In the text

	Fig. 10. Power spectral density of the light curves at viewing angle αobs = 11.2°, for the simulations with V0/c = 0.05, as a function of the frequency f. The black dashed line shows the red-noise scaling f−2, down to a spectral break located close to 0.1 c/rg. The shaded areas indicate the error bars on the power spectra.
In the text