Issue 
A&A
Volume 650, June 2021



Article Number  A163  
Number of page(s)  11  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202040158  
Published online  25 June 2021 
Synthetic gammaray light curves of Kerr black hole magnetospheric activity from particleincell simulations^{⋆}
^{1}
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
email: benjamin.crinquand@univgrenoblealpes.fr
^{2}
School of Mathematics, Trinity College Dublin, Dublin 2, Ireland
^{3}
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
Received:
17
December
2020
Accepted:
27
March
2021
Context. The origin of ultrarapid flares of very highenergy 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. Our aim is to bridge the gap between simulations and observations by synthesizing gammaray light curves in order to characterize the activity of a black hole magnetosphere, using kinetic simulations.
Methods. We performed global axisymmetric 2D generalrelativistic particleincell simulations of a Kerr black hole magnetosphere. We included a selfconsistent treatment of radiative processes and plasma supply, as well as a realistic magnetic configuration, with a largescale equatorial current sheet. We coupled our particleincell code with a raytracing algorithm in order to produce synthetic light curves.
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 quasisteady Waldlike configuration. The dissipated energy is mostly converted to gammaray photons. The light curves at low viewing angle (faceon) mainly trace the spark gap activity and exhibit high variability. On the other hand, no significant variability is found at high viewing angle (edgeon), 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 is detected in active galactic nuclei. These flares could result from the variation in parameters external to the black hole.
Key words: black hole physics / magnetic fields / acceleration of particles / plasmas / radiation mechanisms: nonthermal / methods: numerical
Movies are available at https://www.aanda.org
© B. Crinquand et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Groundbased Cherenkov telescopes have shown that active galactic nuclei (AGN) can be highly variable sources of very highenergy (VHE) emission (> 100 GeV) (Aharonian et al. 2007; Albert et al. 2007; Aleksić et al. 2014). Variability timescales can be shorter than the horizon lightcrossing 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 lightcrossing 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 lowluminosity 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 GoldreichJulian 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 BlandfordZnajek (BZ) process by providing the plasma necessary to establish a quasiforcefree 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 particleincell (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 selfconsistent 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 subhorizon variability.
However, in the case of isolated magnetospheres, more realistic configurations with a largescale 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, largescale 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 gammaray 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 postprocess 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 KerrSchild 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
In this equation α is the lapse function (the redshift of a FIDO with respect to the coordinate time t), β is the shift vector (the 3velocity of a FIDO with respect to the coordinate grid), whereas h_{ij} denotes the spatial 3metric 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 MaxwellAmpère and MaxwellFaraday equations, so they can be solved by classic finitedifference timedomain schemes, with additional steps due to the introduction of the auxiliary fields. Since the MaxwellGauss 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 4velocity u_{μ}, moving in an electromagnetic field with 4potential 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. 2018, 2019; Parfrey et al. 2019):
where F = q(D + (u/Γ)×B) is the Lorentz force, is the FIDOmeasured Lorentz factor of the particle, v is its 3velocity with respect to the grid, and u/Γ its FIDOmeasured 3velocity. 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 selfconsistently 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 highenergy 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 monoenergetic (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 crosssection. At every time step a lepton can scatter off a background soft photon to produce a highenergy photon by IC scattering, whereas a highenergy 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 highenergy 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 paircreation 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 highenergy 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_{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 by , with B_{0} 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 pairproduced particles have high Lorentz factors, and that primary particles can be accelerated to energies above the paircreation threshold. The magnetosphere is initially filled with highenergy 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 GoldreichJulian 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.
2.2. Magnetic configuration
We simulate a generic magnetic configuration with largescale 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)
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 antiparallel 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 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. 
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 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 longterm 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 forcefree 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 framedragging 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 steadystate quantities. Left panel: normalized electric current density J^{r}/ecn_{GJ}. Right panel: normalized toroidal component H_{φ}/B_{0}r_{g}. 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 r_{g}/c and 110 r_{g}/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 paircreation threshold and particle density, both in logarithmic scale.
Fig. 3.
Snapshot of the simulation with τ_{0} = 50 and V_{0}/c = 0.05 at t = 126 r_{g}/c. Left panel: logarithm of the normalized density of photons above the threshold, normalized by n_{GJ}. Poloidal magnetic field lines are represented by white solid lines. Right panel: logarithm of the normalized plasma density, normalized by n_{GJ}. 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^{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 GoldreichJulian density, and the outflowing plasma is highly magnetized. Pair creation is almost quenched near the rotation axis. We checked that in this zone the 4current is null, although it is spacelike near the horizon at intermediate latitudes. In addition, the acceleration of particles in the Xpoints of the current sheet triggers pair creation and highenergy 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 Ypoint 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 forcefree 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 σ = B^{2}/8πnm_{e}c^{2}Γ 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. Longterm 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 Ypoint 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 Waldlike. This is reminiscent of pulsar magnetospheres where the Ypoint naturally migrates toward an equilibrium position at the light cylinder (Spitkovsky 2006). This configuration is close to forcefree, except in the current sheet.
Fig. 5.
Evolution of the magnetic flux through the upper hemisphere of the event horizon for different V_{0} and τ_{0}. 
Fig. 6.
Snapshot of the simulation with τ_{0} = 50 and V_{0} = 0 at t = 611 r_{g}/c. Left panel: logarithm of the normalized photon density, normalized by n_{GJ}. Poloidal magnetic field lines are represented by white solid lines. Right panel: logarithm of the normalized plasma density, normalized by n_{GJ}. 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 paircreation 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 forcefree, 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 highenergy 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 frozenin 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 (JacqueminIde 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} ≠ 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 Ypoint 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 V_{0}/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
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 V_{0}/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 3metric, the integration is performed over the event horizon ℬ at r = r_{h}, and
is the radial component of the Poynting vector, i.e., the flux of electromagnetic energyatinfinity (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 BlandfordZnajek process is activated (Blandford & Znajek 1977; Tchekhovskoy et al. 2010). In the case of a pure splitmonopole 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 timeaveraged luminosity corresponds to the total energy extracted from the black hole by the BlandfordZnajek 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
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:
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
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
If the global reconnection rate is timeindependent, the magnetic flux decreases exponentially with a characteristic decay time
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^{−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 KerrSchild 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 paircreation 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 raytracing 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 paircreation 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 postprocessing. 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_{0}/c = 0.05, along with the highest opacity monopole simulation presented in C20 (which had , and τ_{0} = 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 timeaveraged 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 paircreation 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 V_{0}/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 timeaveraged value of L_{BZ} 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 faceon (line of sight close to the spin axis, low α_{obs}) or edgeon (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 highenergy 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 edgeon. This stems from the fact that the formation of plasmoids in the current sheet, and the subsequent emission of highenergy 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 rednoise 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 τ_{0}. 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 V_{0}/c = 0.05, as a function of the frequency f. The black dashed line shows the rednoise scaling f^{−2}, down to a spectral break located close to 0.1 c/r_{g}. 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 Ypoint. 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 lowluminosity 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 highenergy 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 steadystate 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 rednoise 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 monoenergetic. Using a more realistic power law would have reduced the variability in the gap screening because the paircreation threshold would not have been defined at a single photon energy. This makes it even less likely that realistic gammaray 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 Xray 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 largerscale 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 2018).
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.
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
 Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, Science, 325, 444 [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2006, Science, 314, 1424 [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2007, ApJ, 664, L71 [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, ApJ, 695, L40 [Google Scholar]
 Albert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ, 669, 862 [Google Scholar]
 Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2014, Science, 346, 1080 [Google Scholar]
 Aliu, E., Arlen, T., Aune, T., et al. 2012, ApJ, 746, 141 [Google Scholar]
 Bacchini, F., Ripperda, B., Chen, A. Y., & Sironi, L. 2018, ApJS, 237, 6 [Google Scholar]
 Bacchini, F., Ripperda, B., Porth, O., & Sironi, L. 2019, ApJS, 240, 40 [Google Scholar]
 Bicak, J., & Janis, V. 1985, MNRAS, 212, 899 [Google Scholar]
 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [Google Scholar]
 Carter, B. 1968, Phys. Rev., 174, 1559 [Google Scholar]
 Cerutti, B., Uzdensky, D. A., & Begelman, M. C. 2012, ApJ, 746, 148 [Google Scholar]
 Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147 [Google Scholar]
 Cerutti, B., Philippov, A., Parfrey, K., & Spitkovsky, A. 2015, MNRAS, 448, 606 [Google Scholar]
 Cerutti, B., Philippov, A. A., & Spitkovsky, A. 2016, MNRAS, 457, 2401 [Google Scholar]
 Chen, A. Y., & Yuan, Y. 2020, ApJ, 895, 121 [Google Scholar]
 Christie, I. M., Petropoulou, M., Sironi, L., & Giannios, D. 2019, MNRAS, 482, 65 [Google Scholar]
 Crinquand, B., Cerutti, B., Philippov, A., Parfrey, K., & Dubus, G. 2020, Phys. Rev. Lett., 124, 145101 [Google Scholar]
 Cunningham, C. T., & Bardeen, J. M. 1973, ApJ, 183, 237 [Google Scholar]
 Dexter, J., & Agol, E. 2009, ApJ, 696, 1616 [Google Scholar]
 Dodin, I. Y., & Fisch, N. J. 2010, Phys. Plasmas, 17, 112118 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019, ApJ, 875, L1 [Google Scholar]
 Font, J. A., Ibáñez, J. M., & Papadopoulos, P. 1999, MNRAS, 305, 920 [Google Scholar]
 Gould, R. J., & Schréder, G. P. 1967, Phys. Rev., 155, 1404 [Google Scholar]
 Gralla, S. E., Lupsasca, A., & Strominger, A. 2018, MNRAS, 475, 3829 [Google Scholar]
 GRAVITY Collaboration (Abuter, R., et al.) 2018, A&A, 618, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 H.E.S.S. Collaboration (Abdalla, H., et al.) 2017, A&A, 598, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hirotani, K., & Pu, H.Y. 2016, ApJ, 818, 50 [Google Scholar]
 Hughes, S. A., Keeton, C. R., Walker, P., et al. 1994, Phys. Rev. D, 49, 4004 [CrossRef] [Google Scholar]
 Huppenkothen, D., Bachetti, M., Stevens, A. L., et al. 2019, ApJ, 881, 39 [Google Scholar]
 JacqueminIde, J., Lesur, G., & Ferreira, J. 2021, A&A, 647, A192 [EDP Sciences] [Google Scholar]
 Kagan, D., Sironi, L., Cerutti, B., & Giannios, D. 2015, Space Sci. Rev., 191, 545 [NASA ADS] [CrossRef] [Google Scholar]
 Katsoulakos, G., & Rieger, F. M. 2018, ApJ, 852, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Kisaka, S., Levinson, A., & Toma, K. 2020, ApJ, 902, 80 [Google Scholar]
 Komissarov, S. S. 2004, MNRAS, 350, 427 [Google Scholar]
 Komissarov, S. S., & McKinney, J. C. 2007, MNRAS, 377, L49 [Google Scholar]
 Levinson, A. 2000, Phys. Rev. Lett., 85, 912 [Google Scholar]
 Levinson, A., & Cerutti, B. 2018, A&A, 616, A184 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Levinson, A., & Rieger, F. 2011, ApJ, 730, 123 [Google Scholar]
 Lubow, S. H., Papaloizou, J. C. B., & Pringle, J. E. 1994, MNRAS, 267, 235 [NASA ADS] [Google Scholar]
 MacDonald, D., & Thorne, K. S. 1982, MNRAS, 198, 345 [Google Scholar]
 Mahlmann, J. F., Levinson, A., & Aloy, M. A. 2020, MNRAS, 494, 4203 [Google Scholar]
 McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083 [Google Scholar]
 Mehlhaff, J. M., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2020, MNRAS, 498, 799 [CrossRef] [Google Scholar]
 Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69 [NASA ADS] [Google Scholar]
 Parfrey, K., Giannios, D., & Beloborodov, A. M. 2015, MNRAS, 446, L61 [Google Scholar]
 Parfrey, K., Philippov, A., & Cerutti, B. 2019, Phys. Rev. Lett., 122, 035101 [Google Scholar]
 Rieger, F. M., & Aharonian, F. 2012, Mod. Phys. Lett. A, 27, 1230030 [NASA ADS] [CrossRef] [Google Scholar]
 Ripperda, B., Bacchini, F., & Philippov, A. A. 2020, ApJ, 900, 100 [CrossRef] [Google Scholar]
 Spitkovsky, A. 2006, ApJ, 648, L51 [Google Scholar]
 Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2010, ApJ, 711, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Uzdensky, D. A., Loureiro, N. F., & Schekochihin, A. A. 2010, Phys. Rev. Lett., 105, 235002 [Google Scholar]
 Wald, R. 1974, Phys. Rev. D, 10, 1680 [Google Scholar]
 Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840 [Google Scholar]
 Yee, K. 1966, IEEE Trans. Antennas Propag., 14, 302 [Google Scholar]
Appendix A: Light curve reconstruction
Let us denote by p^{μ} the 4momentum of a photon with initial coordinates (t_{0}, r_{0}, θ_{0}, φ_{0}). 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 rescaled quantities l = L/ℰ and q^{2} = Q/ℰ^{2}. The photon trajectory is entirely determined by its initial position, the parameters l and q^{2}, 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_{0}, θ_{0}, t_{0}, φ_{0}, 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 3velocity can be reconstructed (Cunningham & Bardeen 1973; Gralla et al. 2018):
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 colatitude α 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
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)
The total time delay, from emission to collection by an observer, is t_{delay} = t_{f} − t_{0} + δ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_{0} + t_{delay} and in the direction θ_{f}:
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 KerrSchild coordinates (t, r, θ, φ), whereas geokerr uses spherical BoyerLindquist coordinates . The coordinate transformation from BoyerLindquist to KerrSchild reads (Komissarov 2004; Font et al. 1999)
along with and . This yields, in integrated form,
The t and φ coordinates of emission are transformed from KerrSchild to BoyerLindquist before being dumped by Zeltron. The covariant components x_{t}, x_{θ}, and x_{φ} of any oneform 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_{0} 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^{2} + a^{2} − al)^{2} − (r^{2} + a^{2} − 2r)(q^{2} + (a − l)^{2}) (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^{−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.
Movies
Movie 1 associated with Fig. 3 (density_fig_3) (Access here)
Movie 2 associated with Fig. 4 (sigma_fig_4) (Access here)
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 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. 

In the text 
Fig. 2.
Maps of the steadystate quantities. Left panel: normalized electric current density J^{r}/ecn_{GJ}. Right panel: normalized toroidal component H_{φ}/B_{0}r_{g}. 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 r_{g}/c and 110 r_{g}/c. 

In the text 
Fig. 3.
Snapshot of the simulation with τ_{0} = 50 and V_{0}/c = 0.05 at t = 126 r_{g}/c. Left panel: logarithm of the normalized density of photons above the threshold, normalized by n_{GJ}. Poloidal magnetic field lines are represented by white solid lines. Right panel: logarithm of the normalized plasma density, normalized by n_{GJ}. 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 σ = B^{2}/8πnm_{e}c^{2}Γ 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 V_{0} and τ_{0}. 

In the text 
Fig. 6.
Snapshot of the simulation with τ_{0} = 50 and V_{0} = 0 at t = 611 r_{g}/c. Left panel: logarithm of the normalized photon density, normalized by n_{GJ}. Poloidal magnetic field lines are represented by white solid lines. Right panel: logarithm of the normalized plasma density, normalized by n_{GJ}. 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 V_{0}/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 V_{0}/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 V_{0}/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 timeaveraged value of L_{BZ} 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 V_{0}/c = 0.05, as a function of the frequency f. The black dashed line shows the rednoise scaling f^{−2}, down to a spectral break located close to 0.1 c/r_{g}. The shaded areas indicate the error bars on the power spectra. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.