Analytical and Numerical Methods and Test Calculations of One-dimensional Force-free Magnetodynamics on Arbitrary Magnetic Surfaces across Horizons of Spinning Black Holes

The FFMD equations are based on the Maxwell equations with force-free conditions (Komissarov 2001, 2004). Here we assume that the magnetic field is so strong that we can neglect the plasma inertia. We also assume that the electric resistivity and self gravity of the electromagnetic field vanish. The line element in the spacetime x^μ = (t, x¹, x², x³) is written by ds² = g_μνdx^μdx^ν. Here Greek subscripts like μ and ν run from 0 to 3, while Roman subscripts like i and j run from 1 to 3. We employ the natural unit system, where the speed of light c, gravitational constant G, and magnetic permeability and electric permittivity in the vacuum μ₀, ₀ are unity. We also often set the black hole mass unity, M = 1. The covariant forms of the Maxwell equations are

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{F}^{\mu \nu }=\displaystyle \frac{1}{\sqrt{-g}}\displaystyle \frac{\partial }{\partial {x}^{\mu }}\left(\sqrt{-g}{F}^{\mu \nu }\right)=-{J}^{\nu },\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }^{* }{F}^{\mu \nu }=\displaystyle \frac{1}{\sqrt{-g}}{\partial }_{\mu }({\sqrt{-g}}^{* }{F}^{\mu \nu })=0,\end{eqnarray} \tag{ 2 }$$

where ∇_μ is the covariant derivative, F_μν is the electromagnetic field tensor, ${}^{* }{F}_{\mu \nu }={\epsilon }^{\mu \nu \rho \sigma }{F}_{\rho \sigma }/2$ is the dual tensor of F_μν, and J^μ = (ρ_e, J¹, J², J³) is the four-current density (ρ_e is the electric charge density). Here ${\epsilon }^{\mu \nu \rho \sigma }={\eta }^{\mu \nu \rho \sigma }/\sqrt{-g}$ is the Levi–Civita tensor, g = det(g_μν) is the determinant of (g_μν), and η^μνρσ is the totally asymmetric symbol.³ We use the electric field E_μ and magnetic field B^μ given by

$$\begin{eqnarray}&&\begin{array}{l}{E}_{\mu }={F}_{\mu 0},\ {B}^{\mu }=\,* {F}^{0\mu }=\displaystyle \frac{1}{2}{\epsilon }^{0\mu \rho \sigma }{F}_{\rho \sigma }.\end{array}\end{eqnarray} \tag{ 3 }$$

The Maxwell equations read the conservation law of the electromagnetic energy and momentum,

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{T}^{\mu \nu }=-{f}_{\nu }^{{\rm{L}}},\end{eqnarray} \tag{ 4 }$$

where ${T}^{\mu \nu }={{F}^{\mu }}_{\sigma }{F}^{\nu \sigma }$ − $\tfrac{1}{4}{g}^{\mu \nu }{F}^{\lambda \kappa }{F}_{\lambda \kappa }$ is the four-electromagnetic energy-momentum tensor and f^ν_L = J_μF^μν is the four-Lorentz force density. In the case of the force-free condition ${f}_{\nu }^{{\rm{L}}}$ = 0, Equation (4) yields

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{T}^{\mu \nu }=0.\end{eqnarray} \tag{ 5 }$$

To perform FFMD numerical simulations, we use Equation (4) instead of Equation (1)

The degeneracy of the force-free electromagnetic field is imposed on the FFMD supplementarily. This degeneracy is due to the condition of the magnetospheric plasma, where the charged particles are plentiful enough to support a strong electric current and screen the electric field (Goldreich & Julian 1969). To derive the degeneracy, we consider the low-inertia limit of relativistic magnetohydrodynamics (RMHD). The relativistic Ohm's law is

$$\begin{eqnarray}&&{u}^{\nu }{F}_{\mu \nu }=\eta ({J}_{\mu }-{\rho }_{{\rm{e}}}^{{\prime} }{u}_{\mu }),\end{eqnarray} \tag{ 6 }$$

where u^μ is the four-velocity of the plasma, η is the resistivity, and ${\rho }_{{\rm{e}}}^{{\prime} }$ is the proper electric density. Equation (6) yields

$$\begin{eqnarray}&&{E}_{\tilde{i}}+{\epsilon }_{{ijk}}{v}^{\tilde{j}}{B}^{\tilde{k}}={\displaystyle \frac{\eta }{\gamma }}_{{\rm{L}}}({J}_{\tilde{i}}-{\rho }_{{\rm{e}}}^{{\prime} }{u}_{\tilde{i}})=\eta {\acute{J}}_{\tilde{i}},\end{eqnarray} \tag{ 7 }$$

where ${\gamma }_{{\rm{L}}}={u}^{\tilde{0}}$, ${v}^{\tilde{i}}={u}^{\tilde{i}}/{\gamma }_{{\rm{L}}}$, and ${\acute{J}}_{\tilde{i}}\equiv \tfrac{1}{{\gamma }_{{\rm{L}}}}({J}_{\tilde{i}}$ − ${\rho }_{{\rm{e}}}^{{\prime} }{u}_{\tilde{i}})$ is similar to the current density observed by the plasma rest frame but not exactly. We have

$$\begin{eqnarray}\begin{array}{rcl}{}^{* }{F}^{\rho \sigma }{F}_{\rho \sigma } & = & {}^{* }{F}^{\tilde{\rho }\tilde{\sigma }}{F}_{\tilde{\rho }\tilde{\sigma }}=-4{E}_{\tilde{i}}{B}^{\tilde{i}}\\ & = & \ -4\displaystyle \frac{\eta }{{\gamma }_{{\rm{L}}}}({J}_{\tilde{i}}-{\rho }_{{\rm{e}}}^{{\prime} }{u}_{\tilde{i}}){B}^{\tilde{i}},\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\begin{array}{l}{F}^{\rho \sigma }{F}_{\rho \sigma }={F}^{\tilde{\rho }\tilde{\sigma }}{F}_{\tilde{\rho }\tilde{\sigma }}=2({B}^{\tilde{i}}{B}_{\tilde{i}}-{E}^{\tilde{i}}{E}_{\tilde{i}})\\ \quad =\ 2\left[\displaystyle \frac{1}{{\gamma }_{{\rm{L}}}^{2}}{B}_{\tilde{i}}{B}^{\tilde{i}}+{\left({v}_{\tilde{i}}{B}^{\tilde{i}}\right)}^{2}-\eta \left\{\eta {J}_{\tilde{i}}{J}^{\tilde{i}}+\displaystyle \frac{2}{{\gamma }_{{\rm{L}}}}{\epsilon }_{{ijk}}{v}^{\tilde{i}}{J}^{\tilde{j}}{B}^{\tilde{k}}\right\}\right].\end{array}\end{eqnarray} \tag{ 9 }$$

When η vanishes, Equations (8) and (9) yield the degeneracy

$$\begin{eqnarray}&&{}^{* }{F}^{\rho \sigma }{F}_{\rho \sigma }=0,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{F}^{\rho \sigma }{F}_{\rho \sigma }\gt 0.\end{eqnarray} \tag{ 11 }$$

In the case of the resistive plasma, the degeneracy is not always guaranteed. Because, for example, when η and ${J}_{\tilde{i}}{J}^{\tilde{i}}$ are finite and ${u}_{\tilde{i}}$ vanishes, ${}^{* }\,{F}^{\rho \sigma }{F}_{\rho \sigma }\ne 0$ from Equation (8). Furthermore, when ${B}^{\tilde{i}}$ vanishes, Equation (9) yields ${F}^{\rho \sigma }{F}_{\rho \sigma }=-2{\eta }^{2}{\acute{J}}^{2}\lt 0$ where the degeneracy is broken.

We review the 3+1 formalism of the FFMD equations derived from the covariant Equations (1), (2), and (4). In order to introduce the 3+1 formalism, we use the local coordinate frame called the "normal observer frame" ${x}^{\tilde{\mu }}$, in which the line element is written by

$$\begin{eqnarray}&&{{ds}}^{2}=-d{\tilde{t}}^{2}+{\gamma }_{{ij}}{{dx}}^{\tilde{i}}{{dx}}^{\tilde{j}}.\end{eqnarray} \tag{ 12 }$$

Here we treat γ_ij = g_ij as the elements of a 3 × 3 matrix (γ_ij) and γ = det(g_ij). Then γ^ij indicates the elements of the inverse of the matrix (g_ij), that is, ${\gamma }^{{ik}}{\gamma }_{{kj}}={\delta }_{j}^{i}$. When we define the lapse function α and shift vector βⁱ as g_i0 = g_0i = g_ijβ^j (or β_i = g_i0), α² = −g₀₀ + g_ijβⁱβ^j, we write the line element as

$$\begin{eqnarray}&&{{ds}}^{2}=-{\alpha }^{2}+{g}_{{ij}}({{dx}}^{i}+{\beta }^{i}{dt})({{dx}}^{j}+{\beta }^{j}{dt}).\end{eqnarray} \tag{ 13 }$$

Comparing Equations (12) and (13), we obtain $d\tilde{t}=\alpha {dt}$, ${{dx}}^{\tilde{i}}={{dx}}^{i}+{\beta }^{i}{dt}$, and we have ${\partial }_{\tilde{t}}={\alpha }^{-1}({\partial }_{t}-{\beta }^{i}{\partial }_{i})$, ${\partial }_{\tilde{i}}={\partial }_{i}$.

This local coordinate system is not always orthonormal in the space. Then the components of the vectors and tensors in the local reference frame are not intuitive in general. In the Boyer–Lindquist (BL) coordinates, the local frame of space is already orthogonal, g_ij = ${h}_{i}^{2}{\delta }_{{ij}}$, and then we can orthonormalize easily as

$$\begin{eqnarray}\begin{array}{rcl}{{ds}}^{2} & = & {\eta }_{\mu \nu }{{dx}}^{\hat{\mu }}{{dx}}^{\hat{\nu }}=-d{\hat{t}}^{2}+\displaystyle \sum _{i=1}^{3}{\left({{dx}}^{\hat{i}}\right)}^{2},\\ & = & \ -d{\tilde{t}}^{2}+\displaystyle \sum _{i=1}^{3}{h}_{i}{h}_{j}{\delta }_{{ij}}{{dx}}^{\tilde{i}}{{dx}}^{\tilde{j}},\end{array}\end{eqnarray} \tag{ 14 }$$

where we set $d\hat{t}=d\tilde{t}$, ${{dx}}^{\hat{i}}={h}_{i}{{dx}}^{\tilde{i}}$. We have ${\partial }_{\hat{t}}={\partial }_{\tilde{t}}$, ${\partial }_{\hat{i}}={h}_{i}^{-1}{\partial }_{\tilde{i}}$. Here η_μν is the metric of Minkowski spacetime. This orthonormal local reference frame ${x}^{\hat{\mu }}$ is called the zero angular momentum observer (ZAMO) frame.

The 3+1 formalism of the Maxwell equations is

$$\begin{eqnarray}&&\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial {x}^{i}}(\sqrt{\gamma }{B}^{\tilde{i}})=0,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{B}^{\tilde{i}}=-{\epsilon }^{{ijk}}\displaystyle \frac{\partial }{\partial {x}^{j}}\left[\alpha ({E}_{\tilde{k}}-{\epsilon }_{{kpq}}{N}^{p}{B}^{\tilde{q}})\right],\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial {x}^{i}}(\sqrt{\gamma }{E}^{\tilde{i}})={\tilde{\rho }}_{{\rm{e}}},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{E}^{\tilde{i}}+\alpha ({J}^{\tilde{i}}+{\tilde{\rho }}_{{\rm{e}}}{N}^{i})={\epsilon }^{{ijk}}\displaystyle \frac{\partial }{\partial {x}^{j}}[\alpha ({B}_{\tilde{k}}+{\epsilon }_{{kmn}}{N}^{m}{E}^{\tilde{n}})],\end{eqnarray} \tag{ 18 }$$

where ^ijk = α^0ijk = γ^−1/2η^0ijk is the three-dimensional Levi–Civita tensor, Nⁱ = −βⁱ/α. Note that N^μ = (1/α, Nⁱ) is the four-velocity of the normal observer.

The 3+1 formalism of the force-free condition is

$$\begin{eqnarray}&&{\tilde{J}}^{i}{\tilde{E}}_{i}=0,\,{\widetilde{\rho }}_{{\rm{e}}}{\tilde{E}}_{i}+{\epsilon }_{ijk}{\tilde{J}}^{j}{\tilde{B}}^{k}=0.\end{eqnarray} \tag{ 19 }$$

The conservation laws of energy and momentum of the electromagnetic field are given by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\tilde{u}+\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial {x}^{i}}[\alpha \sqrt{\gamma }({\tilde{S}}^{i}+{N}^{i}\tilde{u})]+\displaystyle \frac{\partial \alpha }{\partial {x}^{i}}{\tilde{S}}^{i}\\ \ +\ \left[{\gamma }_{{jk}}\displaystyle \frac{\partial }{\partial {x}^{i}}(\alpha {N}^{k})+\displaystyle \frac{1}{2}\alpha {N}^{k}\displaystyle \frac{\partial }{\partial {x}^{k}}{\gamma }_{{ij}}\right]{T}^{\tilde{i}\tilde{j}}=-\alpha {\tilde{J}}^{i}{\tilde{E}}_{i},\end{array}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}{\tilde{S}}_{i}+\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial {x}^{k}}[\alpha \sqrt{\gamma }({T}_{\tilde{i}}^{\tilde{k}}+{N}^{k}{\tilde{S}}_{i})]+\displaystyle \frac{\partial \alpha }{\partial {x}^{i}}\tilde{u}\\ \ +\ \displaystyle \frac{\partial }{\partial {x}^{i}}(\alpha {N}^{k}){\tilde{S}}_{k}-\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {x}^{i}}{\gamma }_{{jk}}{T}^{\tilde{j}\tilde{k}}=-\alpha ({\tilde{\rho }}_{{\rm{e}}}{\tilde{E}}_{i}+{\epsilon }_{{ijk}}{\tilde{J}}^{j}{\tilde{B}}^{k}),\end{array}\end{eqnarray} \tag{ 21 }$$

where ${S}_{\tilde{i}}={\epsilon }_{{ijk}}{E}^{\tilde{j}}{B}^{\tilde{k}}$, $\tilde{u}=({E}^{\tilde{i}}{E}_{\tilde{i}}+{B}^{\tilde{i}}{B}_{\tilde{i}})/2$, ${B}_{\tilde{i}}={\gamma }_{{ij}}{B}^{\tilde{j}}$, and ${E}^{\tilde{i}}={\gamma }^{{ij}}{E}_{\tilde{j}}$ (Appendix A).

Here we derive the equation of energy conservation from Equation (4) with the Killing vector of temporal boost transformation ξ^μ = (1, 0, 0, 0). Using the Killing equation ∇_μξ_ν + ∇_νξ_μ = 0, we obtain

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }({T}^{\mu \nu }{\xi }_{\nu })=\displaystyle \frac{1}{\sqrt{-g}}{\partial }_{\mu }(\sqrt{-g}{T}^{\mu \nu }{\xi }_{\nu })=0.\end{eqnarray} \tag{ 22 }$$

When we use the energy flux density S^μ = T^μνξ_ν and the energy-at-infinity density ${e}^{\infty }\equiv {S}^{\tilde{0}}$, we have the equation of energy conservation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{e}^{\infty }+\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial {x}^{i}}(\alpha \sqrt{\gamma }{S}^{i})=0.\end{eqnarray} \tag{ 23 }$$

Here we introduce flux coordinates for a given axisymmetric magnetic surface along an arbitrary axisymmetric surface imaged by Figure 1. Using the vector potential A^μ for the given axisymmetric electromagnetic field, the magnetic surface is described by Ψ = A_ϕ = constant. The flux coordinates are given by x^μ = (t, r, Ψ, ϕ). Here t, r, and ϕ are the time, radial, and azimuthal coordinates, respectively. Using the flux coordinates, we have B^Ψ = 0. For a radial magnetic surface as shown in Figure 2, we use the colatitude coordinate θ as the coordinate Ψ − Ψ₀ + θ₀, where the radial magnetic surface is given by θ = θ₀. Note that, in this case, we use the unit system so that $\tfrac{\partial {\rm{\Psi }}}{\partial \theta }=1$. In general, the flux coordinates are not orthogonal, while the coordinates are orthogonal at the equatorial plane.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When we use the flux coordinates (r, Ψ, ϕ) (${B}^{\tilde{{\rm{\Psi }}}}=0$), we have the following three conditions with respect to ${E}_{\tilde{\phi }}$, ${E}_{\tilde{r}}$, and ${J}^{\tilde{{\rm{\Psi }}}}$ in general. It is noted that the BL and Kerr–Schild (KS) coordinates are used for the flux coordinates of radial magnetic surfaces around the equatorial plane. In the case of a radial magnetic surface, Ψ is directly given by θ without the addition of Ψ₀ − θ₀.

The Ψ-component of Equation (16) yields

$$\begin{eqnarray*}&&\displaystyle \frac{\partial {B}^{\tilde{{\rm{\Psi }}}}}{\partial t}=-\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}[\alpha ({E}_{\tilde{\phi }}+\sqrt{\gamma }{N}^{{\rm{\Psi }}}{B}^{\tilde{r}})]=0,\end{eqnarray*}$$

and we find that $\alpha ({E}_{\tilde{\phi }}+\sqrt{\gamma }{N}^{{\rm{\Psi }}}{B}^{\tilde{r}})$ is uniform. Then we have

$$\begin{eqnarray}&&{E}_{\tilde{\phi }}=-\sqrt{\gamma }{N}^{{\rm{\Psi }}}{B}^{\tilde{r}}\end{eqnarray} \tag{ 24 }$$

because $\alpha {E}_{\tilde{\phi }}+\sqrt{\gamma }{N}^{{\rm{\Psi }}}{B}^{\tilde{r}}$ should vanish at infinity. The azimuthal component of the force-free condition (Equation (19)) reads

$$\begin{eqnarray*}&&{\tilde{\rho }}_{{\rm{e}}}{E}_{\tilde{\phi }}+{\epsilon }_{\phi {\rm{\Psi }}r}{J}^{\tilde{{\rm{\Psi }}}}{B}^{\tilde{r}}=-\sqrt{\gamma }({J}^{\tilde{{\rm{\Psi }}}}+{\tilde{\rho }}_{{\rm{e}}}{N}^{{\rm{\Psi }}}){B}^{\tilde{r}}=0.\end{eqnarray*}$$

Then we have

$$\begin{eqnarray}&&{J}^{\tilde{{\rm{\Psi }}}}=-{\tilde{\rho }}_{{\rm{e}}}{N}^{{\rm{\Psi }}}\end{eqnarray} \tag{ 25 }$$

because ${B}^{\tilde{r}}$ is finite.

The degeneracy of the field in FFMD (Equation (10)) and Equation (24) yields ${E}_{\tilde{i}}{B}^{\tilde{i}}={E}_{\tilde{r}}{B}^{\tilde{r}}+{E}_{\tilde{\phi }}{B}^{\tilde{\phi }}$ = ${B}^{\tilde{r}}({E}_{\tilde{r}}-\sqrt{\gamma }{N}^{{\rm{\Psi }}}{B}^{\tilde{\phi }})\,=0$. Then we have

$$\begin{eqnarray}&&{E}_{\tilde{r}}=\sqrt{\gamma }{N}^{{\rm{\Psi }}}{B}^{\tilde{\phi }}.\end{eqnarray} \tag{ 26 }$$

The relationship given by Equations (24)–(26) is important for the following calculations.

We derive a general form of 1D FFMD equations for an arbitrary axisymmetric magnetic surface of the force-free electromagnetic field. Equations (15)–(18), (21), and (23) yield the following "1D" equations:

$$\begin{eqnarray}&&\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}(\sqrt{\gamma }{B}^{\tilde{r}})=0,\ \displaystyle \frac{\partial }{\partial t}{B}^{\tilde{r}}=0,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{E}^{\tilde{{\rm{\Psi }}}}=-\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}[\alpha ({B}_{\tilde{\phi }}+\sqrt{\gamma }{N}^{r}{E}^{\tilde{{\rm{\Psi }}}}-\sqrt{\gamma }{N}^{{\rm{\Psi }}}{E}^{\tilde{r}})],\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{B}^{\tilde{\phi }}=-\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}\left[\alpha \sqrt{\gamma }({E}_{\tilde{{\rm{\Psi }}}}-\sqrt{\gamma }{N}^{\phi }{B}^{\tilde{r}}+\sqrt{\gamma }{N}^{r}{B}^{\tilde{\phi }})\right],\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}{S}_{\tilde{r}}=-\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}[\alpha \sqrt{\gamma }({T}_{\tilde{r}}^{\tilde{r}}+{N}^{r}{\tilde{S}}_{r})]-\displaystyle \frac{\partial \alpha }{\partial r}\tilde{u}\\ -\ \displaystyle \frac{\partial }{\partial r}(\alpha {N}^{r}){S}_{\tilde{r}}-\displaystyle \frac{\partial }{\partial r}(\alpha {N}^{\phi }){S}_{\tilde{\phi }}-\displaystyle \frac{\partial }{\partial {\rm{\Psi }}}(\alpha {N}^{{\rm{\Psi }}}){S}_{\tilde{{\rm{\Psi }}}}\\ \qquad +\ \displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial r}{\gamma }_{{jk}}{T}^{\tilde{j}\tilde{k}},\end{array}\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{S}_{\tilde{\phi }}=-\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}[\alpha \sqrt{\gamma }({T}_{\tilde{\phi }}^{\tilde{r}}+{N}^{r}{S}_{\tilde{\phi }})],\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{e}^{\infty }=-\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}(\alpha \sqrt{\gamma }{S}^{r}).\end{eqnarray} \tag{ 32 }$$

Note that other equations from Equations (17), (18), (20), (21), and (23) yield non-1D equations that include the Ψ-derivative as shown in Appendix B. Equation (27) shows $\sqrt{\gamma }{B}^{\tilde{r}}$ is uniform and constant. In fact, we have ${B}^{\tilde{r}}{=}^{* }{F}^{\tilde{0}\tilde{r}}\,=\tfrac{1}{\sqrt{\gamma }}({\partial }_{{\rm{\Psi }}}{A}_{\phi }-{\partial }_{\phi }{A}_{{\rm{\Psi }}})=\tfrac{1}{\sqrt{\gamma }}$. The poloidal component of the magnetic field ${\tilde{B}}_{{\rm{p}}}$ is defined by ${\tilde{B}}_{{\rm{p}}}=\sqrt{{\gamma }_{{rr}}}{B}^{\tilde{r}}=\sqrt{{\gamma }_{{rr}}/\gamma }$.

The primitive variable ${E}^{\tilde{{\rm{\Psi }}}}$ is given by conservation quantities ${B}^{\tilde{r}}$, ${B}^{\tilde{\phi }}$, ${S}_{\tilde{r}}$, and ${S}_{\tilde{\phi }}$ as

$$\begin{eqnarray}&&{E}^{\tilde{{\rm{\Psi }}}}=\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{{S}_{\tilde{r}}{B}^{\tilde{\phi }}-{S}_{\tilde{\phi }}{B}^{\tilde{r}}}{{\left(\tilde{B}\right)}^{2}},\end{eqnarray} \tag{ 33 }$$

where ${B}_{\tilde{r}}={\gamma }_{{rr}}{B}^{\tilde{r}}+{\gamma }_{r\phi }{B}^{\tilde{\phi }}$, ${B}_{\tilde{\phi }}={\gamma }_{\phi r}{B}^{\tilde{r}}+{\gamma }_{\phi \phi }{B}^{\tilde{\phi }}$, and ${(\tilde{B})}^{2}={B}^{\tilde{r}}{B}_{\tilde{r}}+{B}^{\tilde{\theta }}{B}_{\tilde{\theta }}$. We have ${E}_{\tilde{{\rm{\Psi }}}}=\tfrac{1}{{\gamma }^{{\rm{\Psi }}{\rm{\Psi }}}}({E}^{\tilde{{\rm{\Psi }}}}-\sqrt{\gamma }{N}^{{\rm{\Psi }}}{\gamma }^{{\rm{\Psi }}r}{B}^{\tilde{\phi }})$ using ${E}^{\tilde{{\rm{\Psi }}}}={\gamma }^{{\rm{\Psi }}i}{E}_{\tilde{i}}$ with Equations (24) and (26) and γ^Ψϕ = 0. We can write

$$\begin{eqnarray}&&{T}_{\tilde{j}}^{\tilde{i}}=-{E}^{\tilde{i}}{E}_{\tilde{j}}-{B}^{\tilde{i}}{B}_{\tilde{j}}+\tilde{u}{\delta }_{j}^{i},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\tilde{u}=\displaystyle \frac{1}{2}({B}^{\tilde{r}}{B}_{\tilde{r}}+{B}^{\tilde{\phi }}{B}_{\tilde{\phi }}+{\left(\tilde{E}\right)}^{2}),\end{eqnarray} \tag{ 35 }$$

where ${(\tilde{E})}^{2}={E}^{\tilde{i}}{E}_{\tilde{i}}$.

Using Equations (29)–(31) and (33)–(35) with ${B}^{\tilde{r}}={B}_{0}/\sqrt{\gamma }$ (B₀ is constant), we perform 1D FFMD numerical simulations for an arbitrary magnetic surface. When we use the flux coordinates, we have B₀ = 1.

Here we show the constants for the steady state in 1D FFMD. We assume that the force-free field is stationary and axisymmetric. First, the equation of energy conservation (Equation (32)) reads the energy flux constant,

$$\begin{eqnarray}&&P=\alpha \sqrt{\gamma }{S}^{r}.\end{eqnarray} \tag{ 36 }$$

Here ${e}^{\infty }=-{\chi }_{\nu }{T}^{\nu 0}=\alpha (\tilde{u}+{N}^{i}{S}_{\tilde{i}})$, ${S}^{r}=-{\chi }_{\nu }{T}^{\nu r}=\tfrac{1}{\alpha }{\epsilon }^{{rjk}}{E}_{j}{B}^{k}\,=\tfrac{1}{\alpha \sqrt{\gamma }}{E}_{{\rm{\Psi }}}{B}_{\phi }$, where χ^μ = (1, 0, 0, 0) is the Killing vector with respect to the time boost symmetry and ${E}_{{\rm{\Psi }}}=\alpha ({E}_{\tilde{{\rm{\Psi }}}}+{\epsilon }_{{\rm{\Psi }}{jk}}{N}^{j}{B}^{\tilde{k}})$ =$\alpha ({E}_{\tilde{{\rm{\Psi }}}}+\sqrt{\gamma }{N}^{r}{B}^{\tilde{\phi }})$, ${B}_{\phi }=\alpha ({B}_{\tilde{\phi }}-{\epsilon }_{\phi {jk}}{N}^{j}{B}^{\tilde{k}})\,$=$\,\alpha ({B}_{\tilde{\phi }}+\sqrt{\gamma }{N}^{r}{E}^{\tilde{{\rm{\Psi }}}})$. We have

$$\begin{eqnarray}&&P={E}_{{\rm{\Psi }}}{B}_{\phi }={\alpha }^{2}({E}_{\tilde{{\rm{\Psi }}}}+\sqrt{\gamma }{N}^{r}{B}^{\tilde{\phi }})({B}_{\tilde{\phi }}+\sqrt{\gamma }{N}^{r}{E}^{\tilde{{\rm{\Psi }}}}).\end{eqnarray} \tag{ 37 }$$

Next, the vanishing right-hand side of Equation (29) provides the constant,

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{F}}}=-\alpha ({E}_{\tilde{{\rm{\Psi }}}}-\sqrt{\gamma }{N}^{\phi }{B}^{\tilde{r}}+\sqrt{\gamma }{N}^{r}{B}^{\tilde{\phi }}).\end{eqnarray} \tag{ 38 }$$

Using the constant $\sqrt{\gamma }{B}^{\tilde{r}}=1$ in the flux coordinates, we also have

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{F}}}=-\displaystyle \frac{\alpha }{\sqrt{\gamma }{B}^{\tilde{r}}}({E}_{\tilde{{\rm{\Psi }}}}-\sqrt{\gamma }{N}^{\phi }{B}^{\tilde{r}}+\sqrt{\gamma }{N}^{r}{B}^{\tilde{\phi }}).\end{eqnarray} \tag{ 39 }$$

Here Ω_F represents the angular velocity of the magnetic field line. Last, Equation (28) in the situation with ${J}^{\tilde{{\rm{\Psi }}}}=-{\tilde{\rho }}_{{\rm{e}}}{N}^{{\rm{\Psi }}}$ yields the constant

$$\begin{eqnarray}&&I=\alpha ({B}_{\tilde{\phi }}+\sqrt{\gamma }{N}^{r}{E}^{\tilde{{\rm{\Psi }}}}-\sqrt{\gamma }{N}^{{\rm{\Psi }}}{E}^{\tilde{r}}).\end{eqnarray} \tag{ 40 }$$

Here I is the current inside of the magnetic surface of Ψ = Ψ₀: Ψ ≤ Ψ₀.

Using Equations (39) and (40), we obtain the analytic solution of the steady-state force-free field with I and Ω_F,

$$\begin{eqnarray}&&{B}^{\tilde{\phi }}=\displaystyle \frac{1}{\alpha \sqrt{\gamma }D}\left[\sqrt{\gamma }I-{\rm{\Upsilon }}{{\rm{\Omega }}}_{{\rm{F}}}-\alpha ({\gamma }_{\phi r}+{\rm{\Upsilon }}{N}^{\phi })\right],\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{E}_{\tilde{{\rm{\Psi }}}}=\displaystyle \frac{1}{\alpha D}\left[{\gamma }_{\phi \phi }{{\rm{\Omega }}}_{{\rm{F}}}-\sqrt{\gamma }{N}^{r}I+\alpha ({\gamma }_{\phi \phi }{N}^{\phi }+{\gamma }_{\phi r}{N}^{r})\right],\end{eqnarray} \tag{ 42 }$$

where ϒ = γ(N^rγ^ΨΨ − N^Ψγ^rΨ) and D = γ_ϕϕ − ϒN^r. The I and Ω_F of the steady state are given in Section 3. In the BL and KS coordinates, we have ϒ = 0 and ${\gamma }_{\mathrm{KS}}{N}_{\mathrm{KS}}^{r}{\gamma }_{\mathrm{KS}}^{\theta \theta }$, respectively.

Equation (37) yields the relation

$$\begin{eqnarray}&&P=-I{{\rm{\Omega }}}_{{\rm{F}}}.\end{eqnarray} \tag{ 43 }$$

In the BL coordinates, these constants become

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{F}}}=\omega -\displaystyle \frac{\alpha {E}_{\hat{\theta }}}{{{RB}}^{\hat{r}}},\ I=\alpha {{RB}}_{\hat{\phi }}.\end{eqnarray} \tag{ 44 }$$

In the KS coordinates, they are

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Omega }}}_{{\rm{F}}} & = & -\displaystyle \frac{\alpha }{\sqrt{\gamma }{B}^{\tilde{r}}}({E}_{\tilde{\theta }}+\sqrt{\gamma }{N}^{r}{B}^{\tilde{\phi }}),\\ I & = & \alpha ({B}_{\tilde{\phi }}+\sqrt{\gamma }{N}^{r}{E}^{\tilde{\theta }}).\end{array}\end{eqnarray} \tag{ 45 }$$

To describe the spacetime around the spinning black hole with mass M and angular momentum J, we use the two coordinate systems, BL and KS, defined below. Here we use a = (J/J_max)M, where J_max = M² is the maximum angular momentum of a spinning black hole with mass M. The metric of the BL coordinates x^μ = (t, r, θ, ϕ) is given by

$$\begin{eqnarray}\begin{array}{rcl}{{ds}}^{2} & = & {g}_{\mu \nu }{{dx}}^{\mu }{{dx}}^{\nu }=-\left(1-\displaystyle \frac{2{Mr}}{{\rho }^{2}}\right){{dt}}^{2}+\displaystyle \frac{{\rho }^{2}}{{\rm{\Delta }}}{{dr}}^{2}\\ & & +\ {\rho }^{2}d{\theta }^{2}+{R}^{2}d{\phi }^{2}-2\times \displaystyle \frac{2{Mar}}{{\rho }^{2}}{\sin }^{2}\theta d\phi {dt},\end{array}\end{eqnarray} \tag{ 46 }$$

where ρ² = r² + a² cos²θ, Δ = r² − 2Mr + a², Σ² = (r² + a²)²−a²Δ sin²θ, ${R}^{2}=\tfrac{{{\rm{\Sigma }}}^{2}}{{\rho }^{2}}{\sin }^{2}\theta$. The horizon is given by Δ = 0, which yields $r=M\pm \sqrt{{M}^{2}-{a}^{2}}\equiv {r}_{{\rm{H}}}$. In this coordinate system, $\alpha =\sqrt{\tfrac{{\rm{\Delta }}{\rho }^{2}}{{{\rm{\Sigma }}}^{2}}}$, ${\beta }^{\phi }=-\tfrac{2{aMr}}{{{\rm{\Sigma }}}^{2}}\equiv -\omega$, g = −ρ⁴ sin²θ.

The metric of the KS coordinates x^μ = (t, r, θ, ϕ) is given by

$$\begin{eqnarray}\begin{array}{rcl}{{ds}}^{2} & = & {g}_{\mu \nu }{{dx}}^{\mu }{{dx}}^{\nu }=-\left(1-\displaystyle \frac{2{Mr}}{{\rho }^{2}}\right){{dt}}^{2}\\ & & +\ \left(1+\displaystyle \frac{2{Mr}}{{\rho }^{2}}\right){{dr}}^{2}+{\rho }^{2}d{\theta }^{2}+{R}^{2}d{\phi }^{2}\\ & & -\ 2\times a\left(1+\displaystyle \frac{2{Mr}}{{\rho }^{2}}\right){\sin }^{2}\theta {drd}\phi +2\\ & & \times \ \displaystyle \frac{2{Mr}}{{\rho }^{2}}{dtdr}-\displaystyle \frac{2{Mar}}{{\rho }^{2}}{\sin }^{2}\theta {dtd}\phi .\end{array}\end{eqnarray} \tag{ 47 }$$

At the horizon r = r_H, all of the metrics of the KS coordinates are finite, while some metric of the BL coordinates, for example, g_rr, is infinity. In the KS coordinates, we have ${\beta }^{r}=\tfrac{2{Mr}}{{\rho }^{2}+2{Mr}},{\beta }^{\theta }=0,{\beta }^{\phi }=0$, γ = ρ² (ρ² + 2Mr) sin²θ, g =−ρ⁴ sin²θ. The BL and KS coordinates are used for the base of the flux coordinates. In the cases of radial magnetic surfaces as shown in Figure 2, we directly use these coordinates as the flux coordinates.

It is noted that the BL ${x}_{{\rm{B}}{\rm{L}}}^{\mu }$ = (t_BL, r_BL, θ_BL, ϕ_BL) and KS ${x}_{\mathrm{KS}}^{\mu }$ = (t_KS, r_KS, θ_KS, ϕ_KS) coordinates are related as

$$\begin{eqnarray}\begin{array}{rcl}{{dt}}_{\mathrm{KS}} & = & {{dt}}_{\mathrm{BL}}+\displaystyle \frac{2{Mr}}{{\rm{\Delta }}}{{dr}}_{\mathrm{BL}},\ d{\phi }_{\mathrm{KS}}=d{\phi }_{\mathrm{BL}}+\displaystyle \frac{a}{{\rm{\Delta }}}{{dr}}_{\mathrm{BL}},\\ {r}_{\mathrm{KS}} & = & {r}_{\mathrm{BL}},\ {\theta }_{\mathrm{KS}}={\theta }_{\mathrm{BL}}.\end{array}\end{eqnarray} \tag{ 48 }$$

Then we note that at the horizon and infinitely far from the black hole, the times on the KS and BL coordinates are infinitely different. That is, the finite time on the BL coordinates corresponds to the infinite past on the KS coordinates. Conversely, at an infinitely far point from the black hole, the finite time on the BL coordinates corresponds to the infinite future on the KS coordinates.

We introduce the flux coordinates x^μ = (t, r, Ψ, ϕ) around an arbitrary axisymmetric magnetic surface Ψ = Ψ₀, where Ψ₀ is a constant, which is described by θ = θ₀(r), where θ₀(r) is function of r. When the magnetic surface Ψ = Ψ around the magnetic surfaces Ψ = Ψ₀ is described by θ = θ(r, Ψ), the Taylor expansion of θ(r, Ψ) with an infinitesimally small variable Ψ − Ψ₀ yields $\theta (r,{\rm{\Psi }})=\theta (r,{{\rm{\Psi }}}_{0})$ + ${\left(\tfrac{\partial \theta (r,{\rm{\Psi }})}{\partial {\rm{\Psi }}}\right)}_{{\rm{\Psi }}={{\rm{\Psi }}}_{0}}$ $({\rm{\Psi }}-{{\rm{\Psi }}}_{0})+\cdots$. Then we have

$$\begin{eqnarray}&&{\rm{\Psi }}\approx {{\rm{\Psi }}}_{0}+b(r)\left(\theta -{\theta }_{0}(r)\right),\end{eqnarray} \tag{ 49 }$$

where $b(r)={\left[{\left(\tfrac{\partial \theta (r,{\rm{\Psi }})}{\partial {\rm{\Psi }}}\right)}_{{\rm{\Psi }}={{\rm{\Psi }}}_{0}}\right]}^{-1}$ and θ₀(r) = θ(r, Ψ₀). The function b(r) determines the flared shape of the magnetic surfaces. The metric of the flux coordinates in the Kerr spacetime is given by the metric of KS coordinates ${x}_{\mathrm{KS}}^{\mu }$,

$$\begin{eqnarray}\begin{array}{rcl}{\gamma }_{{rr}} & = & {\gamma }_{{rr}}^{\mathrm{KS}}+{K}^{2}{\gamma }_{\theta \theta }^{\mathrm{KS}},{\gamma }_{{\rm{\Psi }}{\rm{\Psi }}}=\displaystyle \frac{1}{{b}^{2}}{\gamma }_{\theta \theta }^{\mathrm{KS}},{\gamma }_{\phi \phi }={\gamma }_{\phi \phi }^{\mathrm{KS}},\\ {\gamma }_{r{\rm{\Psi }}} & = & \displaystyle \frac{K}{b}{\gamma }_{\theta \theta }^{\mathrm{KS}},{\gamma }_{{\rm{\Psi }}\phi }=0,{\gamma }_{r\phi }={\gamma }_{r\phi }^{\mathrm{KS}},\\ {\alpha }^{2} & = & {\alpha }_{\mathrm{KS}}^{2},{\beta }^{r}={\beta }_{\mathrm{KS}}^{r},{\beta }^{{\rm{\Psi }}}={bK}{\beta }_{\mathrm{KS}}^{r},{\beta }^{\phi }=0,\end{array}\end{eqnarray} \tag{ 50 }$$

where $K\equiv {\theta }_{0}^{{\prime} }(r)-\tfrac{{b}^{{\prime} }(r)}{b{\left(r\right)}^{2}}$. We also have

$$\begin{eqnarray}\begin{array}{rcl}{\gamma }^{{rr}} & = & {\gamma }_{\mathrm{KS}}^{{rr}},{\gamma }^{{\rm{\Psi }}{\rm{\Psi }}}={b}^{2}({\gamma }_{\mathrm{KS}}^{{\rm{\Psi }}{\rm{\Psi }}}+{K}^{2}{\gamma }^{{rr}}),{\gamma }^{\phi \phi }={\gamma }_{\mathrm{KS}}^{\phi \phi },\\ {\gamma }^{r{\rm{\Psi }}} & = & {Kb}{\gamma }_{\mathrm{KS}}^{{rr}},{\gamma }^{{\rm{\Psi }}\phi }=-{Kb}{\gamma }_{\mathrm{KS}}^{\phi r},{\gamma }^{r\phi }={\gamma }_{\mathrm{KS}}^{r\phi }.\end{array}\end{eqnarray} \tag{ 51 }$$

To show the validity of the 1D FFMD analytic solution, we derive the analytic solutions of the steady-state force-free field around a extremely slowly spinning black hole (a ≪ M) given by Blandford & Znajek (1977). Blandford & Znajek (1977) resorted to a perturbation method in which they expanded on the powers of a/M. Such a technique can only be of use when the change in the poloidal field caused by spinning up a nonrotating field configuration (supported by currents in an equatorial disk) can be regarded as small. They wrote an exact axisymmetric vacuum solution for the magnetic field in a Schwarzschild metric by Ψ (r, θ) = X(r, θ) used as the unperturbed function. According to the perturbation method, the perturbed variables, the electromagnetic angular frequency and the current I are written by ${{\rm{\Omega }}}_{{\rm{F}}}=\tfrac{a}{{M}^{2}}W(r,\theta )$ and $I=\tfrac{a}{{M}^{2}}Y(r,\theta )$, respectively. They showed two examples of the perturbation technique: (a) the radial magnetic field (of opposite polarity in the two hemispheres) and (b) a force-free magnetosphere in which the magnetic field lines lie on paraboloidal surfaces (cutting an equatorial disk). We derive the expression of W and Y using the magnetic surfaces for both cases as follows.

3.1.1. The Case of the Radial Magnetic Surface

The vector potential

$$\begin{eqnarray}&&{\rm{\Psi }}(r,\theta )=X(r,\theta )=-C\cos \theta \,(0\leqslant \theta \leqslant \pi /2)\end{eqnarray} \tag{ 61 }$$

with a constant C is an exact solution of the vacuum Maxwell equations in a Schwarzschild metric, which describes the unperturbed radial magnetic field (Equation (6.1) in Blandford & Znajek 1977). Blandford & Znajek (1977) used this solution as the unperturbed vector potential Ψ of the split monopole field around a extremely slowly spinning black hole. It is noted that in the case of the split monopole, the solution is used except on an equatorial disk containing a toroidal surface current density I_sm = 2C/r². Equation (61) yields $b\,={[\partial \theta (r,{\rm{\Psi }})/\partial {\rm{\Psi }}]}^{-1}=C\sin \theta$. Then we have ${b}_{{\rm{H}}}\sin {\theta }_{{\rm{H}}}\,={b}_{\infty }\sin {\theta }_{\infty }$ because ${\rm{\Psi }}=-C\cos {\theta }_{{\rm{H}}}=-C\cos {\theta }_{\infty }$ reads ${\theta }_{{\rm{H}}}\,={\theta }_{\infty }$. Equation (58) yields ${{\rm{\Omega }}}_{{\rm{F}}}=\tfrac{a}{4{{Mr}}_{{\rm{H}}}-{a}^{2}{\cos }^{2}{\theta }_{{\rm{H}}}}$. Using a ≪ M, we have ${{\rm{\Omega }}}_{{\rm{F}}}=\tfrac{a}{8{M}^{2}}=\tfrac{a}{M}W$. Equation (59) yields $I=-\tfrac{{aC}}{8{M}^{2}}{\sin }^{2}\theta \,=\tfrac{a}{{M}^{2}}Y$. Then we obtain

$$\begin{eqnarray}&&W=\displaystyle \frac{1}{8},Y=-\displaystyle \frac{C}{8}{\sin }^{2}\theta .\end{eqnarray} \tag{ 62 }$$

These W and Y are the same expressions given by Equation (6.5) in Blandford & Znajek (1977) in the radial magnetic surface case.

The power per solid angle from the horizon is given by

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{BZ}} & = & \displaystyle \frac{{dL}}{d{\rm{\Omega }}}=\displaystyle \frac{2\pi d{\rm{\Psi }}}{d{\rm{\Omega }}}\displaystyle \frac{{dL}}{2\pi d{\rm{\Psi }}}=\displaystyle \frac{2\pi d{\rm{\Psi }}}{2\pi \sin \theta d\theta }P\\ & = & \ \displaystyle \frac{{b}_{{\rm{H}}}}{\sin {\theta }_{{\rm{H}}}}P,\end{array}\end{eqnarray} \tag{ 63 }$$

where L is the power radiated by the black hole in the region between the pole and the magnetic surface Ψ and Ω is the solid angle. We have

$$\begin{eqnarray}&&{P}_{\mathrm{BZ}}={\left(\displaystyle \frac{{a}^{2}C}{8{M}^{2}}\right)}^{2}{\sin }^{2}\theta .\end{eqnarray} \tag{ 64 }$$

3.1.2. The Case of the Paraboloidal Magnetic Surface

The vector potential of the paraboloidal magnetic surface is given by Equation (7.1) in Blandford & Znajek (1977) as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Psi }}(r,\theta ) & = & X(r,\theta )=\displaystyle \frac{1}{2}[r(1-\cos \theta )\\ & & +\ 2M(1+\cos \theta )(1-\mathrm{ln}(1+\cos \theta ))].\end{array}\end{eqnarray} \tag{ 65 }$$

Equation (65) yields $b={[\partial \theta (r,{\rm{\Psi }})/\partial {\rm{\Psi }}]}^{-1}\,$=$\,\tfrac{C}{2}\sin \theta [r\,+2M\mathrm{ln}(1+\cos \theta )]$. With a ≪ M, we have

$$\begin{eqnarray*}&&\begin{array}{l}{{\rm{\Psi }}}_{{\rm{H}}}=\tfrac{C}{2}[2M(1-\cos {\theta }_{{\rm{H}}})\\ \ +\ 2M(1+\cos {\theta }_{{\rm{H}}})(1-\mathrm{ln}(1+\cos {\theta }_{{\rm{H}}}))],\\ {{\rm{\Psi }}}_{\infty }=\tfrac{C}{2}\left[\tfrac{1}{2}{r}_{\infty }{\theta }_{\infty }^{2}+4M(1-\mathrm{ln}2)\right].\end{array}\end{eqnarray*}$$

Using ${{\rm{\Psi }}}_{{\rm{H}}}={{\rm{\Psi }}}_{\infty }$, we get $\tfrac{1}{2}{r}_{\infty }{\theta }_{\infty }^{2}=2M[2\mathrm{ln}2-(1+\cos {\theta }_{{\rm{H}}})\mathrm{ln}(1+\cos {\theta }_{{\rm{H}}})]$. We calculate as

$$\begin{eqnarray*}\begin{array}{rcl}{b}_{{\rm{H}}}\sin {\theta }_{{\rm{H}}} & = & {CM}{\sin }^{2}{\theta }_{{\rm{H}}}[1+\mathrm{ln}(1+\cos {\theta }_{{\rm{H}}})],{b}_{\infty }\sin {\theta }_{\infty }\\ & \approx & \ {{Cr}}_{\infty }(1-\cos {\theta }_{\infty })\approx C\displaystyle \frac{1}{2}{r}_{\infty }{\theta }_{\infty }^{2}.\end{array}\end{eqnarray*}$$

Using Equations (58) and (59), we obtain

$$\begin{eqnarray*}&&{{\rm{\Omega }}}_{{\rm{F}}}=\displaystyle \frac{a}{{M}^{2}}\displaystyle \frac{N}{D},I=-\displaystyle \frac{a}{{M}^{2}}\displaystyle \frac{N}{D}{b}_{\infty }\sin {\theta }_{\infty },\end{eqnarray*}$$

where $N=\tfrac{1}{4}{\sin }^{2}{\theta }_{{\rm{H}}}[1+\mathrm{ln}(1+\cos {\theta }_{{\rm{H}}})]$ and $D=4\mathrm{ln}2\,+{\sin }^{2}{\theta }_{{\rm{H}}}$ + $\{{\sin }^{2}{\theta }_{{\rm{H}}}-2(1+\cos {\theta }_{{\rm{H}}})\}$ $\mathrm{ln}(1+\cos {\theta }_{{\rm{H}}})$. Then we have

$$\begin{eqnarray}\begin{array}{rcl}W & = & \displaystyle \frac{N}{D},\\ Y & = & 2{MC}[2\mathrm{ln}2-(1+\cos {\theta }_{{\rm{H}}})\mathrm{ln}(1+\cos {\theta }_{{\rm{H}}})]\displaystyle \frac{N}{D}.\end{array}\end{eqnarray} \tag{ 66 }$$

The expression W is identified with the perturbed solution given by Equation (7.5) in Blandford & Znajek (1977).

Using Equation (63), we obtain the power per unit solid angle,

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{BZ}} & = & 2{\left(\displaystyle \frac{a}{M}\right)}^{2}{C}^{2}\{1+\mathrm{ln}(1+\cos {\theta }_{{\rm{H}}})\}\\ & & \{2\mathrm{ln}2\,-\,(1+\cos {\theta }_{{\rm{H}}})\mathrm{ln}(1+\cos {\theta }_{{\rm{H}}})\}{\left(\displaystyle \frac{N}{D}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 67 }$$

We show the 1D FFMD equations for different types of magnetic surfaces along the equatorial plane for numerical calculations. In the KS coordinates, using x^μ = (t, r, θ, ϕ) = (t, r, Ψ, ϕ), we have the Maxwell equation, and the conservation law of momentum at the equatorial plane as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{B}^{\tilde{\phi }}=-\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}\left[\alpha \sqrt{\gamma }({E}_{\tilde{\theta }}+\sqrt{\gamma }{N}^{r}{B}^{\tilde{\phi }})\right],\end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}{S}_{\tilde{r}} & = & -\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}[\alpha \sqrt{\gamma }({T}_{\tilde{r}}^{\tilde{r}}+{N}^{r}{S}_{\tilde{r}})]-\displaystyle \frac{\partial \alpha }{\partial r}\tilde{u}\\ & & -\ \displaystyle \frac{\partial }{\partial r}(\alpha {N}^{r}){S}_{\tilde{r}}-\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial r}{\gamma }_{{jk}}{T}^{\tilde{j}\tilde{k}},\end{array}\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{S}_{\tilde{\phi }}=-\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}[\alpha \sqrt{\gamma }({T}_{\tilde{\phi }}^{\tilde{r}}+{N}^{r}{S}_{\tilde{\phi )}}].\end{eqnarray} \tag{ 70 }$$

Here we use N^θ = 0 and N^ϕ = 0 at the equatorial plane in the KS coordinates.

In the BL coordinates, we have N^r = 0, N^θ = 0, αN^ϕ = ω, and ${g}_{{ij}}={h}_{i}^{2}{\delta }_{{ij}}$. We use the ZAMO frame ${x}^{\hat{i}}$ easily and obtain the 3+1 formalism equation with the ZAMO frame along the equatorial plane as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{B}^{\hat{\phi }}=-\displaystyle \frac{1}{{h}_{r}{h}_{\theta }}\displaystyle \frac{\partial }{\partial r}\left[{h}_{\theta }(\alpha {E}_{\hat{\theta }}-{h}_{\phi }\omega {B}^{\hat{r}})\right],\end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray}\,\begin{array}{r}\begin{array}{rcl}\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}t}{S}_{\hat{r}} & = & -\displaystyle \frac{1}{{\alpha }^{2}\sqrt{\gamma }}\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}r}\left[\displaystyle \frac{{\alpha }^{2}\sqrt{\gamma }}{{h}_{r}}{T}_{\hat{r}}^{\hat{r}}\right]-\displaystyle \frac{1}{{h}_{r}}\displaystyle \frac{{\rm{\partial }}\alpha }{{\rm{\partial }}r}(\hat{u}-{T}_{\hat{r}}^{\hat{r}})\\ & & -\,\displaystyle \sum _{i\ne r}\displaystyle \frac{\alpha }{{h}_{r}{h}_{j}}\displaystyle \frac{{\rm{\partial }}{h}_{i}}{{\rm{\partial }}r}{T}^{\hat{j}\hat{j}}-\sigma {S}_{\hat{\phi }},\end{array}\end{array}\,\end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{S}_{\hat{\phi }}=-\displaystyle \frac{1}{{h}_{\phi }\sqrt{\gamma }}\displaystyle \frac{\partial }{\partial r}\left[\displaystyle \frac{\alpha \sqrt{\gamma }{h}_{\phi }}{{h}_{r}}{T}_{\hat{\phi }}^{\hat{r}}\right],\end{eqnarray} \tag{ 73 }$$

where σ = (h_ϕ/h_r)∂ (αN^ϕ)/∂r. For the numerical calculation of 1D FFMD, we use Equations (68)–(70) and (71)–(73) for the numerical calculations in the KS and BL coordinates, respectively.

The flux coordinates for the three types of magnetic surfaces along the equatorial plane as shown in Figure 3 are given as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\theta }_{0}(r) & = & \displaystyle \frac{\pi }{2},\\ b(r) & = & \left\{\begin{array}{cc}{\left[1+0.25{\left(r-{r}_{{\rm{H}}}\right)}^{2}\right]}^{m} & (r\gt {r}_{{\rm{H}}})\\ 1 & (r\leqslant {r}_{{\rm{H}}})\end{array}\right..\end{array}\end{eqnarray} \tag{ 74 }$$

The incurvature-flared, radial, and excurvature-flared magnetic surfaces are given by Equation (74) with m > 0, m = 0, and m < 0, respectively.

Zoom In Zoom Out Reset image size

We use the Lax–Wendroff scheme and the simplified total variation diminishing scheme (Davis 1984) for the 1D FFMD numerical calculations. The mesh number of the standard case is 1000 (and 10,000, if necessary). In the KS coordinates, we use the even interval mesh with the range r_min ≤ r ≤ r_max (r_min = 0.9r_H, r_max = (40 ∼ 50)r_H).

In the BL coordinates, we use the even interval mesh on the (pseudo-)tortoise coordinate for the radial coordinate r_* such as ${{dr}}_{* }\propto \tfrac{1}{r-{r}_{0}}{dr}$ (r₀ < r_min is a constant, and we typically set r₀ = 1.9999 for the a = 0.01M case). When r₀ = r_H, the coordinate r_* becomes the exact tortoise coordinate. We use the equations of the Faraday law and the conservation equation of energy flux ${S}_{\tilde{i}}$ (i = r, ϕ).

Here we show some test calculations of 1D FFMD simulations for variously shaped magnetic surfaces around the equatorial plane, as shown in Figure 3. In the case of the radial magnetic surface, we can use the BL or KS coordinates as the flux coordinates. The 1D FFMD calculations along the equatorial plane are listed in Table 1.

4.3.1. Test Calculations of the Blandford–Znajek Solution

First we calculate the Blandford–Znajek solution on the monopole radial magnetic field as a test calculation of our numerical code in the BL and KS coordinates. We use the Blandford–Znajek solutions for the initial condition of time-development numerical simulations to check the 1D FFMD code due to the solutions being stationary. Note that we assume that the spin parameter of the black hole a_* ≡ a/M is much smaller than unity (a_* = 0.01). In the BL coordinates, the solution is given by

$$\begin{eqnarray}\begin{array}{rcl}{B}^{\hat{r}} & = & \displaystyle \frac{{B}_{0}}{{h}_{\theta }R}\left[1-{a}_{* }^{2}f(r)\right]\approx \displaystyle \frac{{B}_{0}}{{h}_{\theta }R},\\ {B}^{\hat{\theta }} & = & 0,{B}^{\hat{\phi }}=-\displaystyle \frac{{{aB}}_{0}}{8\alpha {{RM}}^{2}},\end{array}\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}&&{E}_{\hat{\theta }}=-\displaystyle \frac{R}{\alpha }({{\rm{\Omega }}}_{{\rm{F}}}-\omega ){B}^{\hat{r}},\ {E}_{\hat{r}}={E}_{\hat{\phi }}=0,\end{eqnarray} \tag{ 76 }$$

where Ω_F = a/(8 M²) and f(r) is a complex function, which is given in Section 6 of Blandford & Znajek (1977). Here we use the flux coordinates (t, r, Ψ, ϕ) = (t, r, θ, ϕ); thus, we must set B₀ = 1. In the Blandford–Znajek solution, the slowly spinning black hole limit condition (a_* ≪ 1) is used, and we neglect the term of f(r).

In the KS coordinates, the Blandford–Znajek solution is given by

$$\begin{eqnarray}\begin{array}{rcl}{B}^{\tilde{r}} & = & \displaystyle \frac{{B}_{0}}{\sqrt{\gamma }},\ {B}^{\tilde{\theta }}=0,\\ {B}^{\tilde{\phi }} & = & -\displaystyle \frac{a\alpha {B}_{0}}{2{{Mr}}^{3}}\left(1+\displaystyle \frac{r}{4M}\right),\end{array}\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}&&{E}_{\tilde{\theta }}=-\displaystyle \frac{{{aB}}_{0}}{8{M}^{2}\alpha }-\sqrt{\gamma }{N}^{r}{B}^{\tilde{\phi }},{E}_{\tilde{r}}={E}_{\hat{\phi }}=0.\end{eqnarray} \tag{ 78 }$$

In both the BL and KS coordinates, the numerical results demonstrate that the solutions are stationary, as shown in Figures 4 and 5, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3.2. Test Calculation of the Electromagnetic Field in a Vacuum

In the force-free condition, we have the solution of the electromagnetic field in a vacuum. It is obtained from the Maxwell equations with ρ_e = 0 and ${\boldsymbol{J}}=0$. In the BL coordinates, we have

$$\begin{eqnarray}&&{B}^{\hat{r}}=\displaystyle \frac{{B}_{0}}{{h}_{\theta }R},\ {B}^{\hat{\theta }}=0,{B}^{\hat{\phi }}=0,\end{eqnarray} \tag{ 79 }$$

$$\begin{eqnarray}&&{E}_{\hat{\theta }}={h}_{\theta }{N}^{\phi }{B}^{\hat{r}}=\displaystyle \frac{R\omega }{\alpha }{B}^{\hat{r}},{E}_{\hat{r}}={E}_{\hat{\phi }}=0.\end{eqnarray} \tag{ 80 }$$

In the KS coordinates, the solution of the electromagnetic field in a vacuum is given by

$$\begin{eqnarray}&&{B}^{\tilde{r}}=\displaystyle \frac{{B}_{0}}{\sqrt{\gamma }},\ {B}^{\tilde{\theta }}=0,{B}^{\tilde{\phi }}=-\displaystyle \frac{{{aB}}_{\tilde{r}}}{{\rm{\Delta }}(1-2M/r){r}^{2}},\end{eqnarray} \tag{ 81 }$$

$$\begin{eqnarray}&&{E}_{\tilde{\theta }}=-\sqrt{\gamma }{B}^{\tilde{r}}{B}^{\tilde{\phi }},{E}_{\tilde{r}}={E}_{\hat{\phi }}=0.\end{eqnarray} \tag{ 82 }$$

Note that ${B}^{\tilde{\phi }}$ diverges at the horizon in the KS coordinates, whose metrics are all finite. This indicates that divergence of the solution is physical, and we cannot use this solution.

We perform a test calculation with the vacuum solution in the BL coordinates, as shown in Figure 6. The result shows that I, P, and Ω_F are negligibly small, and, as expected, the solution is nearly stationary. The numerical calculations in Sections 4.3.1 and 4.3.2 indicate the reliability of the 1D FFMD code.

Zoom In Zoom Out Reset image size

4.3.3. Test Calculations of the Dynamic Process toward the Blandford–Znajek Solution with Small a_*

Here we perform numerical test calculations for the emergence of the Blandford–Znajek mechanism around a slowly spinning black hole of a_* = 0.01 with both BL and KS coordinates. Note that the initial condition I = 0 (P = 0) is used in the cases of both the BL and KS coordinates. First, we show the result of the BL coordinates. When I = 0 (P = 0) in the BL coordinates of a_* = 0.01, we have ${B}^{\hat{\phi }}=0$. Figure 7 shows that no quantity ever changes at the horizon due to the time-freezing at the horizon. A wave is caused around the ergosphere, and the front propagates outward like a "tsunami" and inward to the horizon. The constants of the stationary state in the wave are realized as the Blandford–Znajek solution with Ω_F = ω_H/2 = a_*/8M = 0.00125. The energy flux is provided from the front of the inward wave near the horizon. In the front of the inward wave, the energy at infinity becomes negative. To observe this effect, we check the above statement near the black hole, as shown in Figure 8. We find that the profile of P at t = 4M has its maximum value at r = 2.2M, and its value at the horizon is zero. Comparing the profile of P at t = 4M and 8M indicates that the profile spreads both outward and inward. The outward spread is shown in Figure 7. Note that the inward spread continues to approach to the horizon but never actually reaches it. Then the energy flux emerges from the ergosphere. As shown in Figure 7, the energy at infinity ${e}^{\infty }$ becomes negative at t = 40M.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the KS coordinates, we have the initial condition I = 0 (P = 0) and ${E}_{\tilde{\theta }}=0$. The results are shown in Figure 9. The quantities at the horizon rapidly converge to the values of the Blandford–Znajek solution, and the region similar to the Blandford–Znajek solution spreads outward. It appears that the horizon behaves like a boundary, while "causality" prohibits outward propagation of the information through the horizon. This result also demonstrates that the outward energy flux at the horizon increases spontaneously to reach the value of the Blandford–Znajek solution even in the case of initially no current and energy flux.

Zoom In Zoom Out Reset image size

4.3.4. Test Calculations of the Blandford–Znajek Mechanism with Finite a_*

We perform 1D FFMD numerical simulations with the steady-state solution for finite a_* at the equatorial plane in the flux coordinates based on the BL and KS coordinates. In the force-free steady state, we have constants I, P, and Ω_F, and the azimuthal magnetic field component and colatitude electric field component are given by Equations (41) and (42), respectively.

When the magnetic surface locates at the equatorial plane (θ₀ (r) = π/2 and K = 0), using Equations (58) and (59) with ${\theta }_{{\rm{H}}}={\theta }_{\infty }=\pi /2$, we obtain the following constants:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{H}}}=\displaystyle \frac{a}{(({b}_{\infty }/{b}_{{\rm{H}}}){r}_{{\rm{H}}}+2M){r}_{{\rm{H}}}},\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&I=-\displaystyle \frac{{{ab}}_{\infty }}{(({b}_{\infty }/{b}_{{\rm{H}}}){r}_{{\rm{H}}}+2M){r}_{{\rm{H}}}}.\end{eqnarray} \tag{ 84 }$$

Figures 10 and 11 show the simulation results with the steady-state solution (Equations (83) and (84)) for the case of a radial magnetic surface (${b}_{\infty }={b}_{{\rm{H}}}$), where the flux coordinates are the BL and KS coordinates, respectively, for a_* = 0.95. Here the horizon locates at r = 1.31M. This clearly confirms that the solution gives the steady state of the electromagnetic field with the outward power predicted by $P={{\rm{\Omega }}}_{{\rm{F}}}^{2}=0.048$, where we use the unit system such that b_H = 1. Hereafter, we use the unit system with b_H = 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3.5. Test Calculations of Pulse Propagation on the Blandford–Znajek Solution with Finite a_*

To investigate causality, we follow the propagation of the pulse initially added to the steady-state field around a rapidly spinning black hole (a_* = 0.95). Here we initially add the perturbation of the pulse to the azimuthal component of the magnetic field. The additional azimuthal component of the magnetic field of the pulse to the stationary fields at r = r₀ with width δ and amplitude A is given by

$$\begin{eqnarray}&&\begin{array}{l}{B}_{\mathrm{pul}}^{\tilde{\phi }}\\ \ =\ \left\{\begin{array}{cc}\ \displaystyle \frac{A}{2}\left(1+\cos \pi \displaystyle \frac{r-{r}_{0}}{\delta /2}\right)\ & \left({r}_{0}-\displaystyle \frac{\delta }{2}\leqslant r\leqslant {r}_{0}+\displaystyle \frac{\delta }{2}\right)\\ 0\ & \left(r\lt {r}_{0}-\displaystyle \frac{\delta }{2},r\gt {r}_{0}+\displaystyle \frac{\delta }{2}\right).\end{array}\right.\end{array}\end{eqnarray} \tag{ 85 }$$

Figure 12 shows the propagation of the pulse initially at r = 10M on the background steady-state electromagnetic field as obtained by Equations (83) and (84) for a_* = 0.95. Two pulses propagate outward and inward with velocity 0.7 and 1.2, respectively. This demonstrates that the information can propagate both inward and outward at nearly the speed of light outside of the ergosphere.

Zoom In Zoom Out Reset image size

Figure 13 shows the propagation of the pulse initially located at r = 1.8M inside the ergosphere, where the background is the same as that in Figure 11. As shown, two pulses also propagate outward and inward; however, the inward pulse runs fast, and the outward pulse propagates very slowly.

Zoom In Zoom Out Reset image size

Figure 14 shows the propagation of the pulse initially located at r = 0.95M inside the horizon, where the background is the same as that in Figure 11. Here we find two pulses at t = 0.1M (dotted line). One pulse at r = 0.82M propagates inward rapidly. At t = 0.2M, this pulse passes through the horizon. The pulse at r = 0.945 propagates inward very slowly and appears to be nearly at rest. This confirms that the energy can be transported outward even inside of the horizon, but information is never transported outward. This result is consistent with the notion of causality around a black hole.

Zoom In Zoom Out Reset image size

4.3.6. Test Calculations of the Emerging Blandford–Znajek Process with Finite a_*

To investigate the emergence process of the Blandford–Znajek mechanism, we run several simulations with the KS coordinates and observe the time-dependent process of the force-free field inside, at, and outside the horizon. To find the details of the process between the static limit and the horizon, we put the spin parameter of the black hole as a_* = 0.95. We also perform simulations with the BL coordinates to investigate the dependence on the coordinate systems. Here we use the initial condition with ${B}^{\tilde{\phi }}=0$ and ${E}_{\tilde{\theta }}=0$ around the spinning black hole.

Figure 15 shows the time development of the power P, current I, and angular velocity of the magnetic field lines Ω_F on the BL coordinates. The power and angular velocity of the magnetic field lines are zero at the initial condition. The results obtained at t = 20M (dotted line) and t = 40M (solid line) show that finite regions of P and Ω_F spread outward at the speed of light. Note that P and Ω_F converge to the values expected by the steady state, i.e., P = Ω_F² = 0.048 with Equations (83) and (84). Then it seems that the steady-state outward energy flux is supplied in the vicinity of the black hole. Figure 16 shows an enlargement of Figure 15 near the horizon. Here we observe inward propagation of the tsunami of P and Ω_F in the BL coordinates, which is caused by a time lapse of the coordinates near the horizon.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 17 shows the time development of P, I, and Ω_F on the KS coordinates. The power and angular velocity of the magnetic field lines are zero in the initial condition. The results at t = 30M (dotted line) and t = 60M (solid line) show that finite regions of P and Ω_F spread outward at the speed of light. Note that P and Ω_F converge to the values expected by the steady state, i.e., P = Ω_F² = 0.048 with Equations (84) and (83). Then it seems that the steady-state outward energy flux comes from the vicinity of the black hole.

Zoom In Zoom Out Reset image size

4.3.7. Test Calculations of the Emergence of the Blandford–Znajek Mechanism along the Incurvature-flared Magnetic Surface at the Equatorial Plane

Here we show a numerical calculation with an incurvature-flared magnetic surface given by Equation (74) with m = 0.25 imaged in Figure 3(c).

Figure 18 shows the steady-state solution of P, I, and Ω_F along an incurvature-flared magnetic surface at the equatorial plane around a spinning black hole with a_* = 0.95.

Zoom In Zoom Out Reset image size

Figure 19 shows the time development of the P, I, and Ω_F of the force-free electromagnetic field along an incurvature-flared magnetic surface at the equatorial plane around a spinning black hole with a_* = 0.95. Here the horizontal solid red lines show the steady-state solution. In the incurvature-flared magnetic surface case, the electromagnetic field tends to converge very gradually to the steady state (details are provided by Imamura & Koide 2019).

Zoom In Zoom Out Reset image size

4.3.8. Test Calculations of the Emergence of the Blandford–Znajek Mechanism along the Excurvature-flared Magnetic Surface at the Equatorial Plane

Here we show a numerical calculation with an excurvature-flared magnetic surface given by Equation (74) with m = −0.25 imaged in Figure 3(b). Figure 20 shows the time development of the P, I, and Ω_F of the force-free electromagnetic field along an excurvature-flared magnetic surface at the equatorial plane around a spinning black hole with a_* = 0.95. The horizontal solid red lines show the steady-state solution. In the excurvature-flared magnetic surface case, the electromagnetic field does not appear to converge to the steady state (details are given by Imamura & Koide 2019).

Zoom In Zoom Out Reset image size