# Two-dimensional Particle-in-Cell simulations of axisymmetric black hole
magnetospheres

Kouichi Hirotani Institute of Astronomy and Astrophysics, Academia Sinica,
Taipei 10617, Taiwan, R.O.C.  Ruben Krasnopolsky Institute of Astronomy and
Astrophysics, Academia Sinica, Taipei 10617, Taiwan, R.O.C.  Hsien Shang （尚賢）
Institute of Astronomy and Astrophysics, Academia Sinica, Taipei 10617,
Taiwan, R.O.C.  Ken-ichi Nishikawa Department of Physics, Chemistry and
Mathematics, Alabama A&M University, Huntsville, AL 35811, USA  Michael Watson
Department of Physics, Fisk University, 1000 17th Ave North, Nashville, TN
37208, USA

###### Abstract

We investigate the temporal evolution of an axisymmetric magnetosphere around
a rapidly rotating, stellar-mass black hole, applying a two-dimensional
particle-in-cell simulation scheme. Adopting a homogeneous pair production,
and assuming that the mass accretion rate is much less than the Eddington
limit, we find that the black hole’s rotational energy is preferentially
extracted from the middle latitudes, and that this outward energy flux
exhibits an enhancement that lasts approximately 160 dynamical time scales. It
is demonstrated that the magnetohydrodynamic approximations cannot be
justified in such a magnetically-dominated magnetosphere, because the Ohm’s
law completely breaks down, and because the charge-separated electron-positron
plasmas are highly non-neutral. An implication is given regarding the
collimation of relativistic jets.

acceleration of particles — magnetic fields — methods: analytical — methods:
numerical stars: black holes

##  1 Introduction

It is commonly accepted that every active galaxy harbors a super-massive black
hole (BH) in its center (Ferrarese et al., 1996; Larkin & McLaughlin, 2016;
Kormendy & Ho, 2013). A likely mechanism for powering their relativistic jets
is the release of the rotational energy of the BHs (Blandford & Znajek, 1977),
which is referred to as the Blandford-Znajek (BZ) process, as demonstrated by
general-relativistic (GR) magnetohydrodynamic simulations (Koide et al., 2002;
McKinney et al., 2012). In the polar region of a BH magnetosphere, centrifugal
force prevents accretion toward the rotation axis and a high vacuum is
maintained (Hirose et al., 2004). In this nearly vacuum BH magnetosphere,
electrons and positrons ($e^{\pm}$’s) are supplied via the collisions of MeV
photons emitted from the equatorial, advection-dominated accretion flow (ADAF)
(Ichimaru, 1977; Narayan & Yi, 1994). Particularly, when the mass accretion
rate is much less than the Eddington rate, the collisions can no longer
sustain a force-free magnetosphere (Levinson & Rieger, 2011). In such a
charge-starved magnetosphere, an electric field appears along the magnetic
field lines. In this vacuum gap, a portion of the BZ flux is dissipated as
particle acceleration and radiation (Beskin et al., 1992; Hirotani & Okamoto,
1998; Neronov & Aharonian, 2007; Hirotani & Pu, 2016). Such a highly vacuum BH
magnetosphere has been investigated by the particle-in-cell (PIC) scheme one-
dimensionally along a radial magnetic field line (Levinson & Cerutti, 2018;
Chen et al., 2018; Chen & Yuan, 2020; Kisaka et al., 2020), and two-
dimensionally in the poloidal plane (Parfrey et al., 2019; Crinquand et al.,
2020). In the present paper, adopting a fixed magnetic field in the poloidal
plane, we perform a two-dimensional, axisymmetric PIC simulation around a
stellar-mass BH. We will examine the temporal evolution of the poloidal
components of the electric field, toroidal component of the magnetic field,
and the distribution functions of the electrons and positrons that are created
and accelerated within the BH mangetosphere.

After describing the background spacetime in § 2.1, we consider the initial
conditions of the PIC simulation in § 2. Then we present the results of our
2-D GR PIC simulation in § 3, demonstrating that the magnetohydrodynamic (MHD)
approximations totally break down in the vicinity (i.e., at the jet-launching
region) of rotating BHs that accrete plasmas at much lower rate than the
Eddington rate. In the present paper we focus on a set of astrophysical cases,
as described in § 3.5. Code testing cases, including comparison with exact
solutions (e.g., development of the electromagnetic field in a plasma-free
spacetime, or propagations of charged particles in a fixed electromagnetic
field) have been performed, but will be deferred to subsequent papers.
Finally, in § 4 we discuss the implication of astrophysical applications to
the very-long-baseline-interferometric (VLBI) observations of supermassive BHs
in the center of low-luminosity active galactic neclei.

##  2 Stationary magnetosphere

Let us begin with the stationary solution of the electromagnetic fields in a
magnetosphere of a rotating black hole. We will use this stationary solution
as the initial condition of time-dependent particle-in-cell (PIC) simulations.
Throughout the present paper, we consider two-dimentional structure of the BH
magnetosphere in which electrons and positrons are created and accelerated,
assuming axisymmetry with respect to the rotation axis of the BH.

###  2.1 Background geometry

Around a rotating, non-charged BH, the background geometry is described by the
Kerr metric (Bardeen, 1970). In the Boyer-Lindquist coordinates (Boyer &
Lindquist, 1967), the line element can be expressed as

$ds^{2}=g_{tt}dt^{2}+2g_{t\varphi}dtd\varphi+g_{\varphi\varphi}d\varphi^{2}+g_{rr}dr^{2}+g_{\theta\theta}d\theta^{2},$
(1)

where

$g_{tt}\equiv-\frac{\Delta-a^{2}\sin^{2}\theta}{\Sigma},\qquad
g_{t\varphi}\equiv-\frac{2Mar\sin^{2}\theta}{\Sigma},$ (2)
$g_{\varphi\varphi}\equiv\frac{A\sin^{2}\theta}{\Sigma},\qquad
g_{rr}\equiv\frac{\Sigma}{\Delta},\qquad g_{\theta\theta}\equiv\Sigma;$ (3)

$\Delta\equiv r^{2}-2Mr+a^{2}$, $\Sigma\equiv r^{2}+a^{2}\cos^{2}\theta$,
$A\equiv(r^{2}+a^{2})^{2}-\Delta a^{2}\sin^{2}\theta$. In equations (1)–(3),
we adopt the geometrized unit, putting $c=G=1$, where $c$ and $G$ denote the
speed of light and the gravitational constant, The horizon radius, $r_{\rm
H}\equiv M+\sqrt{M^{2}-a^{2}}$, is obtained by $\Delta=0$, where $M$
corresponds to the gravitational radius. The spin parameter becomes $a=M$ for
a maximally rotating BH, and $a=0$ for a non-rotating BH.

###  2.2 The Zero Angular Momentum Observer (ZAMO)

As a fiducial observer, let us introduce the Zero Angular Momentum Observer
(ZAMO), who is static in the poloidal plane ($r$,$\theta$) but rotates around
the BH at the same angular frequency as the spacetime frame dragging
frequency, $\omega\equiv-g_{t\varphi}/g_{\varphi\varphi}$. The tetrad of ZAMO
reads

$\tilde{\mbox{\boldmath$e$}}_{(\hat{t})}=\alpha^{-1}(\mbox{\boldmath$e$}_{(t)}+\omega\mbox{\boldmath$e$}_{(\varphi)}),$
(4)
$\tilde{\mbox{\boldmath$e$}}_{(\hat{r})}=\sqrt{\frac{\Delta}{\Sigma}}\mbox{\boldmath$e$}_{(r)},$
(5)
$\tilde{\mbox{\boldmath$e$}}_{(\hat{\theta})}=\frac{1}{\sqrt{\Sigma}}\mbox{\boldmath$e$}_{(\theta)},$
(6)
$\tilde{\mbox{\boldmath$e$}}_{(\hat{\varphi})}=\frac{1}{\sqrt{g_{\varphi\varphi}}}\mbox{\boldmath$e$}_{(\varphi)},$
(7)

where

$\alpha\equiv\frac{d\tau}{dt}=\frac{\rho_{\rm
w}}{\sqrt{g_{\varphi\varphi}}}=\sqrt{\frac{\Delta\Sigma}{A}}$ (8)

denotes the lapse; $d\tau$ refers to the ZAMO proper time. The tilde (
$\tilde{}$ ) represents a ZAMO-measured quantity and the caret ( $\hat{}$ )
does that the component is projected on an orthnormal basis. At the horizon we
have $\alpha=0$, while away from the BH we have $\alpha=1$. The coordinate
bases are defined as

$\mbox{\boldmath$e$}_{(t)}=\partial_{t},\quad\mbox{\boldmath$e$}_{(r)}=\partial_{r},\quad\mbox{\boldmath$e$}_{(\theta)}=\partial_{\theta},\quad\mbox{\boldmath$e$}_{(\varphi)}=\partial_{\varphi}.$
(9)

We use the ZAMO to solve the electromagnetic fields in a stationary
magnetosphere (§ 2.3), as well as to present the current densities (figs. 8 &
9) and the particle distribution functions (figs. 10 & 11) in a non-stationary
magnetosphere (§ 3).

###  2.3 Gauss’s and Biot-Savart laws

To describe the stationary electromagnetic field, we should simultaneously
solve the Gauss’s law and the Biot & Savart law. The expressions of these two
laws are derived in the Appendix B. For convenience, we replace the
independent variable $r$ with the so-called “tortoise coordinate”, $r_{\ast}$.
It is defined by

$\frac{dr_{\ast}}{dr}=\frac{r^{2}+a^{2}}{\Delta}.$ (10)

In this coordinate, the horizon corresponds to $r_{\ast}\rightarrow-\infty$.
At large distances, $dr_{\ast}/dr\rightarrow 1$.

What is more, to avoid the sigular behaviour due to the $\csc\theta$ factors
in the Gauss’s law (B17) and the Biot and Savart law (B19) at the poles (i.e.,
at $\theta=0$ and $\pi$), we introduce a new meridional variable, $y\equiv
1-\cos\theta$. Adopting this $y$ coordinate, we obtain

$\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}=\frac{\partial}{\partial
y}$ (11)

To avoid the change of the type of the second-order partial differential
equation (specifically, the Biot and Savart law) at the static limit, we
replace the electro-static potential $A_{t}$ with the ZAMO-measured value,
$A_{\hat{t}}$. The tetrad of ZAMO (eqs. [4]–[7]) gives

$A_{t}=\alpha A_{\hat{t}}-\omega A_{\varphi}$ (12)

and

$A_{\varphi}=\sqrt{g_{\varphi\varphi}}A_{\hat{\varphi}}.$ (13)

Using $A_{\hat{t}}$, $A_{\varphi}$, $r_{\ast}$ and $y$, the Biot and Savart
law (eq. [B19]) can be rewritten as

$\displaystyle-\hbox
to0.0pt{$\displaystyle\frac{\Delta\Sigma}{A}\frac{\partial^{2}A_{\varphi}}{\partial
r_{\ast}{}^{2}}+C_{1}\frac{\partial A_{\varphi}}{\partial
r_{\ast}}-\frac{\Delta^{2}\Sigma\sin^{2}\theta}{(r^{2}+a^{2})^{2}A}\frac{\partial^{2}A_{\varphi}}{\partial
y^{2}}$\hss}$ (14)
$\displaystyle+\left(\frac{\Delta}{r^{2}+a^{2}}\right)^{2}\frac{4Ma^{2}r\sin^{2}\theta\cos\theta}{\Sigma
A}\frac{\partial A_{\varphi}}{\partial y}+C_{0}A_{\varphi}$
$\displaystyle+\frac{2Mar\sin^{2}\theta}{\Sigma}\frac{\partial^{2}}{\partial
r_{\ast}{}^{2}}(\alpha A_{\hat{t}})+C_{2}\frac{\partial}{\partial
r_{\ast}}(\alpha A_{\hat{t}})$
$\displaystyle+\frac{2Mar\Delta\sin^{2}\theta}{(r^{2}+a^{2})^{2}\Sigma}\frac{\partial^{2}}{\partial
y^{2}}(\alpha A_{\hat{t}})$
$\displaystyle+\frac{4Mar\Delta\sin^{2}\theta\cos\theta}{(r^{2}+a^{2})\Sigma^{2}}\frac{\partial}{\partial
y}(\alpha A_{\hat{t}})$
$\displaystyle=4\pi\left(\frac{\Delta}{r^{2}+a^{2}}\right)^{2}\Sigma\sin^{2}\theta\cdot
J^{\varphi},$

where

$\displaystyle
C_{0}\equiv-\frac{2Mar\Delta\sin^{2}\theta}{(r^{2}+a^{2})^{2}\Sigma}\left\\{\right.\Delta(\partial_{r}{}^{2}\omega)$
(15)
$\displaystyle+\left[\frac{(r^{2}-a^{2}\cos^{2}\theta)\Delta}{r\Sigma}-\frac{2M(r^{2}-a^{2})}{r^{2}+a^{2}}\right](\partial_{r}\omega)$
$\displaystyle+\sin^{2}\theta(\partial_{y}{}^{2}\omega)$
$\displaystyle+\frac{2(r^{2}+a^{2})\cos\theta}{\Sigma}(\partial_{y}\omega)\left.\right\\},$
$\displaystyle C_{1}\equiv\frac{2\Delta}{(r^{2}+a^{2})^{2}\Sigma^{2}}$ (16)
$\displaystyle\times\left[\right.-(r-M)(r^{2}+a^{2})\Sigma+\frac{(r^{2}+a^{2})\Sigma}{A}C_{3}\qquad\qquad$
$\displaystyle\qquad+(\Delta-a^{2}\sin^{2}\theta)ra^{2}\sin^{2}\theta+\frac{4M^{2}a^{4}r^{3}\sin^{4}\theta}{A}$
$\displaystyle\qquad+2Mar\sin^{2}\theta(r^{2}+a^{2})\Sigma(\partial_{r}\omega)\left.\right],$
$\displaystyle
C_{2}\equiv\frac{2Ma\sin^{2}\theta}{(r^{2}+a^{2})^{2}\Sigma^{2}}$ (17)
$\displaystyle\times\left[\right.-r^{6}+a^{2}(-2+\cos^{2}\theta)r^{4}+4Ma^{2}\sin^{2}\theta
r^{3}$
$\displaystyle\qquad+a^{4}(1-2\sin^{2}\theta)r^{2}+a^{6}\cos^{2}\theta\left.\right],$
$\displaystyle C_{3}\equiv(r-M)A+a^{2}\sin^{2}\theta$
$\displaystyle\times\left[-r^{3}-3Mr^{2}+(2M^{2}-a^{2}\cos^{2}\theta)r+a^{2}M\cos^{2}\theta\right].$

It follows that all the four highest-order derivative terms of equation (14)
have definite signs; thus, the Biot-Savart law now become an elliptic type in
the entire simulation region, because we adopt a physical observer, ZAMO.

The Gauss’s law (B17) can also be re-written with respect to $\alpha
A_{\hat{t}}$ and $A_{\varphi}$. Multiplying
$\left[\Delta/(r^{2}+a^{2})\right]^{2}$ on both sides, we obtain

$\displaystyle\frac{A}{\Sigma^{2}}\frac{\partial^{2}(\alpha
A_{\hat{t}})}{\partial r_{\ast}^{2}}+D_{1}\frac{\partial(\alpha
A_{\hat{t}})}{\partial r_{\ast}}-\frac{\Delta
A}{(r^{2}+a^{2})\Sigma^{2}}\left(\partial_{r}\omega\right)\frac{\partial
A_{\varphi}}{\partial r_{\ast}}$ $\displaystyle+\frac{\Delta
A\sin^{2}\theta}{(r^{2}+a^{2})^{2}\Sigma^{2}}\frac{\partial^{2}(\alpha
A_{\hat{t}})}{\partial\xi^{2}}$
$\displaystyle+\frac{2\Delta\cos\theta}{(r^{2}+a^{2})^{2}\Sigma^{3}}\left[(r^{2}+a^{2})A-a^{2}\Delta\Sigma\sin^{2}\theta\right]\frac{\partial(\alpha
A_{\hat{t}})}{\partial\xi}$
$\displaystyle+\frac{2\Delta\sin^{2}\theta}{(r^{2}+a^{2})^{2}\Sigma^{2}}\left[-A(\partial_{\xi}\omega)+\Delta
a^{2}\cos\theta\cdot\omega\right]\frac{\partial A_{\varphi}}{\partial\xi}$
$\displaystyle+D_{0}A_{\varphi}=4\pi\rho\left(\frac{\Delta}{r^{2}+a^{2}}\right)^{2},$
(19)

where

$\displaystyle D_{0}\equiv\frac{\Delta
A}{(r^{2}+a^{2})^{2}\Sigma^{2}}\biggl{\\{}\biggr{.}-\Delta(\partial_{r}{}^{2}\omega)-\sin^{2}\theta(\partial_{\xi}{}^{2}\omega)$
(20) $\displaystyle+\left[2(r-M)-\Delta\partial_{r}\left(\ln\frac{\Delta
A}{\Sigma}\right)\right](\partial_{r}\omega)$
$\displaystyle-\frac{2\cos\theta}{\Sigma}\left[r^{2}+a^{2}-\frac{\Delta\Sigma}{A}a^{2}\sin^{2}\theta\right](\partial_{\xi}\omega)\biggl{.}\biggr{\\}}\quad$

and

$\displaystyle D_{1}\equiv-\frac{A}{(r^{2}+a^{2})\Sigma^{2}}$ (21)
$\displaystyle\times\left[2(r-M)-\Delta\partial_{r}\left(\ln\frac{(r^{2}+a^{2})A}{\Sigma}\right)\right]$

Note that $\partial^{2}A_{\varphi}/\partial r_{\ast}{}^{2}$ and
$\partial^{2}A_{\varphi}/\partial\xi^{2}$ terms vanish in equation (19).

###  2.4 Boundary conditions

We search for the stationary solution that satisfy the Gauss’s law (19) and
the Biot and Savart law (14). To this end, we must impose boundary conditions
$\alpha A_{\hat{t}}$ and $A_{\varphi}$, or equivalently on $A_{t}$ and
$A_{\varphi}$. In the present paper, we solve these two equations in the first
and fourth quadrants of the poloidal plane ($r$,$\theta$). The region is
bordered by
$\bullet$ the polar boundary at $\theta=0$ (north polar axis) and at
$\theta=\pi$ (south polar axis),
$\bullet$ the inner boundary at $r=r_{\rm H}$, and
$\bullet$ the outer oundary at $r=r_{\rm out}$, where $r_{\rm out}\gg r_{\rm
g}\equiv GMc^{-2}=M$.
In this subsection, we describe the boundary conditions at these three
boundaries.

At the northern and southern polar boundaries, we impose
$E_{\theta}=\partial_{\theta}A_{t}=0$, that is, a Neumann condition on
$A_{t}$. We also impose $B^{\theta}\propto F_{\varphi
r}=-\partial_{r}A_{\varphi}=0$. Thus, we put $A_{\varphi}=0$ at $\theta=0$ and
$\pi$ and measure the magnetic flux function $A_{\varphi}$ from the rotation
axis.

At the inner boundary, we impose
$\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}=0$ and $F_{\varphi r}=0$, where
$F^{r\theta}=0$ is assumed at $t=0$. For instance, in ZAMO, these conditions
constrain that both the radial component of the electric field and the
meridional component of the magnetic field vanish at the inner boundary.

At the outer boundary ($r=r_{\rm out}\gg M$), we impose
$\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}\propto(\partial_{r}A_{t})(\partial_{\theta}A_{\varphi})-(\partial_{\theta}A_{t})(\partial_{r}A_{\varphi})=0$
and $\partial_{r}A_{\varphi}=0$, the latter of which comes from the assumption
of a split-monopole magnetic field, $J^{\varphi}{}_{\rm eq}\propto r^{-4}$.
Thus, in the present paper we impose the Neumann condition,
$\partial_{r}A_{t}=0$. However, in general, if we impose the magnetic field
direction, $\partial_{r}A_{\varphi}/\partial_{\theta}A_{\varphi}$ (e.g., if we
adopt a paraboloidal magnetic field, see BZ77 for details),
$\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}=0$ constrains the direction of
the $A_{t}={\rm constant}$ surface at the outer boundary.

###  2.5 Disk toroidal current

We assume that the plasmas in an ADAF produce a toroidal current
$J^{\varphi}{}_{\rm eq}=C_{\rm eq}r^{-4}$ near the equator all the way to the
horizon within the colatitudes $87^{\circ}<\theta<93^{\circ}$; outside of this
equatorial region, $J^{\varphi}{}_{\rm eq}=0$ is assumed. For a slowly
rotating BH, this disk current produces a split-monopole magnetic field
(BZ77).

The normalization factor $C_{\rm eq}$ is adjusted so that the meridionally
averaged, ZAMO-measured radial magnetic field strength at $r=2M$,

$\langle
B^{r}(2M)\rangle\equiv\frac{\int_{0}^{\pi}\tilde{B}^{\hat{r}}(2M,\theta)\sqrt{A}\sin\theta
d\theta}{\int_{0}^{\pi}\sqrt{A}\sin\theta d\theta},$ (22)

may match a fraction of the equipartition field strength (Yuan & Narayan,
2014),

$B_{\rm eq}(r)=9.7\times
10^{7}\left(\frac{\dot{m}}{M_{1}}\right)^{1/2}\left(\frac{r}{2M}\right)^{-5/4}\mbox{
G},$ (23)

which is obtained if there is an equipartition between the magnetic field
energy density and the plasma internal energy density; the alpha viscous
parameter is assumed to be $0.3$. The dimensionless accretion rate $\dot{m}$
is defined by

$\dot{m}\equiv\frac{\dot{M}}{\dot{M}_{\rm Edd}},$ (24)

where $\dot{M}$ denotes the mass accretion rate. The Eddington accretion rate
is defined by

$\dot{M}_{\rm Edd}\equiv\frac{L_{\rm Edd}}{\eta_{\rm eff}c^{2}}=1.39\times
10^{19}M_{1}{\rm g\ s}^{-1},$ (25)

where $L_{\rm Edd}$ denotes the Eddington luminosity; the conversion
efficiency is assumed to be $\eta_{\rm eff}=0.1$.

In the present paper, we adopt a relatively small mass accretion rate,
$\dot{m}=2.25\times 10^{-4}$, and $\langle B^{r}(2M)\rangle=B_{\rm eq}(2M)$.

###  2.6 Stationary solution

In the present paper, we consider a ten solar-mass BH, $M_{1}\equiv
M/(10M_{\odot})=1$, and solve equations (14) and (19) iteratively in the
region $r_{\rm H}<r\leq 20M$ and $0\leq\theta\leq\pi$. The solved
electromagnetic fields are presented in figure 1. In each panel, thin black
curves show the equi-$A_{\varphi}$ contours, which indicate the poloidal
magnetic field lines for a distant static observer if $a=0$. We superpose
$\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}/[B_{\rm eq}(r=2M)]^{2}$, whose
values are indicated in the color code. The left and right panels correspond
to the cases of $a=0$ and $a=0.9M$, respectively.

In the left panel, $\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}=0$ holds
everywhere; thus, the background color is entirely white. It follows that the
poloidal magnetic field becomes radial for a slowly rotating BH when
$J^{\varphi}{}_{\rm eq}\propto r^{-4}$, which is consistent with the
analytical conclusion (Blandford & Znajek, 1977; McKinney & Narayan, 2007).
Because of $a=0$, there exists no magnetic-field-aligned electric field,
$E_{\parallel}$; thus, $A_{\varphi}$ is solved only from the Biot-Savart law.

However, as the BH spins up (i.e., if $a\neq 0$), the right panel shows that
there arises a non-vanishing $E_{\parallel}$ due to the relative rotational
motion of the magnetic field lines with respect to the space time. At the same
time, equi-$A_{\varphi}$ lines deviates from a radial shape, as the solid
curves indicate. For a rapidly rotating case, $a=0.9M$, the magnetic field
lines are laterally pushed toward the rotation axis in the ergosphere (Tanabe
& Nagataki, 2008; Tchekhovskoy et al., 2010). Note that the magnetosphere is
assumed to be vacuum in the present analysis, while it is force-free (i.e.,
plasma-filled) in Blandford & Znajek (1977); Tanabe & Nagataki (2008);
Tchekhovskoy et al. (2010). To examine how the magnetic field lines are
actually deformed for a physical observer, we adopt ZAMO and plot in figure 2
the radial ($\tilde{B}^{\hat{r}}$, left panel) and the meridional
($\tilde{B}^{\hat{\theta}}$, right panel) components of the magnetic field,
where their explicit expressions are given by equations (A12) and (A13). It
shows that $|\tilde{B}^{\hat{\theta}}|$ is kept below $|\tilde{B}^{\hat{r}}|$
except for the equatorial region where $\tilde{B}^{\hat{r}}$ vanishes by
symmetry. To compare with slowly rotating cases, we adopt the same power-law
in the disk current, $J^{\varphi}{}_{\rm eq}\propto r^{-4}$, for all the cases
of $a$. To avoid a substantial $|\tilde{B}^{\hat{\theta}}|$ in the lower-
latitude ergosphere, we could adopt another functional form of
$J^{\varphi}{}_{\rm eq}$ for $a\neq 0$; however, such a fine-tuning is out of
the scope of the present paper.

We assume a positive $J^{\varphi}{}_{\rm eq}$; thus, $F_{\theta\varphi}$
(i.e., the radial component of the magnetic field) is positive (or negative)
in the northern (or southern) hemisphere. Accordingly, a positive (or a
negative) sign of $\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}$ indicates that
the electric field points outwards in the northern (or southern) hemisphere.
Thus, as plasmas are created (via photon-photon collisions) near the horizon,
$r<2M$, electrons (or positrons) are accelerated inwards (or outwards) in both
hemispheres in this stationary solution. Such accelerated electrons and
positrons produce electric currents in the magnetosphere, whose poloidal
components modifies the poloidal electric field through the Ampere’s law. In
the next section, we will focus on the temporal evolution of the
electromagnetic fields and the particle distribution functions starting from
the initial conditions described in this section.

Figure 1:  Stationary equi-$A_{\varphi}$ contours (solid curves) and the
distribution of $\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}/B_{\rm
eq}(2M){}^{2}$ (color) on the poloidal plane ($r$,$\theta$). The equitorial
current density is assumed to depend on $r$ as $J^{\varphi}{}_{\rm eq}\propto
r^{-4}$. The ordinate represents the distance along the rotational axis of the
BH, while the abscissa does the distance $r\sin\theta$ from the rotation axis,
Both axes show the lengths in the Boyer-Lindquist coordinates, normalized by
the gravitational radius, $r_{\rm g}=GMc^{-2}=M$. The equatorial plane
corresponds to the ordinate of $0$. The black filled circle shows the BH. Left
panel shows the equi-$A_{\varphi}$ contours when $a=0$, in which case the
solid curves denote the magnetic field lines measured by a distant static
observer. For $a=0$, there arises no electric fields along the magnetic field
lines; thus, the background color is entirely white. Right panel shows the
case of $a=0.9M$. The amplitude of
$\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}$ increases with increasing $a$,
because the frame-dragging effect increases with increasing BH’s spin, $a$.
Figure 2:  Poloidal magnetic field measured in ZAMO at $t=0$ for $a=0.9M$. The
abscissa and ordinate are common with figure 1, but the BH vicinity is closed
up. Left and right panels show $B^{\hat{r}}({\rm ZAMO})$ and
$B^{\hat{\theta}}({\rm ZAMO})$, respectively, in Gauss.

##  3 The particle-in-cell (PIC) scheme

Let us look briefly at the collisionless nature of the plasmas in § 3.1,
before turning to a closer examination of the temporal evolution of the BH
magnetosphere in the rest of this section.

###  3.1 Collisionless plasmas

Denoting the density of a pair plasma with $n_{\pm}=\kappa n_{\rm GJ}$, we can
express the collision frequency as

$\nu_{\rm c}\sim\kappa n_{\rm GJ}\sigma c,$ (26)

where $\sigma$ refers to the collisional cross section, and

$n_{\rm GJ}\equiv\omega_{\rm H}B/(4\pi ce)$ (27)

denotes the Goldreich-Julian (GJ) number density, which is rotaionally
induced; $e$ denotes the charge on the electron. If the plasma density is
comparable to the GJ value, $\kappa$ becomes of the order of unity.

Let us evaluate the cross section by $\sigma\sim\pi l_{\rm c}{}^{2}$, where
$l_{\rm c}$ denotes the typical impact parameter. Equating the potential and
kinetic energies, we find $e^{2}/l_{\rm c}\sim(\gamma-1)m_{\rm e}c^{2}$. Then
combining it with equation (26), we obtain

$\nu_{\rm c}\sim\frac{\kappa\pi n_{\rm GJ}e^{4}}{\gamma^{2}m_{\rm
e}{}^{2}c^{3}},$ (28)

where $\gamma\gg 1$ is assumed. On the other hand, the gyro frequency is given
by

$\nu_{\rm B}=\frac{eB}{2\pi\gamma m_{\rm e}c}.$ (29)

We thus obtain

$\frac{\nu_{\rm c}}{\nu_{\rm
B}}\sim\frac{\pi}{4}\frac{\kappa}{\gamma}\frac{a}{M}\frac{r_{0}}{r_{\rm H}},$
(30)

where $r_{0}\equiv e^{2}/(m_{\rm e}c^{2})$ denotes the classical electron
radius. Since $r_{0}/r_{\rm H}\sim 10^{-19}M_{1}{}^{-1}$ holds, we can
conclude that the collision frequency is much less than the gyro frequency,
and that the assumption of collisionless plasmas in the PIC scheme is
justified. This conclusion solely comes from the fact that the GJ density
corresponds to a high vacuum in BH or pulsar magnetospheres.

It should be noted that the Ohm’s law, which is necessary to close the system
of equations in MHD, requires that many collisions should take place within a
single gyration. However, equation (30) shows this assumption cannot be
justified in a black hole magnetosphere unless the plasma density greatly
exceeds the GJ value. In the present paper, we thus construct the electric
current from the actual motion of charged particles, adopting the PIC scheme.
By this method, we can, for instance, incorporate the currents carried by the
drift motion of charged particles, as well as the anisotropic distribution of
particles’ momenta (e.g., along and perpendicular to the magnetic field
lines).

###  3.2 The Maxwell equations

In the present paper, we assume $F_{\varphi t}=0$ throughout the simulations.
Accordingly, together with $\partial_{\varphi}=0$ for all the quantities, we
find $\partial_{t}F_{\theta\varphi}=0$ and $\partial_{t}F_{\varphi r}=0$ from
the Faraday’s law. Thus, in the present paper, we treat both $B^{r}$ and
$B^{\theta}$ are constant in time (and hence unchanged from the initial
condition), and solve the temporal evolution of only $F_{rt}$, $F_{\theta t}$,
and $F^{r\theta}$. This assumption of $F_{\varphi t}=0$ is justifined as long
as the toroidal currents curried by the simulated particles are small compared
to the stationary, equatorial toroidal current that is carried by the
accreting plasmas.

Under this assumption, the Faraday and Ampere’s laws give

$\frac{\partial F^{r\theta}}{\partial
t}=-\frac{\Delta}{\Sigma^{2}}\left(\frac{\partial F_{\theta t}}{\partial
r}-\frac{\partial F_{rt}}{\partial\theta}\right),$ (31) $\frac{\partial
F_{rt}}{\partial
t}=\frac{\Sigma^{2}}{A}\left[\frac{1}{\sqrt{-g}}\frac{\partial}{\partial\theta}\left(\sqrt{-g}F^{r\theta}\right)-4\pi
J^{r}\right],$ (32) $\frac{\partial F_{\theta t}}{\partial
t}=\frac{\Sigma^{2}\Delta}{A}\left[-\frac{1}{\sqrt{-g}}\frac{\partial}{\partial
r}\left(\sqrt{-g}F^{r\theta}\right)-4\pi J^{\theta}\right].$ (33)

Replacing the $r$ derivatives with $r_{\ast}$ derivatives, and introducing the
following dependent variables,

$\displaystyle B$ $\displaystyle\equiv\sqrt{-g}F^{r\theta},$ $\displaystyle D$
$\displaystyle\equiv F_{rt},$ $\displaystyle E$ $\displaystyle\equiv F_{\theta
t}\sin\theta,$ (34)

all of which are well-behaved at the horizon, we can rewrite the three Maxwell
equations (31)–(33) into

$\frac{\partial B}{\partial t}=-c_{1}\frac{\partial E}{\partial
x}+c_{2}\frac{\partial D}{\partial y}$ (35) $\frac{\partial D}{\partial
t}=c_{3}\frac{\partial B}{\partial y}-4\pi\frac{\Sigma^{2}}{A}J^{r}$ (36)
$\frac{\partial E}{\partial t}=-c_{4}\frac{\partial B}{\partial
x}-4\pi\frac{\Sigma^{2}}{A}\Delta\sin\theta\cdot J^{\theta}$ (37)

where

$x\equiv r_{\ast},\qquad y\equiv 1-\cos\theta,$ (38)
$c_{1}\equiv\frac{r^{2}+a^{2}}{\Sigma},\qquad
c_{2}\equiv\frac{\Delta\sin^{2}\theta}{\Sigma},$
$c_{3}\equiv\frac{\Sigma}{A},\qquad c_{4}\equiv\frac{\Sigma(r^{2}+a^{2})}{A}.$
(39)

All the coefficients $c_{1}$, $c_{2}$, $c_{3}$, and $c_{4}$ are positive
definite. Both $c_{1}$ and $c_{4}$ are close to unity in the entire region.
However, $c_{2}\ll 1$ holds near the horizon or near the pole but
$c_{2}\rightarrow 1$ at $r\gg M$ in the lower latitudes. We have $c_{3}\ll 1$
at $r\gg M$ and $c_{3}\approx 0.25$ near the horizon. Note that the singular
behaviour (i.e., the polynomial pole) at $\theta=0$ in equations (32) and (33)
is eliminated by introducing a new meridional coordinate, $y$.

To solve these three Maxwell equations (eqs. [35]–[37]), we must impose
appropriate boundary conditions that are consistent with the initial,
stationary state (§2.4). Along the northern and southern polar axes (i.e., at
$\theta=0$ and $\theta=\pi$), we impose

$B=0,\quad\frac{\partial D}{\partial y}=0,\quad E=0.$ (40)

At the inner boundary, $x=r_{\ast}=-\infty$, we impose

$\displaystyle\frac{\partial B}{\partial x}=0,$
$\displaystyle\quad\tilde{E_{\hat{r}}}\propto D+\omega F_{r\varphi}=0,$ (41)
$\displaystyle\frac{\partial E}{\partial x}=\left(\frac{\partial E}{\partial
x}\right)_{t=0},$

where the quantity within $()_{t=0}$ indicates the initial value at $t=0$. At
the outer boundary, $x=r_{\ast}=r_{\rm out}$, we impose

$\frac{\partial B}{\partial x}=0,\quad D=0,\quad\frac{\partial E}{\partial
x}=\left(\frac{\partial E}{\partial x}\right)_{t=0}.$ (42)

Note that we set $B=0$ in the entire region at $t=0$.

###  3.3 Particle equation of motion

In a highly vacuum BH magnetosphere, charged leptons are deccelerated by the
radiation-reaction forces. To include these forces from the first principles,
we must adopt tiny time steps and consider the force on one part of the charge
by the fields of another part, taking account of retardation within the
particle itself. However, in actual simulations, it is unrealistic to adopt
such tiny timesteps. Thus, as a compromize, 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 as (Thorne & MacDonald,
1982; Hughes et al., 1994; Bacchini et al., 2018)

$\displaystyle\frac{du_{i}}{dt}=$ $\displaystyle-$ $\displaystyle\alpha
u^{t}\partial_{i}\alpha+u_{k}\partial_{i}\beta^{k}-\frac{1}{2}\frac{u_{j}u_{k}}{u^{t}}\partial_{i}g^{jk}$
(43) $\displaystyle+$
$\displaystyle\frac{q}{m}\left(F_{it}+F_{ij}\frac{u^{j}}{u^{t}}\right)+\frac{(F_{\rm
rad})_{i}}{u^{t}},$

where $\alpha$ is defined by equation (8), $\beta^{r}=\beta^{\theta}=0$,
$\beta^{\varphi}=g_{t\varphi}/g_{\varphi\varphi}$, $q/m$ refers to the charge-
to-mass ratio, and

$u^{t}=\frac{\sqrt{1+g^{jk}u_{j}u_{k}}}{\alpha}.$ (44)

The Latin indices run over 1, 2, 3. We can evaluate the radiation-reaction
force by the covariant form, (§ 17.8 of Jackson, 1962)

$F_{\rm
rad}^{j}=\frac{2}{3}\frac{r_{0}}{c}\left(\frac{d^{2}u^{j}}{d\lambda^{2}}+u^{j}\frac{du^{\nu}}{d\lambda}\frac{du_{\nu}}{d\lambda}\right),$
(45)

where $r_{0}$ denotes the classical electron radius, $\lambda$ does the
particle’s proper time, and the Greek indices run over 0, 1, 2, 3. This
radiation reaction force includes the effects of any kind of acceleration
acting on the particles. For instance, photon emissions as a result of the
acceleration in an electro-magnetic field (e.g., via the synchro-curvature
process) and that in a gravitational field, are included (Appendix D).
However, this radiation-reaction force does not include the effect of inverse-
Compton scatterings, which should be considered separately in future works.

Definition of the four velocity $u^{\mu}$ gives

$\frac{dr}{dt}=\frac{u^{r}}{u^{t}},$ (46)
$\frac{d\theta}{dt}=\frac{u^{\theta}}{u^{t}},$ (47)

and

$\frac{d\varphi}{dt}=\frac{u^{\varphi}}{u^{t}}.$ (48)

For presentation purpose, we can convert $u_{i}=dx_{i}/d\lambda$ in terms of
the ZAMO-measured spatial velocity, $u_{\hat{j}}$, as described in Appendix A.

We integrate both equations (43) and (46)–(48) in the phase space with the
global time variable $t$, which corresponds to the proper time of a distant
static observer (i.e., us).

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 $\theta=0$ or $\pi$) will be reflected equator-ward with
opposite meridional velocity. Both the inner and outer boundaries are treated
as particle sinks. Thus, when particles move accross these two radial
boundaries, they are excluded from the simulation.

###  3.4 Plasma supply

In BH magnetospheres, pairs can be supplied via two-photon 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 from
the equatorial ADAF via Bremsstrahlung. In subsequent papers, we will also
consider the collisions between the gap-emitted (inverse-Compton) photons and
the ADAF-emitted (synchrotron) photons.

The pair supply rate (pairs per second per volume) is given by

$\dot{N}_{\pm}=c\sigma_{\gamma\gamma}n_{\gamma}{}^{2},$ (49)

where $\sigma_{\gamma\gamma}$ denotes the total cross section of photon-photon
pair production, and $n_{\gamma}$ 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 (Appendix C)

$\dot{N}_{\pm}\approx 1.0\times
10^{24}\dot{m}^{4}M_{1}{}^{-2}\max\left[\left(\frac{r}{4M}\right)^{-4},1\right]$
(50)

in cgs unit (i.e., in $\mbox{pairs s}^{-1}\mbox{ cm}^{-3}$).

We randomly introduce a macro particle in each cell at every time step with
probability $1/k_{\rm create}=0.1$; that is, particles are injected in each
cell at every $k_{\rm create}=10$ time steps on average. In this case, each
created macro positron or electron has the electric charge

$q_{i}=\pm e\dot{N}_{\pm}k_{\rm create}\Delta_{t}\Delta_{V},$ (51)

where $\Delta_{t}$ denotes the interval of each time step, and $\Delta_{V}$
the invariant volume of each cell. Note that
$\Delta_{t}\Delta_{V}=\sqrt{-g}dtdrd\theta
d\varphi=2\pi\sqrt{-g}\Delta_{t}\Delta_{r}\Delta_{\theta}$ holds, where
$\Delta_{r}$ and $\Delta_{\theta}$ denote the intervals in Boyer-Linquist
radial and meridional coordinates, both of which are non-uniformly gridded.

In the initial state, there are no macro particles in any cell. 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$ for uniform grid intervals in
$x=r_{\ast}$ and $y=1-\cos\theta$ coordinates. In total, there are about
$3\times 10^{8}$ macro particles in the entire simulation region.

###  3.5 Nonstationary magnetosphere

It is checked a posteriori (§ 4.2) that the invariant grid intervals resolve
the skin depth

$l_{\rm p}=\frac{c}{\omega_{\rm p}},$ (52)

at every point at any elapsed time, where the plasma frequency $\omega_{\rm
p}$ is computed by the plasma density $n_{\pm}$ and its mean Lorentz factor
$\langle\gamma\rangle$ as

$\omega_{\rm p}=\sqrt{\frac{4\pi e^{2}n_{\pm}}{m_{\rm
e}\langle\gamma\rangle}}.$ (53)

For stellar-mass BHs, we obtain

$\frac{l_{\rm p}}{r_{\rm g}}=6.84\times
10^{-2}k\gamma_{6}{}^{1/2}\left(\frac{a}{M}\cdot\frac{B}{B_{\rm
eq}}\right)^{-1/2}M_{1}{}^{-1/4},$ (54)

where $\gamma_{6}\equiv\langle\gamma\rangle/10^{6}$, and

$k\equiv\kappa^{-1/2}\dot{m}^{-1/4}\left(\frac{r}{r_{\rm
H}}\right)^{5/8}\left(\frac{m_{\rm e}}{m_{\rm p}/1836}\right)^{1/2};$ (55)

$m_{\rm p}$ denotes the rest mass of the proton.

To resolve the skin depth, $l_{\rm p}$, most PIC simulations are performed in
one or two dimensions with small ion-to-electron mass ratios (e.g., Bohdan et
al. (2019)), covering limited time and length scales. In the present paper,
adopting a heavy electron mass of $m_{\rm e}=m_{\rm p}/20$, we perform a two-
dimensional PIC simulation within about ten gravitational radii over the time
duration of several hundred dynamical time scales. We rescale every expression
according to this modified electron mass. Under this assumption of heavy
charged leptons, we need about $10^{3}$ grid points to fully resolve the skin
depth.

On these grounds, we adopt a radial grid of 960 uniform cells between
$-15.81807<x<12.82443$, which corresponds to $1.46514M<r<13.68548M$, and 1920
uniform cells between $0<y=1-\cos\theta<2$, which corresponds to
$0^{\circ}<\theta<180^{\circ}$. We present the explicit expression of the
discretized Maxwell equations and the particle equation of motion in Appendix
E. To construct the electric currents from particle motion (Appendix F), we
adopt the area weighting (Villasenor & Buneman, 1992).

####  3.5.1 Electromagnetic fields

We start with the magnetic-field-aligned electric field. In figure 3, we
present the distribution of
$\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}/[B_{\rm eq}(2M)]^{2}$. For
$M=10M_{\odot}$ and $\dot{m}=2.25\times 10^{-4}$, we obtain $B_{\rm
eq}(2M)=1.452\times 10^{6}$ G. The left and right panels show the distribution
of $\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}/[B_{\rm eq}(2M)]^{2}$ at $t=0$
and $t=430M$, respectively. In the northern hemisphere, since
$F_{\theta\varphi}>0$ holds, a positive (or a negative) sign of
$\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}$ indicates that an outward (or an
inward) electric field arises along the local magnetic field lines. In the
same manner, in the southern hemisphere, it follows from $F_{\theta\varphi}<0$
that a positive (or negative) sign of
$\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}$ means an inward (or an outward)
electric field.

Figure 3:  Distribution of $\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}/B_{\rm
eq}(2M){}^{2}$ on the poloidal plane ($r$,$\theta$) for $a=0.9M$, $B=B_{\rm
eq}$, and $\dot{m}=2.25\times 10^{-4}$. Left and right panels show the
distribution at $t=0.00$ and $t=430.00GMc^{-3}$, respectively. The abscissa
and ordinate are common with figure 2. Magnetic surfaces (i.e., constant
$A_{\varphi}$ lines) are not depicted, but are common with the right panel of
figure 1 Figure 4:  Radial component of the Blandford-Znajek (BZ) flux as a
function of the elapsed time at radius $r=r_{\rm H}+0.25M$. The ordinate is
normalized by its analytical estimate (see text for details), while the
abscissa by the dynamical time scale. Each panel shows the BZ flux at discrete
colatitudes as labeled.

At $t=0$, there are no poloidal currents in the magnetosphere because of the
lack of plasmas. Thus, the magnetic field has no toroidal component.
Accordingly, the rotational energy of the BH is not being extracted. On the
other hand, because of the relative motion of the magnetic field lines with
respect to the space time, there appears a strong electromagnetic field along
the mangetic field line. The left panel of figure 3 shows that such a
magnetic-field-aligned electric field is exerted in both hemispheres.

As time elapses, $|\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}|$ evolves,
exerting electric currents in the magnetosphere. These currents modify the
poloidal electric field through the Ampere’s law. For example, at $t=430M$ in
the northern hemisphere, the right panel shows that a negative
$\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}$ appears in the higher-middle
latitudes, exerting inwards currents there. In the same manner, in the
southern hemisphere, $E$ points inwards in the higher-middle latitudes outside
the ergosphere. Note that $F_{\theta\varphi}$ changes sign across the equator.

Because of this nearly symmetric distribution of the electric field between
the two hemispheres, electrons (or positrons) are accelerated outwards (or
inwards) in the higher-middle latitudes in both hemispheres. In the lower
latitudes, however, positrons (or electrons) are accelerated outwards (or
inwards) by the strong magnetic-field-aligned electric field in both
hemispheres. As a result, currents flow inwards in the middle latitudes and
flow outwards in the lower latitudes. For the details of the electric
currents, see § 3.5.2. Such a pattern of the magnetospheric currents lead to a
positive toroidal magnetic field ($\propto F^{r\theta}$) in the northern
hemisphere, while a negative $F^{r\theta}$ in the southern hemisphere. Thus,
$F^{r\theta}$ vanishes on the equator.

Now let us consider the Poynting flux, or equivalently the BZ flux,

$T_{\rm em}{}^{r}{}_{t}=\frac{c}{4\pi}F^{r\mu}F_{\mu t}\propto
F^{r\theta}F_{\theta t}$ (56)

in the BH vicinity. If it becomes positive, it means that the BH’s rotational
energy is being extracted electromagnetically. In figure 4, we plot the
temporal evolution of $T_{\rm em}{}^{r}{}_{t}(t,r,\theta)$ at $r=r_{\rm
H}+0.25M$ at four discrete colatitudes as labeled. The ordinate is normalized
by the typical BZ flux, which is analytically given by $F_{\rm
analytical}\equiv L_{\rm BZ}/S_{\rm area}$, where $S_{\rm
area}\equiv\int\int\sqrt{A}\sin\theta d\theta d\phi$. The BZ power (i.e., the
spin-down luminosity) can be estimated to be

$L_{\rm BZ}=\frac{1}{128}\left(\frac{a}{M}\right)^{2}B_{\perp}{}^{2}r_{\rm
H}{}^{2}c$ (57)

in the slow-rotating limit, $|a|\ll M$, where $B_{\perp}$ denotes the averaged
strength of the radial magnetic field. It follows from the figure that the
solution exhibits rapid plasma oscillations, as reported by Levinson & Cerutti
(2018). It also follows that the simulated BZ flux is consistent with its
analytical estimate. See the supplementary material of Crinquand et al. (2020)
for a consistent discussion between the simulated BZ flux and their analytical
estimates.

It should be noted that the BZ flux increases during the elapsed time
$320M<t<450M$ along the middle latitudes, $45^{\circ}\leq\theta\leq
75^{\circ}$. During this flux-enhancement phase, the BZ flux (i.e., the
Poynting flux) fluctuates relatively mildly compared to its amplitude within
$45^{\circ}\leq\theta\leq 60^{\circ}$. In the present magnetically dominated
magnetosphere, in which the magnetic energy density dominates the particles’
rest-mass energy densities, particles’ energy flux is typically less than
$10^{-7}$ compared to the Poynting flux; thus, we neglect their contribution
when we consider the energy flux.

In figure 5, we present the angular dependence of the BZ flux at four elapsed
time as labeled. During the flux-enhancement phase, the BH’s rotational energy
is efficiently extracted from the middle latitudes,
$40^{\circ}<\theta<75^{\circ}$ in the northern hemisphere, and
$105^{\circ}>\theta>135^{\circ}$ in the southern hemisphere.

Figure 5:  Radial component of the Blandford-Znajek flux as a function of the
colatitude angle, $\theta$. The blue dotted, red solid, black dashed, and
green dash-dotted curves show the BZ flux at time $t=330.00GMc^{-3}$,
$380.00GMc^{-3}$, $430.00GMc^{-3}$, and $480.00GMc^{-3}$, respectively.

To smear out the variation, we take a moving average with a period of
$5GMc^{-3}$, and plot the BZ fluxes in figure 6. In this particular figure, we
compare the results for three accretion rates: the top, middle, and bottom
panels show the BZ fluxes at $\dot{m}=0.000250$, $0.000225$, and $0.000200$,
respectively. We find that the flux is enhanced for typically $140\sim 180$
dynamical time scales, and that the flux peaks in the middle latitudes during
the enhancement irrespective of the accretion rate, as the blue dashed
($\theta=50^{\circ}$), the blue solid ($\theta=60^{\circ}$), and the black
dashed ($\theta=70^{\circ}$) curves indicate in each panel. For example, at
$\theta=60^{\circ}$, the FWHMs of the moving-averaged flux become $157M$ and
$149M$ in the northern and southern hemispheres, respectively, when
$\dot{m}=2.50\times 10^{-4}$ (as the top two panels show). They become $136M$
(northern) and $184M$ (southern) when $\dot{m}=2.25\times 10^{-4}$ (middle two
panels), and $141M$ (northern) when $\dot{m}=2.00\times 10^{-4}$ (bottom left
panel). However, it is still not possible to measure the FWHM for the
sourthern hemisphere when $\dot{m}=2.00\times 10^{-4}$ (bottom right panel).
The flux distributes more symmetrically between northern and southern
hemispheres, as the accretion rate increases.

Figure 6:  Moving-averaged BZ fluxes with a period of $5GMc^{-3}=5M$, for
three dimensionless accretion rates, $\dot{m}=2.50\times 10^{-4}$ (top),
$\dot{m}=2.25\times 10^{-4}$ (middle), and $\dot{m}=2.00\times 10^{-4}$
(bottom). Left (or right) panels depict the BZ flux in the northern (or
southern) hemispheres. Each curve denotes the BZ flux at colatitude $\theta$
as labeled.

####  3.5.2 Particle distribution and currents

We next consider the distribution functions of $e^{\pm}$’s, and the resultant
current distribution. Figure 7 shows the densities of electrons (left) and
positrons (middle) at the burst peak, $t=430M$, in log scale. The right panel
shows the GJ value at each point. It follows that both electrons and positrons
have greater densities than the GJ value particularly in the lower latitudes.
In the polar regions, because of the polarization drift caused by the varying
meridional electric field (and the constant radial magnetic field), electrons
migrate meridionally to accumulate in $\theta<15^{\circ}$ and
$\theta>165^{\circ}$, whereas positrons in $15^{\circ}<\theta<20^{\circ}$ and
$165^{\circ}>\theta>160^{\circ}$ at $t\sim 430M$. In the lower latitudes, the
leptonic densities attain $10^{7.7}\mbox{cm}^{-3}$.

The charged leptons carry electric currents, as depicted in figure 8. In the
left and right panels, we present the radial and meridional components of the
electric currents at each point in the ZAMO frame. For a quantity
$f(r,\theta)$, we plot

$F=\mbox{sign}\left\\{\lg\left[\max(|f|,1)\right],f\right\\},$ (58)

where $\mbox{sign}(a,b)=|a|$ if $b\geq 0$ and $=-|a|$ if $b<0$; we set
$f=J^{\hat{r}}$ and $J^{\hat{\theta}}$ for the left and right panels,
respectively. In the left panel, the yellow-red regions show that currents
flow inwards in the middle latitudes, while the blue-violet regions show that
they flow outwards in the lower latitudes. These radial currents are closed by
meridional currents flowing within the ergosphere, as the right panel shows.
For example, in the right panel, the blue (or red) region in the lower-middle
latitudes within $r<2M$ in the northern (or southern) hemisphere shows
equator-ward meridional currents. Because of this current closure, it is
confirmed that the BZ process is, indeed, facilitated.

To grasp the current distribution more easily, we plot the direction and
strength of the poloidal currents as red arrows in figure 9. The current
density is averaged over the area in which we compute the direction and length
of each arrow. In the figure, the length of the arrows indicates the strength
of the current density in logarithmic scale, as indicated by the right panel.
We can confirm the current pattern discussed in the foregoing paragraph. It
also follows that the averaged currents mostly flow in the middle and lower
latitudes; thus, the low density regions in the polar funnels do not
essentially affect the entire structure of the magnetosphere. In this figure,
we also plot the charge density, $(n_{+}-n_{-})/n_{\rm GJ}$, in color, where
$n_{+}$ and $n_{-}$ refer to the positronic and electronic number densities.
Values are plotted using the same method as figure 8 (eq. [58]). It follows
that the real charge density becomes even greater than the GJ value (right
panel of fig. 7), which indicates that the electron-positron plasmas become
highly non-neutral near the BH.

We next consider the distribution functions of the charged leptons. The
dimensionless distribution functions of $e^{\pm}$’s, $n_{\pm}$, are sliced
between the colatitude $\theta_{1}<\theta<\theta_{2}$,

$N_{\pm}(r_{\ast},\gamma)\equiv\frac{1}{n_{\rm
GJ}}\int_{\theta_{1}}^{\theta_{2}}d\theta\frac{\partial^{2}n_{\pm}(r_{\ast},\theta,\gamma)}{\partial\theta\partial\gamma}.$
(59)

In figure 10, we present $N_{-}$ and $N_{+}$ as a function of $r_{\ast}$ and
$\gamma$ (i.e., the Lorentz factor) in the left and right columns,
respectively. The range of $\theta\in[\theta_{1},\theta_{2})$ increases from
upper to lower rows as described in the caption. It shows that the Lorentz
factors are saturated at a terminal value at each point. The terminal value is
determined by the balance between the electrostatic acceleration and the
synchro-curvature radiation drag force. Since
$|\mbox{\boldmath$E$}\cdot\mbox{\boldmath$B$}|$ increases with decreasing
radius significantly at $r_{\ast}<0$ (or equivalently, at $r<4.1M$), particles
gain large kinetic energies inside $r<4M$. Nevertheless, by virtue of the
residual magnetic-field-aligned electric field, particle motion is kept
relativistic in the entire region. Lorentz factors attain $\gamma>10^{4.5}$ in
$r_{\ast}<1M$ (or equivalently, in $r<4.7M$).

To examine the pitch-angle dependence, in figure 11, we plot
$N_{\pm}(r_{\ast},\chi)$, where the pitch angle $\chi$ becomes $+1$ (or $-1$)
when a particle is moving outwards (or inwards) with a very small trans-
magnetic-field momentum, and becomes $0$ when it has no longitudinal momentum.
At $r_{\ast}>0$, particles have small trans-magnetic-field momenta. However,
in the higher-middle latitudes, $\theta\sim 40^{\circ}$, the middle two panels
show that electrons (or positrons) migrate outwards with positive $\cos\chi$
(or inwards with negative $\cos\chi$), indicating an inward poloidal current
there. Figures 10 and 11 show that particles are highly relativistic with non-
thermal and anisotropic distribution functions.

Figure 7:  Densities of electrons (left) and positrons (middle), and the
Goldreich-Julian density (right), $n_{\rm GJ}$, in ${\rm cm}^{-3}$ unit at
$t=430M$. Values are plotted in decadic logarithm. The abscissa and ordinate
are common with figure 2.  Figure 8:  Electric currents measured by ZAMO in
the poloidal plane at $t=430M$. The abscissa and ordinate are common with
figure 2. The left panel shows the radial component, $J^{\hat{r}}$, while the
right one does the meridional component, $J^{\hat{\theta}}$. To plot the
values, we take the decadic lorarithm of the absolute value of a quantity,
then put the same sign as the quantity. For example, in the left panel, the
value of $10^{6}$ (or $-10^{6}$) corresponds to an outward (or an inward)
current density whose absolute value is $10^{6}\mbox{ statampere cm}^{-2}$. In
the same way, in the the right panel, a positive (or a negative)
$J^{\hat{\theta}}$ means an equator-ward current in the northern (or southern)
hemisphere.  Figure 9:  Dimensionless charge denisity $(n_{+}-n_{-})/n_{\rm
GJ}$ (color image) and poloidal electric current (red arrows). The abscissa
and ordinate are common with figure 2. The charge density is plotted in linear
scale. The green (or yellow) regions show positive (or negative) dimensionless
charge densities. The right panel shows four example arrow lengths
corresponding to the indicated strengths of the poloidal current densities in
$\mbox{statampere cm}^{-2}$.  Figure 10:  Lorentz factor dependence of the
distribution of electrons (left) and positrons (right) as a function of the
dimensionless tortoise coordinate, $r_{\ast}/M$, at elapse time $t=430M$. The
color images show the logarithmic value of the dimensionless distribution
functions (see text) along discrete colatitudes. The top, middle, and bottom
rows show the distribution function between the meridional range
$16.606^{\circ}-23.568^{\circ}$, $37.678^{\circ}-41.432^{\circ}$, and
$57.235^{\circ}-60.034^{\circ}$, respectively.  Figure 11:  Similar figure as
figure 10, but the pitch-angle dependence of the distribution functions of
electrons (left) and positrons (right) at $t=430M$. The color images show the
logarithmic value of the dimensionless distribution functions within the same
meridional range as fig. 10.

##  4 Discussion

To sum up, we simulated the evolution of a BH magnetosphere by a PIC scheme,
when a poloidal magnetic field is sustained by a disk toroidal current, and
when the electron-positron pair plasmas are steadily supplied homogeneously
per invariant volume basis. Provided that the mass accretion rate is much less
than the Eddington rate, both the electromagnetic fields and the particle
distribution functions exhibit rapid variability. The rotational energy of the
BH is, indeed, extracted via the Blandford-Znajek process, whose energy flux
concentrates in the middle latitudes, particularly during the flux-enhancement
phase that lasts approximately $160\pm 20$ dynamical time scales. We have
demonstrated that the collision time scale is much longer than the gyration
time scale (i.e., the Ohm’s law cannot be justified), that the pair plasma is
highly non-neutral, that the particles’ energy distribution is non-Maxwellian,
and that the momentum distribution is anisotropic. Thus, we must discard the
MHD approximation when we consider the jet launching region around the BH
whose mass accretion rate is highly sub-Eddington.

In this section, we discuss the dominant radiative process in § 4.1, the
validity of gridding in § 4.2, comparison with other works in § 4.3, and an
implication to the collimation of VLBI jets in § 4.4.

###  4.1 Dominant radiative process

Although the synchro-curvature radiation process is incorporated in the
radiative reaction force, $F_{\rm rad}^{j}$, ICS are not considered as a
radiation drag. Thus, we have to confirm that the ICS process is negligible
compared with the synchro-curvature process around stellar-mass BHs accreting
at $\dot{m}\ll 1$. Here we briefly compare the pure curvature, pure
synchrotron, and ICS processes, and discuss the dominant process. To make a
general discussion, we adopt the actual electron mass $m_{\rm e}=m_{\rm
p}/1836$, instead of $m_{\rm p}/20$, in this subsection.

The magnitude of the radiation drag force due to the pure curvature process is
given by

$\displaystyle F_{\rm curv}$ $\displaystyle=$
$\displaystyle\frac{2}{3}e^{2}\frac{\gamma^{4}}{\rho_{\rm c}^{2}}$ (60)
$\displaystyle=$ $\displaystyle 1.70\times 10^{-4}\left(\frac{\rho_{\rm
c}}{2M}\right)^{-2}M_{1}{}^{-2}\gamma_{7}{}^{4}\mbox{ dyn},\qquad$

where $\rho_{\rm c}$ refers to the curvature radius of the particle’s center
of gyration, $2M=2MGc^{-2}$ the Schwarzschild radius, and
$\gamma_{n}\equiv\gamma/10^{n}$. If this force balances with the electrostatic
force,

$eE_{\parallel}=4.80\times 10^{-5}|E_{\parallel,5}|\mbox{ dyn},$ (61)

we find that the particles saturate at the Lorentz factor,

$\gamma_{\rm curv}=7.28\times 10^{6}\left(\frac{\rho_{\rm
c}}{2M}\right)^{1/2}M_{1}{}^{1/2}E_{\parallel,5}{}^{1/4},$ (62)

where $E_{\parallel,5}\equiv E_{\parallel}/(10^{5}\mbox{ statvolt cm}^{-1})$.

The magnitude of the radiation-drag force due to the pure synchrotron process
is given by

$\displaystyle F_{\rm sync}$ $\displaystyle=$
$\displaystyle\frac{2}{3}r_{0}{}^{2}\gamma^{2}B^{2}\sin^{2}\chi$ (63)
$\displaystyle=$ $\displaystyle 5.29\times
10^{-5}\gamma_{5}{}^{2}B_{6}{}^{2}\frac{\sin^{2}\chi}{0.1}\mbox{ dyn},$

where $r_{0}$ denotes the classical electron radius, and $\chi$ the pitch
angle. Equating equations (61) and (63), we obtain the terminal Lorentz factor

$\gamma_{\rm sync}=9.52\times
10^{4}B_{6}{}^{-1}\left(\frac{\sin^{2}\chi}{0.1}\right)^{-1/2}E_{\parallel,5}^{1/2}.$
(64)

If $\gamma_{\rm sync}<\gamma_{\rm curv}$, the pure synchrotron process
dominates the pure curvature process. If $m_{\rm e}=m_{\rm p}/1836$, We find
that the condition $\gamma_{\rm sync}<\gamma_{\rm curv}$ is satisfied, because

$B_{6}\left(\frac{\sin\chi}{0.1}\right)E_{\parallel,5}^{-1/4}\left(\frac{\rho_{\rm
c}}{2M}\right)^{1/2}>1$ (65)

is satisfied by the vast majority of the particles. In another word, the
synchro-curvature process reduces to the pure synchrotron process for most of
the particles. However, since we assume heavy electrons in this paper, the
increased gyro radius makes the synchrotron process less efficient; thus, the
radiation drag force is given by the synchro-curvature process in general for
heavy electrons.

We next compare the pure synchrotron process with the ICS. At radius $r$ from
the central BH, the number density of the ADAF synchrotron photons are given
by (Mahadevan, 1997)

$N_{\rm ph}=\frac{L_{\rm sync}}{4\pi r^{2}h\nu c},$ (66)

where $h\nu$ refers to the energy of the photons emitted from the ADAF via the
synchrotron process. The synchrotron luminosity is given by

$L_{\rm sync}=1.67\times 10^{36}M_{1}{}^{1/2}\dot{m}^{3/2}\mbox{ ergs
s}^{-1}.$ (67)

The typical photon energy is given by equation (22) of Mahadevan (1997).
Accordingly, we obtain

$N_{\rm ph}=1.55\times 10^{20}T_{\rm e,9}{}^{5}M_{1}\dot{m}\mbox{ photons
cm}^{-3},$ (68)

where $T_{\rm e,9}$ refers to the electron temperature, $T_{\rm e}$,
normalized by $10^{9}K$. The ICS drag force per electron (or positron) is
given by

$F_{\rm ICS}\approx N_{\rm ph}\sigma_{\rm KN}\gamma m_{\rm e}c^{2},$ (69)

where

$\sigma_{\rm KN}\approx\frac{3}{8}\sigma_{\rm T}x^{-1}(\ln 2x+\frac{1}{2}),$
(70)

refers to the Klein-Nishina cross section, and $x\approx\gamma$ holds on
average. We thus have

$F_{\rm ICS}\approx 3.89\times 10^{-10}M_{1}\dot{m}T_{\rm e,9}{}^{5}$ (71)

We thus obtain

$\gamma_{5}{}^{2}B_{6}{}^{2}\left(\frac{\sin\chi}{0.1}\right)^{2}M_{1}^{-1}\dot{m}^{-1}T_{\rm
e,9}{}^{-5}\gg 10^{-5}$ (72)

On these grounds, we obtain $F_{\rm sync}\gg F_{\rm curv}\gg F_{\rm ICS}$ for
the actual electron mass. Note that the Lorentz factor will increase with
decreasing electron mass. However, for heavy electrons, energy transfer
efficiency increases due to their greater mass ratio to protons. Thus, ICS may
not be negligible even for stellar-mass BHs. However, we neglected ICS in this
paper, considering future extension to smaller electron masses.

In short, for the actual electron mass, it is possible that the ICS process is
negligible compared to the synchro-curvature process, when we consider
stellar-mass BHs in a quiescent state. However, for supermssive BHs, we
generally obtain $F_{\rm ICS}>F_{\rm curv}\gg F_{\rm sync}$ because of the
large curvature radius and the weak magnetic field strength (Hirotani et al.,
2016) for the actual electron mass.

###  4.2 Grid interval versus skin depth

Let us compare the invariant grid interval $\epsilon r_{\rm
g}\equiv\Delta_{s}\equiv{\rm
max}\left(\sqrt{g_{rr}}\Delta_{r},\sqrt{g_{\theta\theta}}\Delta_{\theta}\right)$
with the skin depth, $l_{\rm p}$ (eq. 52), where $\Delta_{r}$ and
$\Delta_{\theta}$ denote the interval in $r$ and $\theta$ coordinates,
respectively. Representative values of the dimensionless function
$\epsilon(r,\theta)$ are presented in table 1. Because of a constant
$\sin\theta\Delta_{\theta}$ gridding, $\Delta_{\theta}$ increases with
decreasing $\sin\theta$. Note that the dimensionless grid interval $\epsilon$
is defined independently from the BH mass.

Adopting heavy electrons, $m_{\rm e}=m_{\rm p}/20$, and normalizing the skin
depth with $r_{\rm g}$, we obtain

$\frac{l_{\rm p}}{r_{\rm g}}=3.4\left(\gamma_{5}n_{5}{}^{-1}\right)^{1/2}$
(73)

where $\gamma_{\nu}\equiv\langle\gamma\rangle/10^{\nu}$ and $n_{\nu}\equiv
n_{\pm}/10^{\nu}$; $\langle\gamma\rangle$ denotes the averaged Lorentz factor.
The highest density region appears in the lower latitudes at $r<4M$, i.e., at
$r_{\ast}<-0.25M$. In this region, $n<10^{7.7}\mbox{ cm}^{-3}$ and
$\langle\gamma\rangle>10^{4.7}$; thus, we obtain $l_{\rm p}>0.108r_{\rm
g}>5\Delta_{\rm s}$ there. Other regions have greater skin depth. Accordingly,
the skin depth can be resolved at every point (and in fact at every time step)
during the PIC simulation, although it is marginal when $l_{\rm p}\sim
5\Delta_{\rm s}$ happens at the flux peak.

Table 1
Invariant grid interval, $\epsilon(r,\theta)=\Delta_{s}/r_{\rm g}$ $r/r_{\rm
g}$ $\theta=2.62^{\circ}$ $29.0^{\circ}$ $60.0^{\circ}$ $90.0^{\circ}$ 8.507
.2762 .0262 .0261 .0260 4.099 .1355 .0217 .0215 .0214 2.021 .0714 .0122 .0117
.0114 1.471 .0557 .0036 .0028 .0026

###  4.3 Comparison with previous works

Let us compare the present work with recent two works on 2-D GR PIC
simulations (Parfrey et al., 2019; Crinquand et al., 2020).

First, we discuss the difference with Parfrey et al. (2019). Ignoring
radiative transfer, and considering instead an injection of pairs whose rate
is proportional to the local magnetic-field-aligned electric field, they
performed 2-D GR PIC simulations of the magnetosheres of an extremely rotating
BH with $a=0.999M$. The BH mass is not explicitly specified in their paper,
but presumably supermassive. They adopted an extremely small magnetic field
strength to emphasize the Penrose process, an energy-extraction mechanism from
rotating BHs when particles fall onto the horizon with negative energies
measured at infinity. Particles were created in the reconnecting current sheet
on the equatorial plane. The magnetization parameter was $2000$; that is, the
magnetosphere was still magnetically dominated. Magnetic field configuration
is initially the Wald’s vacuum solution, which is produced by a toroidal ring
current flowing on the equatorial plane at large distances. Then the field
lines are bend back toward the BH to penetrate the horizon. For M87*, their
dimensionless magnetic field strength, $\tilde{B}_{0}=10^{3}$ corresponds to
the actual strength $B=\tilde{B}_{0}(m_{\rm e}c^{2}/er_{\rm g})\sim 10^{-9}$
G, which is about $10^{-11}$ times smaller than what is observed (Event
Horizon Telescope Collaboration et al., 2019). Under this condition, they
found that the charged leptons created in the current sheet plunge onto the
horizon with negative energies as a result of the interaction with the
electromagnetic field, and that the Penrose process contribute in the
extraction of the energy and angular momentum from a maximally rotating BH.

In the present paper, on the other hand, we consider a stellar-mass BH. In
this case, we can resolve the skin depth with a realistic magnetic field
strength, $B=B_{\rm eq}$. Accordingly, the magnetization parameter becomes of
the order of $10^{7}$; thus, we did not mention about the particle
contribution on the energy extraction from a rotating BH. Nevertheless,
particle energy flux does become negative in our simulation as well; thus, the
Penrose process is indeed working, although their contribution is negligible.
We assume a homogeneous and constant plasma supply in the present work, which
forms a constast to Parfrey et al. (2019), who considered pair creation in the
reconnecting current sheet. The magnetic field is assumed to be fixed on the
poloidal plane during the simulation in the present paper.

Second, let us discuss the difference with Crinquand et al. (2020). In their
2D GR PIC code, implementing inverse-Compton scatterings (ICS) and photon-
photon pair production self-consistently, they solved the PIC equations
together with the radiative transfer equation. Electron-positron pair plasmas
are supplied in response to gap opening. They presumed supermassive BHs. In
this case, ICS dominates the pure-curvature radiation (Hirotani et al., 2016),
where the particle’s pitch angles become small enough due to the ICS and
continuous acceleration along the magnetic field lines (Beskin et al., 1992).
Thus, their treatment, which takes account of only the ICS as radiative
processes, can be justified. Crinquand et al. (2020) assumed a true monopole
magnetic field, whose radial component point outwards in the entire
magnetosphre. Accordingly, there appear no current sheets in their solution.
Electric currents flow inwards (or outwards) in the northern (or southern)
hemisphere, leading to a negative $F^{r\theta}$ in the entire magnetosphere.
Since $E_{\theta}=F_{\theta t}>0$ holds in both hemispheres, the resultant
Poynting flux is also positive in both hemispheres. In their Zelton code, they
solved all the six components of the electromagetic fields, computing the
evolution of $F_{\varphi t}$ from the toroidal component of the current
density $J^{\varphi}$ that is constructed from the azimuthal motion of the
charged particles.

In the present paper, on the other hand, instead of solving the radiative
transfer equation, we simply adopt a uniform and constant pair production in
the particle source term. We consider a stellar-mass BH, in which case
synchro-curvature process becomes non-negligible compared to ICS (Hirotani et
al., 2016). We assume a split-monopole magnetic field, whose radial component
is positive (or negative) in the northern (or southern) hemisphere.
Accordingly, there appears a current sheet on the equator. Electric currents
flow inwards in the higher latitudes in both hemispheres, and flow outwards in
the lower latitudes, leading to a negative (or positive) $F^{r\theta}$ in the
northern (or southern) hemisphere, and a positive BZ flux in both hemispheres,
as discussed in §3.5.1.

Let us briefly point out the limitation of the present analysis. Unlike the
Zelton code, we solved only three of the six components of the electromagnetic
fields in the present paper, putting $F_{\varphi t}=0$, and hence
$\partial_{t}F_{\theta\varphi}=\partial_{t}F_{\varphi r}=0$ throughout the
simulation, as the first step. If the local physics (as supposed in recent 1-D
GR PIC simulations) does not significantly affect the entire structure of the
magnetosphere, such an assumption of a fixed poloidal magnetic field could be
justified. However, when we apply the present 2D method to a global
magnetosphere, toroidal currents carried by the charged leptons at each
position can alter the entire magnetic field configuration on the poloidal
plane. Thus, in the same manner as in recent 1-D models Levinson & Cerutti
(2018); Chen & Yuan (2020); Kisaka et al. (2020), we must draw attention to
the limitation of $F_{\varphi t}=0$. In the next paper, we will extend our
analysis to the case $F_{\varphi t}\neq 0$, constructing $J^{\varphi}$ at each
position from the actual motion of the charged leptons, and solving all the
six components of the electromagnetic fields.

###  4.4 Implication for supermassive black holes

Let us finally discuss what can be expected if the present results obtained
for stellar-mass BHs can be applied to supermassive BHs. We have demonstrated
that the BH’s rotational energy is efficiently extracted along the magnetic
field lines that cross the event horizon in the middle latitudes. Let us
discuss an implication of this result on the formation of limb-brightened
jets. It has been revealed that the innermost region of the M87 jet exhibits a
limb-brightened structure by the VLBA observations in 15–43 GHz (Kovalev et
al., 2007; Ly et al., 2007; Junor et al., 1999; Hada et al., 2011, 2013;
Walker et al., 2018). At 86 GHz, this limb-brightened structure is already
well developed at 0.15 mas from the VLBI core, where the corresponding
apparent opening angle becomes approximately $100^{\circ}$ (Hada et al.,
2016). If a relatively large viewing angle of $\theta_{\rm view}\sim
30^{\circ}$ is adopted (Ly et al., 2007; Hada et al., 2016), the jet has a
deprojected opening angle $\chi_{\rm open}\sim 50^{\circ}$ at the deprojected
distance $z=84GMc^{-2}$. However, if a smaller view angle, $\theta_{\rm
view}\sim 17^{\circ}$, is adopted (Biretta et al., 1999; Wang & Zhou, 2009;
Perlman et al., 2011; Nakamura & Meier, 2014; Mertens et al., 2016; Walker et
al., 2018), we obtain $\chi_{\rm open}\sim 30^{\circ}$ at $z=108GMc^{-2}$.

If a jet begins to collimate outside the outer light surface (Camenzind,
1986), we may assume that the jet is radial within the distance
$\zeta\varpi_{\rm LC}$ from the rotation axis, and becomes paraboloidal
outside of it (Asada & Nakamura, 2012; Hada et al., 2013; Nakamura & Asada,
2013; Asada et al., 2016; Nakamura et al., 2018), where $\varpi_{\rm LC}\equiv
c/\Omega_{\rm F}$ denotes the typical radius of the outer light surface.
Figure 12 sketches the geometry of this jet downstream region. Assuming that
we observe the jet at the position ($r$,$\theta$), we can express $\zeta$ in
terms of the observables $z$ and $\theta=\chi_{\rm open}/2$ as
$r(1-\cos\theta)=(\zeta\varpi_{\rm LC}/\sin\theta_{0})(1-\cos\theta_{0})$,
which gives

$\zeta\approx\frac{z}{\varpi_{\rm
LC}}\frac{1-\cos\theta}{\cos\theta}\frac{\sin\theta_{0}}{1-\cos\theta_{0}},$
(74)

where $z=r\cos\theta$. If $\Omega_{\rm F}\approx 0.5\omega_{\rm H}$, $a\approx
0.9M$ gives $\varpi_{\rm LC}\approx 6.4GM/c^{2}$. (Note that $\Omega_{\rm
F}=F_{t\theta}/F_{\theta\varphi}$ is not assumed but solved in the PIC
simulation.) If the magnetic field line with the footpoint angle
$\theta_{0}\approx 60^{\circ}$ (or $75^{\circ}$) is brightened at downstream
$z$ with half opening angle $\theta=\chi_{\rm open}/2$, we obtain
$\zeta\approx 2.4$ (or $\zeta\approx 1.8$) for $\theta_{\rm view}=30^{\circ}$,
and $\zeta\approx 1.0$ (or $\zeta\approx 0.78$) for $\theta_{\rm
view}=17^{\circ}$, using the 86 GHz VLBI observations. On these grounds, if
the BZ flux concentrates in the middle latitudes also around supermassive BHs,
it is possible that the M87 jet begins to collimate slightly outside the outer
light surface, typically within the distance $2.4\varpi_{\rm LC}$ from the
rotation axis, which may become a good target of the Event Horizon Telescope
and GRAVITY.

Figure 12:  Schematic picture of the jet downstream region. The jet shape is
assumed to be radial inside the radius $r_{0}$, and paraboidal outside of it.
The jet is observed with VLBI at radius $r$. Both $\theta$ and $\theta_{0}$
are measured from the jet axis, which coincides with the BH’s rotation axis.

We thank the anonymous referee for careful reading and valuable comments to
this manuscript. The authors acknowledge grant support to the CompAS/Antares
group, and the access to high-performance computing facility from the
Institute of Astronomy and Astrophysics in Academia Sinica (ASIAA). The
authors acknowledge grant support from Ministry of Science and Technology
(MoST) under 105-2119-M-001-044-MY3, 108-2112-M-001-009- and
109-2112-M-001-028-.

##  Appendix A ZAMO-measured quantities

We give the expressions of ZAMO-measured quantities, which are indicated by
tilde ($\,\tilde{}\,$), in this appendix section. Using equations (4)–(7), we
obtain

$u_{\hat{r}}=\sqrt{\frac{\Delta}{\Sigma}}u_{r},$ (A1)
$u_{\hat{\theta}}=\sqrt{\frac{1}{\Sigma}}u_{\theta},$ (A2)
$u_{\hat{\varphi}}=\sqrt{\frac{A}{\Sigma}}\sin\theta u_{\varphi},$ (A3)

Or equivalently, their contravariant components become

$u_{\hat{r}}=u^{\hat{r}}\equiv\sqrt{g_{rr}}u^{r},$ (A4)
$u_{\hat{\theta}}=u^{\hat{\theta}}\equiv\sqrt{g_{\theta\theta}}u^{\theta},$
(A5)
$u_{\hat{\varphi}}=u^{\hat{\varphi}}\equiv\sqrt{g_{\varphi\varphi}}(u^{\varphi}-\omega
u^{t}),$ (A6) $u^{j}\equiv\gamma v^{j},$ (A7)

where

$\gamma\equiv\frac{1}{\sqrt{1-\vec{v}\cdot\vec{v}}}.$ (A8)

In the same way, the tetrad transformation law gives the following expressions
of the elctromagnetic fields in ZAMO:

$\displaystyle\tilde{E}_{\hat{r}}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$F$}(\tilde{\mbox{\boldmath$e$}}_{(\hat{r})},\tilde{\mbox{\boldmath$e$}}_{(\hat{t})})$
$\displaystyle=$ $\displaystyle\frac{\sqrt{(r^{2}+a^{2})^{2}-\Delta
a^{2}\sin^{2}\theta}}{\Sigma}\left(F_{rt}+\omega F_{r\varphi}\right)$
$\displaystyle\tilde{E}_{\hat{\theta}}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$F$}(\tilde{\mbox{\boldmath$e$}}_{(\hat{\theta})},\tilde{\mbox{\boldmath$e$}}_{(\hat{t})})$
$\displaystyle=$ $\displaystyle\frac{\sqrt{(r^{2}+a^{2})^{2}-\Delta
a^{2}\sin^{2}\theta}}{\sqrt{\Delta}\Sigma}\left(F_{\theta t}+\omega
F_{\theta\varphi}\right)$
$\tilde{E}_{\hat{\varphi}}=\mbox{\boldmath$F$}(\tilde{\mbox{\boldmath$e$}}_{(\hat{\varphi})},\tilde{\mbox{\boldmath$e$}}_{(\hat{t})})=\frac{1}{\rho_{\rm
w}}F_{\varphi t}=\frac{1}{\sqrt{\Delta}\sin\theta}F_{\varphi t}$ (A11)
$\tilde{B}_{\hat{r}}={}^{\ast}\\!\mbox{\boldmath$F$}(\tilde{\mbox{\boldmath$e$}}_{(\hat{t})},\tilde{\mbox{\boldmath$e$}}_{(\hat{r})})=F_{\hat{\theta}\hat{\varphi}}=\frac{1}{\sqrt{\Sigma
g_{\varphi\varphi}}}F_{\theta\varphi},$ (A12)
$\tilde{B}_{\hat{\theta}}={}^{\ast}\\!\mbox{\boldmath$F$}(\tilde{\mbox{\boldmath$e$}}_{(\hat{t})},\tilde{\mbox{\boldmath$e$}}_{(\hat{\theta})})=-F_{\hat{r}\hat{\varphi}}=\frac{\sqrt{\Delta}}{\sqrt{\Sigma
g_{\varphi\varphi}}}F_{\varphi r},$ (A13)
$\tilde{B}_{\hat{\varphi}}={}^{\ast}\\!\mbox{\boldmath$F$}(\tilde{\mbox{\boldmath$e$}}_{(\hat{t})},\tilde{\mbox{\boldmath$e$}}_{(\hat{\varphi})})=\sqrt{g_{rr}g_{\theta\theta}}F^{r\theta}=\frac{\Sigma}{\sqrt{\Delta}}F^{r\theta},$
(A14)

where ${}^{\ast}\\!\mbox{\boldmath$F$}$ refers to the Maxwell tensor, which is
the dual of the Faraday tensor $F$. For example, figure 2 depicts the
distribution of equations (A12) and (A13) on the poloidal plane. We could use
equations (A4) and (A5) to compute $J^{\hat{r}}$ and $J^{\hat{\theta}}$, which
are depicted in figures 8–9. However, in actual calculations, we compute
$J^{r}$ and $J^{\theta}$ from equations (46) and (47) and convert them into
the ZAMO-measured quantities, using the transformation law of the one-form
bases that are dual to the ZAMO’s tetrad, (4)–(7).

##  Appendix B Stationary, vacuum magnetosphere

In this appendix section, we formulate the basic equations to solve the
initial, stationary electromagnetic fields.

###  B.1 The Gauss’s law

The inhomogeneous part of the Maxwell equations can be written as

$\nabla_{\mu}F^{\nu\mu}=\frac{4\pi}{c}J^{\nu},$ (B1)

where $\nabla_{\mu}$ denotes the covariant derivative with respect to the
coordinate variable $x^{\mu}$, $F^{\nu\mu}$ the electro-magnetic field
strength tensor, and $J^{\nu}$ the four current density. Putting $\nu=0$, we
obtain the Gauss’s law,

$\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}F^{0\mu})=4\pi\rho,$ (B2)

where $\rho$ refers to the electric charge density, and
$\sqrt{-g}=\Sigma\sin\theta$. In the Kerr spacetime, we obtain

$\displaystyle F^{01}$ $\displaystyle=$ $\displaystyle
F^{tr}=\frac{g^{rr}}{\rho_{\rm
w}{}^{2}}\left(-g_{\varphi\varphi}F_{tr}+g_{t\varphi}F_{\varphi r}\right),$
(B3) $\displaystyle F^{02}$ $\displaystyle=$ $\displaystyle
F^{t\theta}=\frac{g^{\theta\theta}}{\rho_{\rm
w}{}^{2}}\left(-g_{\varphi\varphi}F_{t\theta}+g_{t\varphi}F_{\varphi\theta}\right),$
(B4) $\displaystyle F^{03}$ $\displaystyle=$ $\displaystyle
F^{t\varphi}=\frac{1}{\rho_{\rm w}{}^{2}}F_{\varphi t}$ (B5)

At $t=0$, we assume stationary and axisymmetric magnetospehre,
$\partial_{t}=\partial_{\varphi}=0$, to obtain

$\displaystyle F_{rt}$ $\displaystyle=$ $\displaystyle
A_{t,r}-A_{r,t}=A_{t,r}=\partial_{r}A_{t},$ (B6) $\displaystyle F_{\theta t}$
$\displaystyle=$ $\displaystyle
A_{t,\theta}-A_{\theta,t}=A_{t,\theta}=\partial_{\theta}A_{t},$ (B7)
$\displaystyle F_{\varphi t}$ $\displaystyle=$ $\displaystyle
A_{t,\varphi}-A_{\varphi,t}=0,$ (B8) $\displaystyle F_{\varphi r}$
$\displaystyle=$ $\displaystyle
A_{r,\varphi}-A_{\varphi,r}=-\partial_{r}A_{\varphi},$ (B9) $\displaystyle
F_{\theta\varphi}$ $\displaystyle=$ $\displaystyle
A_{\varphi,\theta}-A_{\theta,\varphi}=\partial_{\theta}A_{\varphi},$ (B10)

where $A_{\mu}$= ($A_{t}$,$A_{r}$,$A_{\theta}$,$A_{\varphi}$) denotes the
vector potential; $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$.
Note that $F_{r\theta,\varphi}+F_{\theta\varphi,r}+F_{\varphi r,\theta}=0$
(i.e., $\nabla\cdot\mbox{\boldmath$B$}=0$) is automatically satisfied.

For an observer whose four velocity is $u^{\nu}$, the electromagnetic field
components are given by

$E_{\mu}=F_{\mu\nu}u^{\nu},$ (B11)

and

$B^{\mu}=\frac{1}{2}\eta_{\nu}{}^{\mu\rho\sigma}F_{\rho\sigma}u^{\nu},$ (B12)

where the completely anti-symmetric, Levi-Civita tensor density is defined by

$\eta_{\nu\mu\rho\sigma}=-\eta_{\mu\nu\rho\sigma}=-\eta_{\rho\mu\nu\sigma}=-\eta_{\sigma\mu\rho\nu}=-\eta_{\nu\rho\mu\sigma}=-\eta_{\nu\sigma\rho\mu}=-\eta_{\nu\mu\sigma\rho}$
(B13)

and

$\eta_{tr\theta\varphi}=-\frac{1}{\eta^{tr\theta\varphi}}=\sqrt{-g}.$ (B14)

Thus, for a distant static observer whose four velocity is
$\xi^{\mu}=(1,0,0,0)$, that is, the time-like Killing vector, the electric
field components are obtained by

$E_{r}=F_{rt},\quad E_{\theta}=F_{\theta t},\quad E_{\varphi}=F_{\varphi t}$
(B15)

in the Boyer-Lindquist coordinate. The magnetic field becomes

$\displaystyle B^{r}$ $\displaystyle=$
$\displaystyle\frac{-g_{tt}F_{\theta\varphi}+g_{t\varphi}F_{\theta
t}}{\sqrt{-g}},\quad$ $\displaystyle B^{\theta}$ $\displaystyle=$
$\displaystyle\frac{-g_{tt}F_{\varphi r}+g_{t\varphi}F_{rt}}{\sqrt{-g}},\quad$
$\displaystyle B_{\varphi}$ $\displaystyle=$
$\displaystyle-\sqrt{-g}F^{r\theta}.\quad$ (B16)

Substituing equations (B6), (B7), (B9) and (B10) into equation (B2), we obtain
the Gauss’s law:

$\displaystyle\frac{1}{\Sigma}\frac{\partial}{\partial
r}\left(\frac{A}{\Sigma}\frac{\partial A_{t}}{\partial r}\right)$
$\displaystyle+$
$\displaystyle\frac{1}{\Delta\Sigma\sin\theta}\frac{\partial}{\partial\theta}\left(\frac{A}{\Sigma}\sin\theta\frac{\partial
A_{t}}{\partial\theta}\right)+\frac{1}{\Sigma}\frac{\partial}{\partial
r}\left(\frac{2Mar}{\Sigma}\frac{\partial A_{\varphi}}{\partial r}\right)$
(B17) $\displaystyle+$
$\displaystyle\frac{1}{\Delta\Sigma\sin\theta}\frac{\partial}{\partial\theta}\left(\frac{2Mar}{\Sigma}\sin\theta\frac{\partial
A_{\varphi}}{\partial\theta}\right)=4\pi\rho.$

It follows that this Poisson equation contains only $A_{t}$ (i.e., the scalar
potential) if the BH is non-rotating (i.e., $a=0$). However, around a rotating
BH (i.e., if $a\neq 0$), $A_{\varphi}$ comes into the Gauss’s law. We thus
need one more differential equation that contains both $A_{t}$ and
$A_{\varphi}$.

###  B.2 The Biot and Savart law

To obtain an independent constraint on $A_{t}$ and $A_{\varphi}$, we consider
the Biot and Savart law. Putting $\nu=\varphi$ in equation (B1), we obtain

$\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}F^{\varphi\mu})=\frac{4\pi}{c}J^{\varphi},$
(B18)

Thus, equations (B3)–(B10) give

$\displaystyle-\frac{1}{\sin\theta}\frac{\partial}{\partial
r}\left(\frac{\Delta-a^{2}\sin^{2}\theta}{\Sigma}\frac{\partial
A_{\varphi}}{\partial
r}\right)-\frac{1}{\Delta}\frac{\partial}{\partial\theta}\left(\frac{\Delta-a^{2}\sin^{2}\theta}{\Sigma\sin\theta}\frac{\partial
A_{\varphi}}{\partial\theta}\right)$ (B19) $\displaystyle+$
$\displaystyle\frac{\partial}{\partial
r}\left(\frac{2Mar\sin\theta}{\Sigma}\frac{\partial A_{t}}{\partial
r}\right)+\frac{1}{\Delta}\frac{\partial}{\partial\theta}\left(\frac{2Mar\sin\theta}{\Sigma}\frac{\partial
A_{t}}{\partial\theta}\right)=4\pi\Sigma\sin\theta J^{\varphi},$

where $c=1$ is used. It follows that this equation contains only $A_{\varphi}$
(i.e., the azimuthal component of the vector potential) if $a=0$. However, if
$a\neq 0$, we must simultaneously solve the two second-order partial
differential equations (B17) and (B19) for $A_{t}$ and $A_{\varphi}$.

It should be noted that the coefficients of the
$\partial_{r}{}^{2}A_{\varphi}$ and the $\partial_{\theta}{}^{2}A_{\varphi}$
terms in equation (B19) change sign at the static limit, where $g_{tt}$
vanishes. Since the coefficients of these two highest-order derivative terms
change sign with respect to those of the other two highest-order derivative
terms (i.e., the $\partial_{r}{}^{2}A_{t}$ and the
$\partial_{\theta}{}^{2}A_{t}$ terms), the solution diverges during iterations
from the inside of the static limit (i.e., within the so-called “ergosphere”)
if we impose boundary conditions in the same way as standard elliptic type
partial differential equations. This ill behaviour is incurred because a
static observer (with respect to the star) becomes unphysical within the
ergosphere.

There are many ways to overcome this ill behaviour. In the present paper, we
solve this issue by adopting the ZAMO as a physical observer. Using the ZAMO’s
tetrad, equations (4)– (7), we can replace $A_{t}$ with the ZAMO-measured
scalar potential, $A_{\hat{t}}$ (eq. [12]), and obtain the two elliptic-type
equations (14) and (19), which describe the initial electromagnetic fields at
$t=0$.

##  Appendix C Pair production rate

Using the self-similar analytical solution of Newtonian ADAF model (Mahadevan,
1997), we find that the photons are mitted at the rate

$q_{\rm ff}=1.8\times 10^{2}\theta_{\rm
e}\dot{m}^{2}M_{9}{}^{-2}\left(\frac{r}{2M}\right)^{-3}\mbox{ergs
s}^{-1}\mbox{cm}^{-3}$ (C1)

by Bremsstrahlung, where $\theta_{\rm e}\equiv kT_{\rm e}/m_{\rm e}c^{2}$
denotes the dimensionless electron temperature. The luminosity of this
Bremsstrahlung emission can be computed as

$L_{\rm ff}(r)=2\pi\int_{\theta_{1}}^{\theta_{2}}\int_{r_{\rm min}}^{r}q_{\rm
ff}(r)r^{2}\sin\theta drd\theta\approx 2\pi\times 1.8\times 10^{2}\theta_{\rm
e}\dot{m}^{2}M_{9}{}^{-2}(2M)^{3}\ln\left(\frac{r}{2M}\right),$ (C2)

where $\theta_{1}\sim 60^{\circ}$ and $\theta_{2}\sim 120^{\circ}$ denote the
upper and lower boundary colatitudes of the ADAF.

In general, if one photon species (e.g., ADAF Bremsstrahlung photons or gap-
emitted gamma-rays) collide with another photon species (e.g., ADAF
Bremsstrahlung photons or ADAF synchrotron photons), the pair proction rate is
given by

$\dot{n}_{\pm}(r)=\int
d\nu_{\gamma}\alpha_{\gamma\gamma}\frac{1}{c}\int\frac{I_{\nu}}{h\nu}d\Omega_{\gamma},$
(C3)

where $d\Omega_{\gamma}$ refers to the photon propagation solid angle, and

$\alpha_{\gamma\gamma}=(1-\mu)\int_{\epsilon_{\rm
th}}^{\infty}\frac{dF_{\gamma}}{d\epsilon_{\gamma}}\sigma_{\gamma\gamma}d\epsilon_{\gamma}\approx(1-\mu)F_{\gamma}\sigma_{\gamma\gamma}$
(C4)

denotes the photon-photon absorption coefficient, and

$F_{\gamma}\equiv\int d\nu_{\gamma}\int\frac{I_{\nu}}{h\nu}d\Omega_{\gamma}$
(C5)

does the photon number flux; $I_{\nu}$, $\nu_{\gamma}$ and $\mu$ show the
specific intensity, photon frequency, and the cosine of the photon collision
angles, respectively.

In the present case, the two species are the same ADAF Bremsstrahlung photons.
Thus, we obtain

$\dot{n}_{\pm}(r)\approx(1-\mu)\sigma_{\gamma\gamma}\frac{F_{\gamma}{}^{2}}{c}$
(C6)

Evaluating the flux with the luminosity by

$F_{\gamma}(r)\approx\frac{L_{\rm ff}(r)}{2\pi r^{2}\epsilon_{\gamma}},$ (C7)

where the photon energy can be estimated to be $\epsilon_{\gamma}\approx
3\theta_{\rm e}m_{\rm e}c^{2}$, and assuming that the pair production rate is
uniform within $r<4M$, we obtain equation (50), where the relative velocity
factor is evaluated as $1-\mu=0.2$, and the pair-production cross section as
the $20\%$ of the Thomson cross section, $\sigma_{\gamma\gamma}=0.2\sigma_{\rm
T}$.

##  Appendix D Covariance of the radiation-reaction force

We explain here that the radiation reaction force contains not only those
obtained in a flat space (Cerutti et al., 2012, 2013),

$P_{\rm
rad,flat}=\frac{2}{3}e^{2}c\frac{q^{2}}{m^{2}}\gamma^{2}\left[-\left(\frac{\mbox{\boldmath$v$}}{c}\cdot\tilde{\mbox{\boldmath$E$}}\right)^{2}+\left(\tilde{\mbox{\boldmath$E$}}+\frac{\mbox{\boldmath$v$}}{c}\times\tilde{\mbox{\boldmath$B$}}\right)^{2}\right],\qquad$
(D1)

in locally homogeneous $\tilde{\mbox{\boldmath$E$}}$ and
$\tilde{\mbox{\boldmath$B$}}$ fields, but also contains the power of radiative
processes that result from any kinds of acceleration. To find this, we write
down equations (43) in a covariat form,

$\frac{du^{\mu}}{d\lambda}=-\Gamma^{\mu}{}_{\rho\sigma}u^{\rho}u^{\sigma}+\frac{q}{m}F^{\mu}{}_{\rho}u^{\rho},$
(D2)

where the friction term is dropped here in the right-hand side to draw the
conclusion from the first principles. Both the inertial forces (e.g., the
centrifugal and the Colioris forces) and the gravitational forces (due to the
spacetime curvature) are included in the connection coefficients,
$\Gamma^{\mu}{}_{\rho\sigma}$. For example, when a charge (or its guiding
center) moves along a curved magnetic field line in a flat space, we can
introduce a local cylindrical (or polar) coordinate system whose origin
resides at the center of the curved path of 3D particle motion. In this
coordinate, the $\Gamma^{\mu}{}_{\rho\sigma}u^{\rho}u^{\sigma}$ term gives the
centripetal acceleration, $a_{\perp}=c^{2}/R_{\rm c}$, where $R_{\rm c}$
refers to the curvature radius of the 3D particle motion. Since equation (45)
is covariant, the effect of such curvature radiation is included in any frame
of refence. Accordingly, we can take account of radiation effects due to any
acceleration. Note that a charged particle does radiate by gravitational
acceleration irrespective of the frame of reference.

In actual program, we can compute the radiation-reaction force $F_{\rm
rad}^{j}$ iteratively. Namely,
(1) we first set $F_{\rm rad}^{j}=0$ and update $du^{j}/d\lambda$ and
$du^{0}/d\lambda=\gamma$ by equation (43),
(2) second, we subsitiute the solved $u^{\mu}$ in equation (45) to update
$F_{\rm rad}^{j}$,
(3) third, use the updated $F_{\rm rad}^{j}$ and solve equation (43) for
$u^{\mu}$ again,
(4) iterate steps (2) and (3) until $F_{\rm rad}^{j}$ saturates. If we use the
$F_{\rm rad}^{j}$ obtained in the previous step for each particle and start
from step (3), the number of iterations can be reduced.

##  Appendix E Discretization of basic equations

We discretize the Maxwell equations (35)–(37) adopting the Yee lattice.
Electric field components, $E$ and $D$, are evaluated at half integer time
steps, e.g., at $t\pm(1/2)\Delta_{t}$, and the magnetic field compoent, $B$,
is evaluated at integer time step at time $t+\Delta_{t}$. Accordingly,
equations (35)–(37) are discretized as follows:

$\displaystyle B^{n+1}{}_{i+1/2,j+1/2}-B^{n}{}_{i+1/2,j+1/2}$ $\displaystyle=$
$\displaystyle-
c_{1}\nu_{x}\left(E^{n+1/2}{}_{i+1,j+1/2}-E^{n+1/2}{}_{i,j+1/2}\right)+c_{2}\nu_{y}\left(D^{n+1/2}{}_{i+1/2,j+1}-D^{n+1/2}{}_{i+1/2,j}\right),$
(E1) $\displaystyle D^{n+1/2}{}_{i+1/2,j}-D^{n-1/2}{}_{i+1/2,j}$
$\displaystyle=$ $\displaystyle
c_{3}\nu_{y}\left(B^{n}{}_{i+1/2,j+1/2}-B^{n}{}_{i+1/2,j-1/2}\right)-4\pi(\Delta_{t})\left.\frac{\Sigma^{2}}{A}\right|_{i+1/2,j}(J^{r})^{n}_{i+1/2,j},$
(E2) $\displaystyle E^{n+1/2}{}_{i,j+1/2}-E^{n-1/2}{}_{i,j+1/2}$
$\displaystyle=$ $\displaystyle-
c_{4}\nu_{x}\left(B^{n}{}_{i+1/2,j+1/2}-B^{n}{}_{i-1/2,j+1/2}\right)-4\pi(\Delta_{t})\left.\frac{\Sigma^{2}\Delta\sin\theta}{A}\right|_{i+1/2,j}(J^{\theta})^{n}_{i,j+1/2},$
(E3)

where the super- and sub-scripts denote the temporal and spatial labels,
respectively; and the Courant numbers are defined by

$\nu_{x}\equiv\frac{\Delta_{t}}{\Delta_{x}},$ (E4)
$\nu_{y}\equiv\frac{\Delta_{t}}{\Delta_{y}}.$ (E5)

Here, $\Delta_{t}$ denotes the time step, $\Delta_{x}$ denotes the constant
interval of the radial tortoise coordinate, and
$\Delta_{y}=\sin\theta\Delta_{\theta}$ denotes the constant azimuthal grid
interval. Denoting the inner and outer boundary positions as $r_{\ast}{}^{\rm
min}$ and $r_{\ast}{}^{\rm max}$, we obtain $\Delta_{x}=(r_{\ast}{}^{\rm
max}-r_{\ast}{}^{\rm min})/(N_{x}-1)$, where $N_{x}$ denotes the number of
radial grids. Since we consider the colatitudes range $0<1-\cos\theta<2$, we
obtain $\Delta_{y}=2/(N_{y}-1)$, where $N_{y}$ denotes the number of
meridional grids.

As for the Hamilton-Jacobi equations (43) and (46)–(48), we evaluate the
particle position ($r$,$\theta$,$\varphi$) at half integer time steps, and the
momentum ($u^{r}$,$u^{\theta}$,$u^{\varphi}$, and $\gamma$) at integer time
steps. Namely, particle’s momentum evolves by (eq. [43])

$\mbox{\boldmath$u$}^{n+1}-\mbox{\boldmath$u$}^{n}=(\Delta_{t})\mbox{\boldmath$F$}(t^{n+1/2},\mbox{\boldmath$x$}^{n+1/2},\mbox{\boldmath$u$}^{n+1/2}),$
(E6)

where $F_{i}$ represents the right-hand side of equation (43).

On the other hand, particle’s position evolves by (eq. [46]–[48])

$r^{n+1/2}=r^{n-1/2}+\left(\frac{u^{r}}{u^{t}}\right)^{n},$ (E7)
$\theta^{n+1/2}=\theta^{n-1/2}+\left(\frac{u^{\theta}}{u^{t}}\right)^{n},$
(E8)
$\varphi^{n+1/2}=\varphi^{n-1/2}+\left(\frac{u^{\varphi}}{u^{t}}\right)^{n}$
(E9)

The Boyer-Lindquist radial coordinate $r$ can be readily converted into
$x=r_{\ast}$ by equation (10), whereas $\theta$ is related to $y$ by
$y=1-\cos\theta$.

To find the initial, stationary solution, we divide the tortoise coordinate
$x$ uniformly between $-15.8180<x<20.0000$, and the meridional coordinate $y$
uniformly between $0<y=1-\cos\theta<2$. Accordingly, the Boyer-Lindquist
coordinates are divided non-uniformly as $r_{1}=1.4651M$, $r_{2}=1.4656M$,
$r_{3}=1.4662M$, $\ldots$, $r_{599}=19.9462M$, $r_{600}=20.0000M$, and
$\theta_{1}=0^{\circ}$, $\theta_{2}=2.6160^{\circ}$,
$\theta_{3}=3.7000^{\circ}$, $\theta=4.5319^{\circ}$, $\ldots$,
$\theta_{959}=89.9104^{\circ}$, $\theta_{960}=89.9701^{\circ}$,
$\theta_{961}=90.0298^{\circ}$, $\ldots$,
$\theta_{1918}=(180.0000-3.7000)^{\circ}$,
$\theta_{1919}=(180.0000-2.6160)^{\circ}$ $\theta_{1920}=180.0000^{\circ}$.

To perform PIC simulations, we must restrict the radial range small compared
to the stationary case, because the Gauss’s and Biot-Savart laws are solved
for the toroidal current distribution within $r<20M$, which makes the solved
poloidal electromagnetic fields be accurate typically only within $r<14M$. The
actual value of the right-hand side (i.e., $14M$ in this case) depends on the
$r$ dependence of $J^{\varphi}{}_{\rm eq}$. For a split-monople case,
$J^{\varphi}{}_{\rm eq}\propto r^{-4}$, or a paraboloidal case,
$J^{\varphi}{}_{\rm eq}\propto r^{-3}$, the outer radial boundary of $r\sim
14M$ is conservatively justified. For PIC simulations, we thus adopt the inner
960 points to restrict the simulation range within
$x_{1}=-15.8180<x<x_{960}=12.8244$, or equivalently,
$r_{1}=1.4646M<r<r_{960}=13.6854M$. For meridional coordinate $y$, and adopt
the same gridding as in the stationary case. We checked that the results
little change when we halve the resolution either in radial or meridional
direction, or when we halve the creation rate $1/k_{\rm create}$ (and hence
the number of particles per cell).

##  Appendix F Current deposit

In our present 2D PIC simulations, we adopt a grid consisting of rectangles of
unit size, and place on it unit rectangular charges, where each unit rectangle
grid cell has an area $\Delta_{x}\Delta_{y}$ (Appendix E). The total amount of
charges is assumed to be the uniformly distributed over its 2D surface. As a
charge moves, each grid cell boundary sweeps a fraction of the particle’s unit
area surface, giving the electric current flowing across the boundary. Using
this ‘area weighting’ (Villasenor & Buneman, 1992), we sum up the currents
carried by individual particles on each grid boundary.

We should notice here that we take account of the particle motion only in the
poloidal plane when we compute the current density, although the particles’
EOM is solved three-dimensionally. Thus, we cannot compute the toroidal
current density, $J^{\varphi}$, in our area weighting method. It follows that
the time evolution of $E_{\varphi}$ cannot be solved from the Ampere’s law.
Accordingly, the evolution of $B^{r}$ or $B^{\theta}$ cannot be solved either,
because $r$ and $\theta$ components of the Faraday law contain
$\partial_{\theta}E_{\varphi}$ and $\partial_{r}E_{\varphi}$, respectively.
Thus, we are able to solve only the three components (i.e., $B$, $D$, and $E$)
of the electromagnetic field by equations (35)–(37) in the present 2D scheme,
unless $J^{\varphi}$ is constructed from particle’s toroidal motion by a
‘volume weighting’ method, which is out of the scope of the present paper.

## References

  * Asada & Nakamura (2012) Asada, K., & Nakamura, M. 2012, ApJ, 745, L28, doi: 10.1088/2041-8205/745/2/L28
  * Asada et al. (2016) Asada, K., Nakamura, M., & Pu, H.-Y. 2016, ApJ, 833, 56, doi: 10.3847/1538-4357/833/1/56
  * Bacchini et al. (2018) Bacchini, F., Ripperda, B., Chen, A. Y., & Sironi, L. 2018, ApJS, 237, 6, doi: 10.3847/1538-4365/aac9ca
  * Bardeen (1970) Bardeen, J. M. 1970, Nature, 226, 64, doi: 10.1038/226064a0
  * Beskin et al. (1992) Beskin, V. S., Istomin, Y. N., & Parev, V. I. 1992, Soviet Ast., 36, 642
  * Biretta et al. (1999) Biretta, J. A., Sparks, W. B., & Macchetto, F. 1999, ApJ, 520, 621, doi: 10.1086/307499
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433, doi: 10.1093/mnras/179.3.433
  * Bohdan et al. (2019) Bohdan, A., Niemiec, J., Pohl, M., et al. 2019, ApJ, 885, 10, doi: 10.3847/1538-4357/ab43cf
  * Boyer & Lindquist (1967) Boyer, R. H., & Lindquist, R. W. 1967, Journal of Mathematical Physics, 8, 265, doi: 10.1063/1.1705193
  * Camenzind (1986) Camenzind, M. 1986, A&A, 162, 32
  * Cerutti et al. (2012) Cerutti, B., Uzdensky, D. A., & Begelman, M. C. 2012, ApJ, 746, 148, doi: 10.1088/0004-637X/746/2/148
  * Cerutti et al. (2013) Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147, doi: 10.1088/0004-637X/770/2/147
  * Chen & Yuan (2020) Chen, A. Y., & Yuan, Y. 2020, ApJ, 895, 121, doi: 10.3847/1538-4357/ab8c46
  * Chen et al. (2018) Chen, A. Y., Yuan, Y., & Yang, H. 2018, ApJ, 863, L31, doi: 10.3847/2041-8213/aad8ab
  * Crinquand et al. (2020) Crinquand, B., Cerutti, B., Philippov, A. e., Parfrey, K., & Dubus, G. 2020, Phys. Rev. Lett., 124, 145101, doi: 10.1103/PhysRevLett.124.145101
  * Event Horizon Telescope Collaboration et al. (2019) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019, ApJ, 875, L5, doi: 10.3847/2041-8213/ab0f43
  * Ferrarese et al. (1996) Ferrarese, L., Ford, H. C., & Jaffe, W. 1996, ApJ, 470, 444, doi: 10.1086/177876
  * Hada et al. (2011) Hada, K., Doi, A., Kino, M., et al. 2011, Nature, 477, 185, doi: 10.1038/nature10387
  * Hada et al. (2013) Hada, K., Kino, M., Doi, A., et al. 2013, ApJ, 775, 70, doi: 10.1088/0004-637X/775/1/70
  * Hada et al. (2016) —. 2016, ApJ, 817, 131, doi: 10.3847/0004-637X/817/2/131
  * Hirose et al. (2004) Hirose, S., Krolik, J. H., De Villiers, J.-P., & Hawley, J. F. 2004, ApJ, 606, 1083, doi: 10.1086/383184
  * Hirotani & Okamoto (1998) Hirotani, K., & Okamoto, I. 1998, ApJ, 497, 563, doi: 10.1086/305479
  * Hirotani & Pu (2016) Hirotani, K., & Pu, H.-Y. 2016, ApJ, 818, 50, doi: 10.3847/0004-637X/818/1/50
  * Hirotani et al. (2016) Hirotani, K., Pu, H.-Y., Lin, L. C.-C., et al. 2016, ApJ, 833, 142, doi: 10.3847/1538-4357/833/2/142
  * Hughes et al. (1994) Hughes, S. A., Keeton, Charles R., I., Walker, P., et al. 1994, Phys. Rev. D, 49, 4004, doi: 10.1103/PhysRevD.49.4004
  * Ichimaru (1977) Ichimaru, S. 1977, ApJ, 214, 840, doi: 10.1086/155314
  * Jackson (1962) Jackson, J. D. 1962, Classical Electrodynamics
  * Junor et al. (1999) Junor, W., Biretta, J. A., & Livio, M. 1999, Nature, 401, 891, doi: 10.1038/44780
  * Kisaka et al. (2020) Kisaka, S., Levinson, A., & Toma, K. 2020, arXiv e-prints, arXiv:2007.02838.  https://arxiv.org/abs/2007.02838
  * Koide et al. (2002) Koide, S., Shibata, K., Kudoh, T., & Meier, D. L. 2002, Science, 295, 1688, doi: 10.1126/science.1068240
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Kovalev et al. (2007) Kovalev, Y. Y., Lister, M. L., Homan, D. C., & Kellermann, K. I. 2007, ApJ, 668, L27, doi: 10.1086/522603
  * Larkin & McLaughlin (2016) Larkin, A. C., & McLaughlin, D. E. 2016, MNRAS, 462, 1864, doi: 10.1093/mnras/stw1749
  * Levinson & Cerutti (2018) Levinson, A., & Cerutti, B. 2018, A&A, 616, A184, doi: 10.1051/0004-6361/201832915
  * Levinson & Rieger (2011) Levinson, A., & Rieger, F. 2011, ApJ, 730, 123, doi: 10.1088/0004-637X/730/2/123
  * Ly et al. (2007) Ly, C., Walker, R. C., & Junor, W. 2007, ApJ, 660, 200, doi: 10.1086/512846
  * Mahadevan (1997) Mahadevan, R. 1997, ApJ, 477, 585, doi: 10.1086/303727
  * McKinney & Narayan (2007) McKinney, J. C., & Narayan, R. 2007, MNRAS, 375, 513, doi: 10.1111/j.1365-2966.2006.11301.x
  * McKinney et al. (2012) McKinney, J. C., Tchekhovskoy, A., & Bland ford, R. D. 2012, MNRAS, 423, 3083, doi: 10.1111/j.1365-2966.2012.21074.x
  * Mertens et al. (2016) Mertens, F., Lobanov, A. P., Walker, R. C., & Hardee, P. E. 2016, A&A, 595, A54, doi: 10.1051/0004-6361/201628829
  * Nakamura & Asada (2013) Nakamura, M., & Asada, K. 2013, ApJ, 775, 118, doi: 10.1088/0004-637X/775/2/118
  * Nakamura & Meier (2014) Nakamura, M., & Meier, D. L. 2014, ApJ, 785, 152, doi: 10.1088/0004-637X/785/2/152
  * Nakamura et al. (2018) Nakamura, M., Asada, K., Hada, K., et al. 2018, ApJ, 868, 146, doi: 10.3847/1538-4357/aaeb2d
  * Narayan & Yi (1994) Narayan, R., & Yi, I. 1994, ApJ, 428, L13, doi: 10.1086/187381
  * Neronov & Aharonian (2007) Neronov, A., & Aharonian, F. A. 2007, ApJ, 671, 85, doi: 10.1086/522199
  * Parfrey et al. (2019) Parfrey, K., Philippov, A., & Cerutti, B. 2019, Phys. Rev. Lett., 122, 035101, doi: 10.1103/PhysRevLett.122.035101
  * Perlman et al. (2011) Perlman, E. S., Adams, S. C., Cara, M., et al. 2011, ApJ, 743, 119, doi: 10.1088/0004-637X/743/2/119
  * Tanabe & Nagataki (2008) Tanabe, K., & Nagataki, S. 2008, Phys. Rev. D, 78, 024004, doi: 10.1103/PhysRevD.78.024004
  * Tchekhovskoy et al. (2010) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2010, ApJ, 711, 50, doi: 10.1088/0004-637X/711/1/50
  * Thorne & MacDonald (1982) Thorne, K. S., & MacDonald, D. 1982, MNRAS, 198, 339, doi: 10.1093/mnras/198.2.339
  * Villasenor & Buneman (1992) Villasenor, J., & Buneman, O. 1992, Computer Physics Communications, 69, 306, doi: 10.1016/0010-4655(92)90169-Y
  * Walker et al. (2018) Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128, doi: 10.3847/1538-4357/aaafcc
  * Wang & Zhou (2009) Wang, C.-C., & Zhou, H.-Y. 2009, MNRAS, 395, 301, doi: 10.1111/j.1365-2966.2009.14463.x
  * Yuan & Narayan (2014) Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529, doi: 10.1146/annurev-astro-082812-141003