Two-dimensional Particle-in-Cell simulations of axisymmetric black hole magnetospheres: angular dependence of the Blandford-Znajek flux

We examine the temporary evolution of axisymmetric magnetospheres around rapidly rotating black holes (BHs), by applying our two-dimensional particle-in-cell simulation code. Assuming a stellar-mass BH, we find that the created pairs fail to screen the electric field along the magnetic field, provided that the mass accretion rate is much small compared to the Eddington limit. Magnetic islands are created by reconnection near the equator and migrate toward the event horizon, expelling magnetic flux tubes from the BH vicinity during a large fraction of time. When the magnetic islands stick to the horizon due to redshift and virtually vanish, a strong magnetic field penetrates the horizon, enabling efficient extraction of energy from the BH. During this flaring phase, a BH gap appears around the inner light surface with a strong meridional return current toward the equator within the ergosphere. If the mass accretion rate is 0.025 percent of the Eddington limit, the BH's spin-down luminosity becomes 16-19 times greater than its analytical estimate during the flares, although its long-term average is only 6 percent of it. We demonstrate that the extracted energy flux concentrates along the magnetic field lines threading the horizon in the middle latitudes. It is implied that this meridional concentration of the Poynting flux may result in the formation of limb-brightened jets from low-accreting BH systems.


INTRODUCTION
The study of nonthermal plasmas in the vicinity of black holes (BHs) is astrophysically interesting in the context of collimated relativistic outflows observed from microquasars and active galactic nuclei. Such relativistic outflows, a.k.a. jets, are believed to be energized by rapidly rotating BHs that are immersed in a globally ordered magnetic field. In particular, when magnetic field (B) lines thread the event horizon, an electromotive force (EMF) is produced across the field lines in the same way as an unipolar inductor. In the direct vicinity of the horizon, this EMF induces a meridional current J , which exerts a counter torque on the BH via J ×B force (Blandford & Znajek 1977;Koide et al. 2002;McKinney et al. 2012). If such a current forms a closed circuit in a global BH magnetosphere, the extracted BH's rotational energy is carried away as the Poynting flux to be dissipated at some electric load located at large distances, e.g., as synchrotron emissions in the jet downstream.
In the case of low-luminosity radio galaxies, their flat spectrum radio emission with high brightness temperature is interpreted to originate in a synchrotronemitting jet (Jones et al. 1974a,b;Blandford & Königl 1979;Marscher 1983;Lobanov 1998). In the case of BH binaries, their flat spectrum radio emission during a hard/quiescent state is also considered to be from a jet (Hjellming & Johnston 1988;Stirling et al. 2001;Dhawan et al. 2000;Fender et al. 2004;Gallo et al. 2005). In such low accreting systems, their small plasma density results in a negligible turbulent diffusion, which prevents equatorial accretion to enter the jet-launching, polar regions in which horizon-penetrating magnetic field lines reside.
Although it is difficult to replenish jet materials by accretion in this way, photon-photon pair production is considered to be a viable mechanism as the source of jet plasmas. Gap models envisage large-scale regions in which a magnetic-field-aligned electric field accelerates charged leptons into ultra-relativistic energies, leading to a pair-production cascade of the gap-emitted γrays in a target, soft photon field (Beskin et al. 1992;Hirotani & Okamoto 1998;Levinson 2000;Levinson & Rieger 2011;Ptitsyna & Neronov 2016;Hirotani et al. , 2017Levinson & Segev 2017;Hirotani et al. 2018;Ford et al. 2018;Katsoulakos & Rieger 2020). Drizzle models, on the other hand, consider near-horizon, transient and small scale regions in which high-energy component of electrons produce MeV photons, which collide each other to produce pairs (Mościbrodzka et al. 2011;Wong et al. 2021). In the present paper, we focus on the former model, considering a jet-launching, low-accreting BH systems, which are realized in the center of low-luminosity radio galaxies or BH binaries in a hard/quiescent state.
When a BH is accreting at a highly sub-Eddington rate, such as in the case of M87* or BH binaries in a hard/quiescent state, plasmas become highly collisionless, which makes it impossible to justify the magnetohydrodynamic (MHD) approximations in the jet-launching regions (Hirotani et al. 2021, hereafter Paper I). Instead, kinetic method, or the particle-in-cell (PIC) method becomes appropriate (e.g., Nishikawa et al. (2021) for a recent review).
Assuming that the gap local physics does not affect the global structure, such as the magnetospheric currents, one-dimensional general relativistic PIC (GRPIC) simulations were performed by Levinson & Cerutti (2018); Chen & Yuan (2020); Kisaka et al. (2020). They consistently solved the radiative-transfer equation (or at least treated the inverse-Compton scatterings and pair production in a realistic way in the last case), together with the motion of the created pairs and the evolution of the electromagnetic fields.
To incorporate the back reaction on the global magnetosphere, however, we should proceed to twodimensional (2D) cases. In this context, Parfrey et al. (2019) first performed 2D GRPIC simulations of BH magnetospheres, assuming the pair injection rate is proportional to the strength of acceleration electric field, instead of solving the radiative transfer equation. Adopting an extremely small magnetic field strength, they demonstrated that the Penrose process contributes to the extraction of energy from a maximally rotating BH. Subsequently, Crinquand et al. (2020) examined 2D GR-PIC simulations, solving the radiative transfer equation, and assuming a fixed monopole magnetic field, which is, indeed, appropriate in a time-averaged sense in the horizon vicinity. They showed that the BH's rotational energy can be electromagnetically extracted via the Blandford-Znajek (BZ) process, and that a highly time-dependent spark gap opens near the inner light surface. In addition, Crinquand et al. (2021) coupled their 2D GRPIC code with a raytracing algorithm by post processing, and examined synthetic γ-ray light curves of the gap activity. Moreover, Bransgrove et al. (2021) applied 2D GRPIC and 3D MHD methods to a stellar-mass BH which is collapsing from a neutron star surrounded by plasma, and demonstrated that the 'no hair' theorem holds in the sense that the stress-energy tensor decays exponentially in time.
In addition to these works, in Paper I, we applied our 2D GRPIC code to stellar-mass BHs, without solving the radiative transfer equation. When solving the Maxwell equations, we solved only three components of the electromagnetic fields, namely the radial and meridional components of the electric field and the toroidal component of the magnetic field. Assuming a radial magnetic-field geometry near the BH, we demonstrated that the BH's rotational energy is preferentially extracted along the magnetic field lines threading the event horizon in the middle latitudes, namely between 60 • and 70 • (or 110 • and 120 • ) from the rotation axis.
Developing the 2D GRPIC method adopted in Paper I, we solve all the six components of the electromagnetic fields in the present paper. In the next section, we describe the basic equations in our GRPIC scheme. Then in § 3, we show that the main conclusion of Paper I -the middle-latitude concentration of the BZ flux -also holds when we solve all the six electromagnetic-field components, and demonstrate that the force-free approximation breaks down when plasmas are less efficiently supplied in the magnetosphere. We finally discuss an implication on the formation of a limbbrightened jet in § 5.

THE PARTICLE-IN-CELL (PIC) SCHEME
We formulate our axisymmetric, 2D GR PIC method in this section.
To avoid singular behaviours at sin θ = 0 (i.e., at the poles) due to the differential operator csc θ∂ θ in the Maxwell equations, we introduce a new meridional variable, y ≡ 1 − cos θ. Adopting this y coordinate, we obtain 1 sin θ Here, y = 0 (or y = 2) corresponds to the rotational axis in the uuper (or lower) hemisphere, and y = 1 denotes the equatorial plane.
In the radial direction, we adopt the so-called "tortoise coordinate" r * , In this coordinate, the event horizon corresponds to r * → −∞. Away from the BH (r M ), it tends to the standard radial coordinate, dr * /dr → 1.

Background electromagnetic fields
Throughout this paper, we assume that there exist stationary electromagnetic fields that are represented by the Wald solution (Wald 1974), which is realized when a ring current flows on the equator at a large enough distance from the BH. The amplitude of the ring current determines the magnetic field strength, B, near the BH. To specify B, we use its equipartition value, B eq , with the accreting plasmas. When accretion takes place as an advection-dominated accretion flow (ADAF) (Ichimaru 1977;Narayan & Yi 1994;Tchekhovskoy et al. 2011;Narayan et al. 2021), we obtain (Yuan & Narayan 2014), Here, the dimensionless accretion rateṁ is defined bẏ whereṀ denotes the mass accretion rate. The Eddington accretion rate is defined bẏ where L Edd denotes the Eddington luminosity. We adopt the conversion efficiency η eff = 0.1. In the present paper, we assume that the background, Wald solution has a field strength B so that its meridional average at r = 2M may match B eq . Since our code cannot simulate a plasma containing protons, we do no solve the normal plasmas consisting an ADAF. Instead, as the footpoint of a pair-dominated jet (Reynolds et al. 1996;Wardle et al. 1998;Hirotani 2005;Kino et al. 2014), we focus on the electron-positron pair plasmas that are created within the magnetosphere. Therefore, we useṁ merely to specify B eq (and hence B near the BH), and to specify the pair-supply rate,Ṅ ± ( § 2.6).

The Maxwell equations
In addition to the Wald solution described in the foregoing section, we consider additional electromagnetic fields produced by the electric currents flowing near the BH. Accordingly, the total electromagnetic fields are given by the superposition of the stationary Wald solution and these non-stationary fields. We describe how to solve the latter fields within our PIC scheme below.
When we solve the Maxwell equations, we should not solve for the Faraday tensor components that are defined in the Boyer-Lindquist coordinates, because numerical instabilities arise inside the ergosphere, in which a non-rotating observer becomes space-like. Instead, we can adopt a physical observer, such as the Zero-Angular-Momentum Observer (ZAMO) to avoid this ill behaviour. In this paper, instead of converting quantities into their ZAMO-measured values, we redefine the r and θ components of the electric field such that where F µν denotes a Faraday tensor component, and ω = 2M ar/A does the frame-dragging angular frequency. These transformations mimic the adoption of ZAMO; see also the arguments after equation (12) of Paper I for physical explanations of these variable transformations. In what follows, we put r g = 1 in equations for simplicity, and recover r g when appropriate.
Using these re-defined electric field (eqs. [9] & [10]), we find the following six time-dependent Maxwell equations: where and if we recover r g for clarity. The electric currents are defined by where the sigma symbol, m n=1 , denotes a summation of the electric currents carried by all the particle crossing the area where we count the current, m designates the number of macro particles in each subdomain, q n does the charge on the macro particle with identity number n, and δV does the invariant three-dimensional volume of each cell around the grid point. The fourvelocity components, u t , u r , u θ , and u ϕ satisfy the definition of the proper time, g µν u µ u ν = −1. Note that u r /u t = (1/c)dr/dt, u θ /u t = (r g /c)dθ/dt, and u ϕ /u t = (r g /c)dϕ/dt are all dimensionless.
It is worth noting that solving equations (11)-(16) for the PIC fields, E i and B i (i = 1, 2, 3), is equivalent with solving the Maxwell equations for the total electromagnetic fields, because the Maxwell equations are linear, and hence additive. For instance, writing the total fields as a summation of the (stationary) Wald fields and the (time-dependent) PIC fields, we find that all the terms containing the Wald fields vanish. Namely, temporal derivatives of the Wald fields vanish by definition. The curl of the Wald electric field vanishes by the Faraday's law. The curl of the Wald magnetic field vanishes by the Ampere's law, because we impose a current-free condition when obtaining the Wald solution. Note that the required equatorial, ring current at a large enough radius comes into the Wald solution only through the outer boundary condition, not through the (electric-current) source terms.
Let us also describe how to construct the electric currents. We divide the particle motion into the poloidal and toroidal components, adopting the area weighting (Villasenor & Buneman 1992) on the poloidal plane, and computing the toroidal component of the current from the azimuthal displacement of individual charged leptons, reflecting the area weighting. However, the area weighting in a curved spacetime results in an accumulation of small numerical errors in each time step, and eventually ends up with non-physical short-wavelength noise in the electromagnetic fields due to failure to satisfy the Gauss's law. Therefore, in the present paper, we apply the Poisson correction to the electric field (e.g., Langdon & Lasinski 1976).
Let us briefly describe the units used in the code. Electromagnetic field components E 1 , E 2 , E 3 , B 1 , B 2 , B 3 , are computed in the cgs gaussian unit. Note that all these six components are well-behaved at the horizon. The charge density ρ e is also measured in the cgs gaussian unit (i.e., statcoulomb cm −3 ). The current components, J 1 , J 2 , and J 3 are in statampere × (r g /c) unit. The macro particle's four-velocity components are dimensionless. Namely, u 0 = dt/dτ denotes the ratio of the elapsed coordinate time dt and the proper time dτ . Thus, u 0 became the Lorentz factor in the special relativistic limit. u 1 = dr/dτ denotes the radial velocity in r (not r * ) coordinate in dτ basis, u 2 = dθ/dτ does the meridional angular velocity, and u 3 = dϕ/dτ does the azimuthal angular velocity.
To solve these six Maxwell equations (eqs. [11]-[16]), we must impose boundary conditions. Along the northern and southern polar axes (i.e., at θ = 0 and θ = π), we impose At the outer boundary, we impose that the radial derivatives of E 1 , E 2 , E 3 , B 1 , B 2 , and B 3 vanish for simplicity. We tested several combinations of other radial boundary conditions, and confirmed that conclusions were unaffected.
At the inner boudary, we impose the "radiative boundary condition" so that the tangential components of the electromagnetic fields may look to a ZAMO like radiation propagating into the stretched horizon, or equivalently, all the components of the electromagnetic fields should be finite for a freely-falling observer ( § III C 4 of Thorne et al. 1986). In ZAMO's orthonormal basis, this condition reads (0, Bθ, Bφ) = (0, Eθ, Eφ ) × (1, 0, 0), which gives the following inner boundary conditions in our current notation.

Initial conditions
We start the PIC simulation from at time t = 0. Note that E i and B i (i = 1, 2, 3) describe the nonstationary fields that are to be superposed on the stationary Wald fields. We assume that the densities of electrons and positrons are 0.1n GJ at t = 0, where the Goldreich-Julian number density is defined by at each point. In each cell, we give 50 initial macro electrons and positrons (i.e., 100 in total) with equal charge weight, distributing their positions randomly in the cell. We assume that these macro particles are static on the poloidal plane and corotating with the space time at t = 0. That is, macro electrons and positrons are moving at the same toroidal velocity initially; accordingly, they do not carry any electric currents at t = 0.

Particle equation of motion
In a highly vacuum BH magnetosphere, charged leptons are deccelerated by the radiation-reaction forces. In the same manner as in Paper I, we include the radiation reaction force as a friction term in the equation of motion (EOM).
With such a friction term, the EOM can be expressed by (Chapter 17 of Jackson 1962) where dτ refers to the particle proper time, e represents the absolute value of the charge on the electron, and m does the mass of the electron. Let us briefly describe the boundary conditions on the motion of electrons and positrons. Due to the symmetry, we assume that the particles moving across the polar axis (at θ = 0 or π) will be reflected toward the equator with opposite meridional velocity. Both the inner and outer boundaries are treated as particle sinks. Thus, when particles move across these two radial boundaries, they are excluded from the simulation.

Plasma supply
In BH magnetospheres, pairs can be supplied via twophoton and/or one-photon (i.e., magnetic) pair production processes. In the present paper, we focus on the former process, and consider the collisions of MeV photons emitted via Bremsstrahlung from an ADAF. The collision rate, or equivalently, the pair supply rate (pairs per second per volume) is given bẏ where σ γγ denotes the total cross section of photonphoton pair production, and n γ does the MeV photon density. Adopting the Newtonian self-similar ADAF model (Mahadevan 1997), and assuming that the most energetic MeV photons are emitted within r = 4M , we obtain (Paper I) We randomly introduce a macro particle in each cell at every time step with probability 1/k create = 0.01; that is, particles are injected in each cell at every k create = 100 time steps on average. In this case, each created macro positron or electron has the electric charge where ∆ t denotes the interval of each time step, and ∆ V the invariant volume of each cell. Note that ∆ t ∆ V = √ −gdtdrdθdϕ = 2π √ −g∆ t ∆ r ∆ θ holds, where ∆ r and ∆ θ denote the intervals in Boyer-Linquist radial and meridional coordinates. In the present paper, r * and y coordinates are uniformly gridded; thus, both r and θ are gridded non-uniformly.
As the PIC simulation proceeds, the number of macro particles increases with t to saturate at a few hundred in each PIC cell on average. Here, the maximum value of the Courant number is set to be 0.5 in r * and y = 1 − cos θ coordinates. In total, there are about 3 × 10 8 macro particles in the entire simulation region.
To solve the temporal evolution of the electromagnetic fields and the particle distribution functions, we adopt a radial grid of 1200 uniform cells (in r * coordinate) between 1.01r H < r < r out = 25.0M , and 1120 uniform cells (in y coordinate) between 0 < y = 1 − cos θ < 2, which corresponds to 0 • < θ < 180 • . It is noteworthy that the plasma skin depth can not be resolved at r > 7M when the plasma density increases enough (particularly in the polar regions). Thus, we adopt only the results obtained in r < 7M as the appropriate PIC solution, interpreting that the solution in r > 7M merely gives the outer boundary condition at r = 7M . Alternatively, we could set some outer boundary conditions on the electromagnetic fields and the particle injection rate at r = 7M . Nevertheless, in the present paper, we infer the boundary conditions at r = 7M by solving the magnetosphere also in 7M < r < r out . By this treatment, the boundary conditions on the electromagnetic fields at r = r out little affect the solution in r < 7M . On the other hand, inward particle flux at r = 7M do depend on the condition imposed on the particle injection rate at r = r out . Thus, we assume no particle injection across r = r out , although particles freely escape outward across r = r out . Since the pair production rate rapidly decreases with radius (eq. [30]), particle inward flux at r = 7M is little affected by the outer boundary position, as long as it is located at r ≥ 25M .
We adopt the cgs Gaussian unit in the code. Thus, in the Maxwell equations (11)-(16), we measure the electric field E i (i = 1, 2, 3) in statvolts cm −1 , and the magentic field B i in gauss. Time and spatial variables are in dimensionless unit. For instance, t is measured in r g c −1 , and r is in r g ; y ≡ 1 − cos θ is dimensionless by definition. In the particle equations of motion (28), the four velocity u µ is measured in c, and the proper time τ in r g c −1 .
It is checked a posteriori that the invariant grid intervals resolve the skin depth at every point at any elapsed time, where the plasma frequency ω p is computed by the plasma density n ± and its mean Lorentz factor γ as (Levinson & Cerutti 2018) ω p = 4πe 2 n ± m γ .
In addition, we adopt a heavy electron mass of m = m p /20. Partly because we adopt this heavy electron/positron mass, and mainly because the leptons are ultra-relativistic ( γ 1), and because the plasma density is very small (n ± ∼ n GJ ), the plasma skin depth, l p , is resolved by the current grid intervals for stellar-mass BHs ( § 5.1).

NONSTATIONARY MAGNETOSPHERE
In this section, we apply the PIC method to a rapidly rotating stellar-mass BH with spin a = 0.9M and mass M = 10M . We adopt a small mass accretion ratė m = 2.5 × 10 −4 throughout this paper. In what follows, electromagnetic fields mean the total fields that are obtained by superposing the stationary Wald fields and the nonstationary PIC fields, unless explicitly mentioned to distinguish them.

Initial electromagnetic fields
Because of the frame dragging, there appears an electric field along the local magnetic field lines. Forṁ = 0.00025, we obtain a non-vanishing E · B/[B eq (2M )] 2 at t = 0 from the Wald solution, as presented in the left panel of figure 1. We also plot equi-A ϕ contours, which represent the magnetic field lines in the poloidal plane in a stationary and axisymmetric magnetosphere (as in the case of the Wald solution), as solid curves. Here, A ϕ denotes the magnetic flux function, and is the ϕ component of the vector potential A µ , which is related to the Faraday tensor by Since the contour interval is taken to be constant, the density of the solid curves shows the magnetic-field strength. In the right panel, we plot the magnetic field lines with smaller contour interval than the left panel, closing up the BH vicinity, < 2M , where (i.e., the abscissa) denotes the distance from the rotation axis measured in the Boyer-Lindquist r coordinate.
We define that the magnetic field direction so that it may point upward in both hemispheres. Thus, the yellow-red (or green-violet) region indicates inward electric field in the upper (or lower) hemisphere. This nonvanishing, magnetic-field-aligned electric field, will accelerate electrons outwards (or positrons inwards) in the higher latitudes, and accelerate electrons inwards (or positrons outwards) in the lower latitudes. The motion of such charges induces electric currents in the magnetosphere once the simulation begins at t = 0.

Electric fields and currents
As time elapses, the electric current carried by the charged leptons alter the electric field through the Ampere's law, and hence the magnetic field through the Faraday's law. Accordingly, the acceleration electric field, E · B/B, evolves to a qualitatively different configuration from the initial one. Figure 2 shows E · B/[B eq (2M )] 2 by the color images at four discrete elapsed times, t = 534M , 540M , 546M , and 552M . We also plot equi-A ϕ contours by black solid curves. It follows that a non-vanishing E · B changes with time, and becomes much greater than B eq 2 in some limited regions. At a small accretion rateṁ = 2.5 × 10 −4 , the ADAF cannot supply enough MeV photons that are capable of materializing as pairs, and cannot sustain the magnetosphere force-free. We also find that E · B peaks around the inner light surface ( fig. 2 of  at t = 540M , and that the magnetic field lines densely reside near the horizon at the same timing, as the top right panel shows. This point will be examined more closely in § 3.6 in relation to the energy extraction from the BH. In the top left panel, magnetic field lines do not seem to exist near the horizon; however, it merely shows that the magnetic field is week there. In figure 3, we plot the real charge density (color image) and the electric currents (red arrow) on the poloidal plane at the same four discrete timings. It shows that the charge density is comparable to the Goldreich-Julian value. Both the charge density and the current evolve within several dynamical timescales. It also follows from the top right panel that a return current is formed towards the equator within the ergosphere at t = 540M . We will examine this point more closely in § 3.6 in relation to the energy extraction from the BH.

Magnetic field
Let us next examine the magnetic field. First, we plot the ZAMO-measured radial component Br = F θϕ /( √ A sin θ) in figure 4. It follows that Br takes a large value in the higher and middle latitudes at t ∼ 540M and gradually decreases with time. By symmetry, Br changes sign between the upper and the lower hemispheres.
Second, we plot the ZAMO-measured meridional component, Bθ = ∆/AF ϕr / sin θ, in figure 5. Except for the horizon vicinity, Bθ takes negative values with smaller amplitude than Br. We can also find largeamplitude oscillations appearing in radial direction in the direct vicinity of the horizon. This is due to the redshift effect and will be discussed again in the final part of this subsection.
Third, we plot the ZAMO-measured toroidal component, Bφ = B 3 / √ ∆ in figure 6. For this horizontal component at the horizon, there appears a divergent factor, ∆ −1/2 , in the right-hand side, because the ZAMO becomes an unphysical observer at the horizon. (Note that we use ZAMO-measured quantities only for presentation purpose. We never adopt the ZAMO in any actual computation in our PIC scheme.) Neglecting this artificial divergence (i.e., the ill behaviour of ZAMO) at the horizon, we find |Bφ| < |Br| in most regions. It is also clear that Bφ generally becomes negative (or positive) in the upper (or lower) hemisphere, because the magnetic field lines tend to be swept back by rotation.
Let us briefly examine the evolution of magnetic-field lines. In figure 7, we present the equi-A ϕ contours near the horizon at four discrete timings as indicated. It is found that magnetic islands appear due to magnetic reconnection within the ergosphere at t ∼ 546M , and migrate towards the BH, being elongated along the horizon due to the gravitational redshift, as the right two panels show. Because of this effect, Bθ alternate the sign in the direct vicinity of the horizon, as figure 5 demonstrates. Note that it takes infinite time for magnetic islands to arrive the horizon in the Boyer-Lindquist coordinates. We thus adopt the tortoise coordinate, r * , to resolve the anti-parallel magnetic field lines accumulated in the direct vicinity of the horizon.

Magnetic reconnection near the horison
It is noteworthy that Crinquand et al. (2021) reported magnetic reconnections taking place on the equator, including outside the static limit (i.e., at r > 2M on the equator). In their simulation, their initial poloidal magnetic field (their eq. [9]) quickly dies out in a few tens of M = r g /c. Accordingly, relatively weak poloidal magnetic field outside the ergosphere allowed X and O points (where the magnetic field forms an X-like or O-like geometry in the poloidal plane) to arise at r > 3M (their fig. 4). On the other hand, in the present analysis, we superpose stationary Wald electromagnetic fields on nonstationary PIC fields, assuming a stationary ring current at some large-enough distance on the equator. As a result, the initial poloidal magnetic field does not die out and prevents X and O points to arise at r > 2M . That is, magnetic reconnection takes place efficiently only inside the static limit.
In figure 8, we plot the equi-A ϕ contours by white curves and the total pair number density, (n − +n + )/n GJ in color. It follows from the online animation that magnetic reconnections subsequently take place at the Xtype point, which is located at r < 2M and θ ∼ π/2 (i.e., inside the ergosphere near the equator). Pair plasmas Figure 1. Left: Initial distribution of E · B/Beq(2M ) 2 (color image) on the poloidal plane (r,θ), when the BH's mass and spin parameter is M = 10M and a = 0.9M , and the dimensionless mass accretion rate isṁ = 2.5 × 10 −4 . Black solid curves denote equi-Aϕ contours, which give magnetic field lines in a stationary and axisymmetric magnetosphere. The magnetic-field strength is chosen to be B(r = 2M ) = Beq(r = 2M ) after averaging over the spherical surface at r = 2M (see text). The electromagnetic fields are given by the Wald solution, which is stationary and axisymmetric. Both the abscissa and ordinate are measured in GM c −2 unit. The event horizon is located at r = 1.435M , and represented by the dashed semicircle. Right: Close up of the magnetic field lines in the BH vicinity.
are ejected from this reconnection region horizontally both rightward and leftward. Plasmas ejected rightward (i.e., into greater r) are compressed outside the ergosphere, 2.0M < r < 2.1M . On the other hand, plasmas ejected leftward (i.e., toward the horizon) are not efficiently compressed around the O-type point (located around 1.7M < r < 1.8M ), because the magnetic-field strength is weak there compared to those at r > 2M , as can be seen from the difference of contour intervals. Instead, in-falling plasmas are compressed in the direct vicinity of the horizon (in the Boyer-Lindquist coordinates); however, the resultant high-density regions (in red color) are hidden under the piled-up magnetic field lines, which appear as a bunch of white curves right above the horizon.

Particle energy distribution
Let us briefly browse the energy dependence of macro particles. In figure 9, we plot the distribution of u t = u 0 = dt/dτ as a function of the distance from the BH. In special relativity, u t refers to the Lorentz factor. However, in the present case, u t also contains the time dilation due to the gravitational redshift and the frame dragging at the particle's position. Nevertheless, for simplic-  . Snap shots of the real charge density, (n+ − n−)/nGJ (color image), and electric currents (red arrows) at four discrete elapsed times as indicated. The charge density is normalized by the Goldreich-Julian value, and is depicted in linear scale as shown by the vertical color bars. For instance, the green-violet (or yellow-red) regions show positive (or negative) dimensionless charge densities. For presentation purpose, the currents are measured in ZAMO and depicted in cgs unit (i.e., in statampere cm −2 ) in logarithmic scale; example arrows are depicted in the right-most panel. For presentation purpose, we uniformly divide the poloidal plane in X = r sin θ (abscissa) and Y = r cos θ (ordinate) directions, average the currents in each square cell, and plot the averaged current by a red arrow from the center of each cell. Both abscissa and ordinates are in rg = M = GM c −2 unit. The horizon resides at r = 1.435M .     ity, we represent u t as the 'Lorentz factor' in this figure. The color image shows the particle distribution function normalized by the Goldreich-Julian number density. For the details of this normalization, see § 3.5.2 of Paper I, where γ in equation (60) of Paper I corresponds to u t in the present notation.
We plot the electron (or positron) energy distributions at t = 540M as the top (or bottom) two panels. The left two panels show the distribution within the colatitude 39.69 • < θ < 41.42 • , and the right ones do within 58.72 • < θ < 60.02 • . It is found that the particle Lorentz factors are kept above 2 × 10 5 at each position. Because a strong acceleration electric field arises near the inner light surface at t ∼ 540M , the Lorentz factor distribution peaks within r < 1.8M for both electrons and positrons at this timing, t = 540M . If we integrate over the particles within the specified latitudinal range at each r at elapsed time t = 540M , we find that the positronic (or electronic) density peaks at ∼ 36n GJ (or ∼ 30n GJ ) with averaged Lorentz factor γ ∼ 8.3 × 10 6 (or ∼ 6.3 × 10 6 ) in the region 2.4M < r < 2.5M and 50 • < θ < 70 • .

The Blandford-Znajek flux
Now let us consider the BZ flux. The radial component of the BZ flux (i.e., the Poynting flux) become, In what follows, we normalize T em r t with its analytically inferred value, F ana BZ (r), where F ana BZ (r) ≡  (1 + a 2 /r 2 ) 2 − (∆a 2 /r 4 )(1 − z 2 )dz. The typical BZ luminosity is estimated by (Tchekhovskoy et al. 2010) where the total magnetic flux threading the horizon is given by and k ≈ 1/6π for a radial magnetic field near the horizon. The present magnetosphere is magnetically dominated, in the sense that the magnetic energy density dominates the particles' rest-mass energy densities. Accordingly, particle contribution is negligible when we discuss the energy flux.
In figure 10, we plot the normalized BZ flux, T em To see the averaged flux, we take a moving average with period 20M . Figure 11 shows the result for the upper hemisphere (left panel) and the lower hemisphere (right panel). It follows that the net flux is positive in both hemispheres, which means that the BH's rotational energy is extracted in the form of the Poynting flux, and that the solution is more or less symmetric between the upper and lower hemispheres. We can also find that the BZ flux concentrates in the middle latitude (red curves) rather than in the higher (green) or lower (blue) latitudes.
It is worth considering the relation between the BZflux variation and the magnetic-island evolution. As figure 7 indicates, when a magnetic island is formed near the equator within the ergosphere, it expels magnetic field lines outside. As time elapses, the magnetic island migrates inwards being elongated along the horizon, and eventually sticks to it and virtually disappear, as the right-most panel of figure 7 shows. After the magnetic islands virtually disappear, magnetic flux tubes return to the horizon vicinity, as the top-right and bottom-left panels of figure 2 show, or the left-most panel of figure 7 shows. During this phase, magnetic field becomes strong as the top-right panel of figure 4 shows, inducing a strong acceleration electric field near the inner light surface, as the top-right panel of figure 2 shows. It also follows from the top-right panel of figure 3 that a strong return current is formed at the same timing within the ergosphere towards the equator. This meridional current acts Lorentz forces on the magnetized plasmas pushing them into negative-energy orbits, and exerts a counter torque on the horizon. Accordingly, the BH's rotational energy is efficiently extracted.
On the other hand, when there exists a giant magnetic island within the ergosphere, the weak magnetic field near the horizon cannot facilitate the BZ process efficiently. Indeed, magnetic islands occupy a good fraction of the ergosphere in a large fraction of time. As a result, the BZ flux increases only sporadically, as figure 10 shows. The characteristic flaring period, ∼ 50M , is regulated by the time scale of reconnection taking place within the ergosphere.
It is worth noting that the BZ process is facilitated by virture of these return currents outside the horizon. Assuming a stationary and axisymmetric BH magnetosphere, Punsly (1996) pointed out the infeasibility of the BZ process in relation to the causality and the plasma's inertia. However, this criticism can be overcome by the formation of these strong, time-dependent return currents within the ergosphere.
Since we are not considering the pair production between the gap-emitted γ-rays and disk-emitted soft photons, we do not find the quasi-periodic gap activities through gap reopening and resultant pair cascade as reported in 1D GRPIC simulations with photon tracking (Kisaka et al. 2020;Chen & Yuan 2020).
To examine the angular dependence of the BZ flux, we plot the BZ flux as a function of the colatitude at five discrete elapsed times in the left panel of figure 12. At t = 540M , the BZ flux peaks at θ ∼ 60 • and 120 • . If we take a time average after t = 480M , we find that the BZ flux does peak in the middle latitude in 40 • < θ < 60 • (upper hemisphere) and in 120 • < θ < 140 • (lower hemisphere), as the right panel indicates.
Let us quickly take a look at the BZ luminosity, by integrating the BZ flux over the entire spherical surface at a fixed r. Whenṁ = 2.5 × 10 −4 , we can evaluate the magnetic-field strength at the horizon by its equipartition value, B eq = 1.53 × 10 6 G at r = 2M , which allows us to analytically infer the luminosity as L ana BZ = 4.16 × 10 34 ergs s −1 . Evaluating the simulated flux at r = 2.335M , and normalizing its integrated value with L ana BZ , we obtain the BZ luminosity as presented in figure 13. We find that the BZ luminosity flares every ∼ 50 dynamical time scales, as expected, and the peak values attain 16 − 19 times of the analytically inferred value.
The time-averaged BZ luminosity is found to be 2.63× 10 33 ergs s −1 , which corresponds to 6.2 % of the analytical value. Therefore, when the accretion rate is as small asṁ = 2.5 × 10 −4 , we can conclude that the spin-down luminosity of the BH exceeds 15 times of the analytical estimate during the flare although its long-term average is kept at only 6 % of it.

Power spectrum of the BZ flux
To look more closely at the time variability properties, let us Fourier-transform the BZ flux. For this purpose, we introduce a normalized power spectral density (PSD), P (f k ), of function f (t) = T em r t /F BZ ana , so that the normalization may be defined by where in the right-hand side f k = k/(N ∆ t ) denotes the frequency. The function (in this particular case, the BZ flux) f (t) is sampled at N points, which span a range of time T = (N − 1)∆ t . The sampling interval is ∆ t = 1.098×10 −3 GM c −3 ; thus, the Nyquist frequency is f c = 4.796 × 10 2 M −1 . We plot the PSD at θ = 45 • in figure 14. It follows that a quasi-periodic oscillation (QPO) appears around 0.02M −1 , and its higher harmonics appear around 0.04M −1 and 0.06M −1 . The fundamental frequency of the QPO, ∼ 0.02M −1 , means that there is a modulation of the amplitude with a period P QPO ∼ 50M , which corresponds to interval of the flares that can be seen in figures 10 and 11.

THE CASE OF GREATER ACCRETION RATE
Let us briefly consider a greater accretion rate,ṁ = 2.8 × 10 −4 . In figure 15, we present moving-averaged BZ fluxes at six discrete colatitudes for this case. Comparing with the left panel of figure 11, we find that the BZ flux shows flaring activity with an interval P QPO ≈ 49M , and that the BZ flux concentrates in the middle latitudes between 30 • and 60 • . Although the QPO frequency, P QPO , little depends onṁ, the peak of the normalized BZ flux decrease below 6, whereas it was between 8 and 11 whenṁ = 2.5 × 10 −4 . This is because the amplitude of the fluctuation reduces owing to the increased plasma density whenṁ (and hence the pair production rate) increases.     As an example, we present E · B/B eq (2M ) 2 at t = 300.00M (i.e., when the BZ flux peaks) in figure 16. Comparing with the top right panel of figure 2, which also shows E · B at a peak, we find that the magneticfield-aligned electric field decreases with increasingṁ. In another word, the magnetosphere becomes highly charge-starved when the pair production rate (eq. [30]) is small enough. Asṁ increases (i.e., as pair production increases), the magnetosphere approaches a force-free solution, in the sense that the electric field is more efficiently screened along the magnetic field lines. However, if we adoptṁ ≥ 3.0 × 10 −4 , the solutions become unstable from the horizon vicinity near the equator, because the variation of the compiled electromagnetic fields there becomes too sharp to be resolved even in the r * coordinates. This is a limitation which we may encounter in the Boyer-Lindquist coordinates. On the contrary, Parfrey et al. (2019) adopted the Kerr-Schild coordinates, which are regular at the horizon, and demonstrated that their PIC solutions approach the force-free solution when they assumed greater pair-production rates than ours.

DISCUSSION
To sum up, we simulated the evolution of a BH magnetosphere by a PIC scheme, when the electron-positron pair plasmas are steadily supplied externally. Provided that the mass accretion rate is as small as 0.025 % of the Eddington limit, the rotational energy of a black hole (BH) is electromagnetically extracted via the Blandford-Znajek process. For a ten-solar-mass BH with the spin parameter a = 0.9M , the extracted energy flux shows flaring activity with a period of 50 dynamical time scales, which is regulated by the magnetic reconnection within the ergosphere. During the flare, strong acceleration electric field appears around the inner light surface with a meridional return current towards the equator inside the static limit. The flare's energy flux concentrate along the magnetic field lines that thread the event horizon in the middle latitudes.

Plasma skin depth
Let us show that the plasma skin depth (eq. [32]) is resolved by the present grid interval. Normalizing the pair density by the GJ value, n ± = κn GJ , we obtain where γ 7 ≡ γ /10 7 and B 6 ≡ B/(10 6 G);B 0 ≡ eBr g /m e c 2 measures the magnetic field strength in dimensionless unit. For a non-rotating BH (i.e., a → 0), we find κ → ∞ so that κn GJ gives the actual pair density. In Parfrey et al. (2019); Crinquand et al. (2021); Bransgrove et al. (2021), they employedB 0 < 5 × 10 5 in their re-scaled formulation, which corresponds to B < 5 × 10 2 G for M = 10M and B < 5 × 10 −5 G for M = 10 9 M . Note that we have B eq ∼ 10 6ṁ −4 1/2 G for M = 10M , and B eq ∼ 10 2ṁ −4 1/2 G for M = 10 9 M , whereṁ −4 ≡ṁ/10 −4 . By virtue of such tiny magnetic field strengths, B B eq , the skin depth was resolved in their works even for supermassive BHs, M ∼ 10 9 M , despite γ 7 1 due to the smallness of B. In the present paper, on the other hand, we adopted a stellar-size BH mass M = 10M and heavy electron/positron mass m = m p /20. Evaluating the mean Lorentz factor, plasma density, and the magnetic field strength at each position at each time, we find γ 7 /(κB 6 ) > 0.001 at r < 3M and γ 7 /(κB 6 ) > 0.006 at r > 3M as the conservative lower limits, where M = r g . Therefore, it follows from the last line of equation (38) that the plasma skin depth is greater than 0.045M at r < 3M , which can be resolved by the present grid interval, < 0.008M there. Also at 3M < r < 7M , we find that l p > 0.11M can be resolved by the grid interval < 0.026M there.

Horizontal magnetic field lines at the exact horizon
As figure 7 shows, the magentic field lines become horizontal at the exact horizon when magnetized plasmas accrete, if we adopt the Boyer-Lindquist coordinates. Similar configuration of the lines of force also appear for the electric field when two oppositely charged particles are separated near the horizon ( fig. 17 of § II D 5 in Thorne et al. 1986). Nevertheless, if we introduce the 'stretched horizon' slightly above the true horizon and consider the physics only outside of the stretched horizon, we obtain horizon-penetrating lines of force at finite elapsed time. After an infinite time elapses, the stretched horizon eventually matches the true horizon. Figure 16. Magnetic-field-aligned electric field, E · B/Beq(2M ) 2 (color image) on the poloidal plane (r,θ), when the BH's mass and spin parameter is M = 10M and a = 0.9M , and the dimensionless mass accretion rate isṁ = 2.8 × 10 −4 . Black solid curves denote equi-Aϕ contours. Both the abscissa and ordinate are measured in GM c −2 unit. The event horizon is located at r = 1.435M .
Although the time-dependent magnetic lines of force is kept horizontal at the true horizon within a finite time, it does not mean that the BH's rotational energy can be extracted only at infinite elapsed time. This is because the magnetized plasmas fall onto the horizon with negative energies, which realizes outward Poynting flux continuously in space. In figure 17, we present the BZ flux measured at much smaller radius r = r H +0.1M = 1.535M . Note that the coordinate time t (i.e., the abscissa) corresponds to the proper time of a distant static observer. By virtue of the time-dependent in-falling motion of magnetized plasmas with negative energies (see e.g., Hirotani et al. 1992;McKinney et al. 2012, for the MHD Penrose proess), comparable amount of energy is carried outward at both r = 1.53M and r = 2.33M along individual magnetic flux tubes.

Formation of limb-brightened jets
Finally, let us briefly discuss how the middle-latitude concentration of the Poynting flux could result in the formation of limb-brightened jets, which were observed from Mrk 501 (Giroletti et al. 2004 In the present paper, we assumed that a nearly cylindrical magnetic field lines are sustained by a strong, equatorial ring current at a large distance from the BH, and superposed the electromagnetic fields created by the magnetospheric currents near the BH, using a PIC method. However, it is unrealistic to consider that such a ring current (that creates the Wald fields) would exist at much larger distance than the jet downstream region.  fig. 10 (i.e., forṁ = 2.5 × 10 −4 ) but measured at a smaller radius, r = 1.535M , and at six discrete colatitudes as labeled.
Thus, it is natural that the external Wald field is maintained only near the BH. Nevertheless, we could concatenate our PIC results (obtained near the BH) with the downstream region (far away from the BH), assuming e.g., that the time-averaged luminosity (i.e., Poynting flux times cross section) carried along each magnetic flux tube is constant as a function of the distance from the BH. Near the BH, we demonstrated that the Poynting flux is tiny in the higher latitudes. Thus, it is possible that a hollow jet is formed in the sense that the Poynting flux is small along the jet axis. In the subsequent paper, we will convert the Poynting flux as a function of the magnetic flux function, A ϕ , near the BH, and quantitatively argue the formation of limb-brightened jets at large distances from the BH.