Near-horizon Structure of Escape Zones of Electrically Charged Particles around Weakly Magnetized Rotating Black Hole

Since the seminal works of Blandford & Znajek (1977) and Blandford & Payne (1982), rotation of the black hole with the surrounding accretion disk and a large-scale magnetic field are considered key ingredients for the effective acceleration and collimation of outflows of matter from the vicinity of accreting black holes. The Blandford–Znajek process operates near the event horizon owing to the frame-dragging effect, which amplifies and twists the magnetic field anchored in the accretion disk. Similarly, the Blandford–Payne process requires a magnetized accretion disk. As a result, the charged matter of the disk may be accelerated outward in jets aligned with the spin axis, transporting away the rotational energy of the black hole. On the other hand, small-scale magnetic fields may be induced in the accretion disk owing to turbulence driven by the magnetorotational instability (Balbus & Hawley 1998). While turbulence is effective in the transport of the angular momentum (i.e., it contributes to the viscosity of the disk), the resulting small-scale field structure seems insufficient for launching the jet (Beckwith et al. 2008; Gold et al. 2017).

The crucial role of the ordered large-scale magnetic fields in powering the outflows (and astrophysical jets, in particular) is being confirmed by numerical experiments. Indeed, in 3D general relativistic megnetohydrodynamic (GRMHD) simulations of accretion processes, ordered bundles of magnetic field lines seem to emerge in a self-consistent manner under rather general conditions (e.g., Penna et al. 2010; Tchekhovskoy 2015; Sadowski 2016). The concept of the Blandford–Znajek or Blandford–Payne process is thus further supported as a plausible mechanism of launching the outflow.

In the given context we investigate the motion of charged particles (e.g., heavy ions, dust grains) near the magnetized rotating black hole. In particular, we discuss under which circumstances the material bound to stable orbits embedded in the equatorial plane may be set on the escaping trajectories, liberated from the attraction of the center and even accelerated to relativistic velocities. We consider a single particle limit, i.e., the particles are considered as noninteracting and hydrodynamical terms are neglected in governing equations. Direct applicability of the model to astrophysical systems is thus limited to the regions where the gas is sufficiently diluted (the mean free path of particles is larger than the characteristic length scale of the system given by the radius of event horizon). Nevertheless, the main motivation of this study lies in conceptual questions concerning the effect of the magnetic field on the stability of orbits near the strongly gravitating center.

Previously, we have demonstrated that the structure of the electromagnetic field is strongly influenced by the relativistic effects of a spinning compact object (Karas & Kopáček 2009; Karas et al. 2014, 2012). In particular, we have shown that close to the event horizon, where the frame dragging dominates, these effects shape the field lines (described by the vacuum solution of Maxwell's equations) in a way that may mimic the behavior typically associated with magnetohydrodynamic plasma. Namely, we observed that the structure of the field in some regions of spacetime may correspond to current sheets or X-type null points of magnetic reconnection. As the field's topology changes profoundly owing to gravity, also the motion of charged matter is strongly affected. As a most striking consequence, the particle dynamics undergoes a transition from regularity to deterministic chaos even in the axisymmetric case of uniform magnetic field aligned along the rotation axis (Kopáček et al. 2010a, 2010b). Moreover, we have seen (Kopáček & Karas 2014) that any deviation from the perfect axisymmetry leads to overall dominance of chaotic dynamics in the given setup.

According to general relativity, the energy density of the electromagnetic field (as well as any other field) stands on the right-hand side in Einstein's field equations for the spacetime structure. Hence, the electromagnetic field contributes to the metric; to derive a self-consistent solution, the coupled Einstein and Maxwell field equations must be solved simultaneously. However, the energy density in astrophysically realistic electromagnetic fields (proportional to the square of intensities, $\propto {E}^{2}+{B}^{2}$) is far too low to influence the spacetime noticeably. In these circumstances, test-field (linearized in E and B terms) solutions are adequate and accurate enough for describing the motion of plasma and particles, while the corresponding exact solutions of the Einstein–Maxwell set of equations are mainly of theoretical rather than astrophysical interest.

In particular, in the present discussion we adopt the axisymmetric test-field solution corresponding to an asymptotically uniform magnetic field aligned with the rotation axis (Wald 1974). This assumption seems astrophysically relevant, and also GRMHD simulations reveal a "magneto-spin alignment" mechanism causing magnetized disks and jets to align with the spin near black holes (McKinney et al. 2013). We have seen (Kopáček et al. 2010a) that in this setup regular and chaotic trajectories of charged particles generally coexist since the magnetic field acts as a nonintegrable perturbation. Besides triggering chaos, the field also changes the topology of the effective potential in such a way that it in principle allows the particle to leave the equatorial plane and escape to the asymptotic region in the narrow corridor aligned with the symmetry axis (see Figures 12–14 in Kopáček et al. 2010a). In this paper we shall investigate this phenomenon systematically and also make some astrophysically motivated assumptions considering the initial conditions of the trajectories. Namely, we focus on initially neutral particles on equatorial co-rotating circular (Keplerian) orbits. Below the innermost stable circular orbit (ISCO), the particle is supposed to freely fall onto the horizon. At some radius r₀ (below or above the ISCO), the particle obtains electric charge (by photoionization or other process), and its dynamics changes owing to magnetic field. Within this scenario we study under which circumstances the escaping trajectories may be realized and discuss the maximal acceleration that can be achieved.

Charged particles escaping the attraction of black holes immersed in magnetic field have already been the subject of several studies. Weakly magnetized Schwarzschild black holes have been considered by Frolov & Shoom (2010) and Al Zahrani et al. (2013). An exact solution describing the static black hole immersed in asymptotically uniform magnetic field (Ernst's spacetime) was employed by Huang et al. (2015). The case of a spinning black hole with a magnetic test field has also been analyzed (Al Zahrani 2014; Shiose et al. 2014). The special case of a slowly rotating black hole was considered by Hussain et al. (2014). Moreover, trajectories of particles escaping from the vicinity of magnetized naked singularities were also examined (Babar et al. 2016). In these studies, the initial setup consists of circular equatorial orbits of charged particles that are kicked from the equatorial plane mechanically by the collision with another particle or photon. The immediate cause of destabilization of a circular orbit that may lead to escaping trajectory is thus an ad hoc introduction of velocity components perpendicular to the equatorial plane. In our scenario, however, no such kick is needed, as we destabilize the equatorial geodesic orbit of initially neutral particles by the charging process itself. Analogical setup was recently considered by Stuchlík & Kološ (2016) with the class of spherical orbits (Wilkins 1972), and the authors demonstrated that charging the neutral particle on such an orbit (with a nonzero velocity component in the direction normal to the equatorial plane) may initiate its escape. Here we discuss in detail the effect of ionization of particles on circular (or inspiraling below ISCO) equatorial orbits, which have no obvious tendency to leave the equatorial plane, and thus their escape seems less probable compared to the above-referenced models.

The paper is organized as follows. In Section 2 we describe the employed model of the weakly magnetized rotating black hole and review the equations of motion for charged particles. Trajectories of particles escaping from the circular orbits in the equatorial plane are discussed in Section 3. Using the effective potential method, we first identify the necessary conditions for the escape analytically in Section 3.1. In the following we numerically integrate the trajectories and discuss the role of the black hole's spin and magnetic field in the emergence of escaping orbits (Section 3.2). Acceleration of escaping particles is investigated quantitatively in Section 3.3. Results are briefly discussed in Section 4 and summarized in Section 5.

Kerr metric describing the geometry of the spacetime around the rotating black hole may be expressed in Boyer–Lindquist coordinates ${x}^{\mu }=(t,r,\theta ,\varphi )$ as follows (Misner et al. 1973):

$$\begin{eqnarray}&&{{ds}}^{2}=-\displaystyle \frac{{\rm{\Delta }}}{{\rm{\Sigma }}}\ {[{dt}-a\sin \theta d\varphi ]}^{2}+\displaystyle \frac{{\sin }^{2}\theta }{{\rm{\Sigma }}}\\ &&\quad \times \,{[({r}^{2}+{a}^{2})d\varphi -a{dt}]}^{2}+\displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}}\,{{dr}}^{2}+{\rm{\Sigma }}d{\theta }^{2},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{\rm{\Delta }}\equiv {r}^{2}-2{Mr}+{a}^{2},\,{\rm{\Sigma }}\equiv {r}^{2}+{a}^{2}{\cos }^{2}\theta .\end{eqnarray} \tag{ 2 }$$

The coordinate singularity at ${\rm{\Delta }}=0$ corresponds to the outer/inner horizon of the black hole, ${r}_{\pm }=M\pm \sqrt{{M}^{2}-{a}^{2}}$. Rotation of the black hole is measured by the spin parameter $a\in \langle -M,M\rangle$. Here we only consider $a\geqslant 0$ without the loss of generality.

We stress that geometrized units are used throughout the paper. Values of basic constants (gravitational constant G, speed of light c, Boltzmann constant k, and Coulomb constant k_C) therefore equal unity, $G=c=k={k}_{C}=1$.

We consider the presence of an organized magnetic field threading the black hole horizon. Black holes do not support their own magnetic field of intrinsic origin, like, e.g., stars can. In order to endow a rotating black hole with a weak, large-scale magnetic component, we employ the test-field solution (Wald 1974) for an asymptotically uniform magnetic field aligned with the rotation axis (i.e., the system remains axially symmetric and stationary). Near the black hole, the field becomes twisted owing to the frame-dragging effect. Moreover, an electric field component is induced here. The vector potential can be expressed in terms of Kerr metric coefficients (1) as follows:

$$\begin{eqnarray}{A}_{t} & = & \tfrac{1}{2}B({g}_{t\phi }+2a\,{g}_{{tt}})\\ {A}_{\phi } & = & \tfrac{1}{2}B({g}_{\phi \phi }+2{{ag}}_{t\phi }),\end{eqnarray} \tag{ 3 }$$

where B is magnetic intensity in the asymptotic region.

Equations of motion for charged particles of charge q and rest mass m may be derived using Hamiltonian ${ \mathcal H }$, defined as

$$\begin{eqnarray}&&{ \mathcal H }=\tfrac{1}{2}{g}^{\mu \nu }({\pi }_{\mu }-{{qA}}_{\mu })({\pi }_{\nu }-{{qA}}_{\nu }),\end{eqnarray} \tag{ 4 }$$

where ${\pi }_{\mu }$ is the generalized (canonical) momentum and ${g}^{\mu \nu }$ denotes contravariant components of the metric tensor. The equations of motion are given as

$$\begin{eqnarray}&&\displaystyle \frac{{{dx}}^{\mu }}{d\lambda }\equiv {p}^{\mu }=\displaystyle \frac{\partial { \mathcal H }}{\partial {\pi }_{\mu }},\quad \displaystyle \frac{d{\pi }_{\mu }}{d\lambda }=-\displaystyle \frac{\partial { \mathcal H }}{\partial {x}^{\mu }},\end{eqnarray} \tag{ 5 }$$

where $\lambda \equiv \tau /m$ is the affine parameter, which becomes dimensionless in geometrized units (τ denotes the proper time). Kinematical four-momentum ${p}^{\mu }$ is given as ${p}^{\mu }={\pi }^{\mu }-{{qA}}^{\mu }$, which follows directly from the first equation of Equation (5). The conserved value of the Hamiltonian is therefore given as ${ \mathcal H }=-{m}^{2}/2$. Moreover, the system is stationary, and momentum ${\pi }_{t}$ is thus an integral of motion expressing (negatively taken) energy of the test particle ${\pi }_{t}\equiv -E$. Axial symmetry of the system ensures that the corresponding component of the angular momentum also remains conserved as another integral of motion, ${\pi }_{\varphi }\equiv L$. Nevertheless, in the following we switch to specific quantities $E/m\to E$, $L/m\to L$, and $q/m\to q$ when describing the test particle.

The effective potential expressing minimal allowed energy of charged test particles in axisymmetric stationary spacetime may be derived as follows (Kopáček et al. 2010a):

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}(r,\theta )=(-\beta +\sqrt{{\beta }^{2}-4\alpha \gamma })/2\alpha ,\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&\alpha =-{g}^{{tt}},\quad \beta =2[{g}^{t\varphi }(L-{{qA}}_{\varphi })-{g}^{{tt}}{{qA}}_{t}],\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\gamma =-{g}^{\varphi \varphi }\\ &&\,\times \,{(L-{{qA}}_{\varphi })}^{2}-{g}^{{tt}}{q}^{2}{A}_{t}^{2}+2{g}^{t\varphi }{{qA}}_{t}(L-{{qA}}_{\varphi })-{m}^{2}.\end{eqnarray} \tag{ 8 }$$

Analysis of the two-dimensional effective potential (6) may be used to localize circular orbits not only in the equatorial plane (Kološ et al. 2015) but also off this plane (Kovář et al. 2008, 2010, 2013), and it allows us to discuss the stability of resulting orbits (Tursunov et al. 2016).

3.1. Charging the Matter and Necessary Conditions for the Escape

We consider a standard (thin-disk) accretion geometry (Novikov & Thorne 1973) in which the electrically neutral matter gradually sinks between co-rotating (prograde, equatorial) Keplerian orbits whose energy ${E}_{\mathrm{Kep}}(r)$ and angular momentum ${L}_{\mathrm{Kep}}(r)$ are given as follows (Bardeen et al. 1972):

$$\begin{eqnarray}{E}_{\mathrm{Kep}} & = & \displaystyle \frac{{r}^{2}-2{Mr}+a\sqrt{{Mr}}}{r\sqrt{{r}^{2}-3{Mr}+2a\sqrt{{Mr}}}},\\ {L}_{\mathrm{Kep}} & = & \displaystyle \frac{\sqrt{M}({r}^{2}+{a}^{2}-2a\sqrt{{Mr}})}{\sqrt{r({r}^{2}-3{Mr}+2a\sqrt{{Mr}})}}.\end{eqnarray} \tag{ 9 }$$

Co-rotating circular geodesics are not allowed below the radius of marginally stable orbit ${r}_{\mathrm{ms}}$, whose position is given as

$$\begin{eqnarray}&&{r}_{\mathrm{ms}}=M(3+{Z}_{2}-\sqrt{(3-{Z}_{1})(3+{Z}_{1}+2{Z}_{2})}),\end{eqnarray} \tag{ 10 }$$

where ${Z}_{1}\equiv 1+{\left(1-\tfrac{{a}^{2}}{{M}^{2}}\right)}^{1/3}\left[{\left(1+\tfrac{a}{M}\right)}^{1/3}+{\left(1-\tfrac{a}{M}\right)}^{1/3}\right]$ and ${Z}_{2}\equiv \sqrt{3{a}^{2}/{M}^{2}+{Z}_{1}^{2}}$.

The ISCO is located at ${r}_{\mathrm{ms}}$. We suppose that below ISCO the geodesics are freely falling onto the horizon of the black hole with the energy and angular momentum corresponding to the marginally stable orbit, i.e., the particles keep $E={E}_{\mathrm{Kep}}({r}_{\mathrm{ms}})$ and $L={L}_{\mathrm{Kep}}({r}_{\mathrm{ms}})$ during their inspiral.

However, neutral elements may undergo a sudden charging process (due to photoionization) at a given radius and obtain nonzero specific charge q. While the change of the rest mass m may be neglected, the particle dynamics changes immediately owing to the presence of the magnetic field (3). Values of integrals of motion are changed as follows:

$$\begin{eqnarray}&&E={E}_{\mathrm{Kep}}-{{qA}}_{t},\quad L={L}_{\mathrm{Kep}}+{{qA}}_{\varphi }.\end{eqnarray} \tag{ 11 }$$

The effective potential of the particle changes accordingly. From Equation (6) we obtain

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}(q)={V}_{\mathrm{eff}}(q=0)-{{qA}}_{t},\end{eqnarray} \tag{ 12 }$$

which ensures that although the charge shifts the energy of the particle, the difference separating the particle from the allowed minimum (defined by the effective potential) remains unchanged. In particular, stable Keplerian orbits found in local minima of ${V}_{\mathrm{eff}}$ remain on the potential curve even after the ionization. However, they are not located in the potential minima anymore; the former circular orbit is disturbed. As a result, it turns into a stable, radially oscillating (elliptic) orbit, or it may become an unstable (plunging) orbit inspiraling onto the horizon of the black hole. Interestingly, it can also be set on an escaping trajectory. To inspect the latter option, we examine the asymptotic behavior of the effective potential Equation (6). For large r we find the dominant term to be given as

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}\propto \displaystyle \frac{| {qB}| r\sin \theta }{2},\end{eqnarray} \tag{ 13 }$$

which shows that particles may only escape to infinity along the axis of symmetry, in which case $r\sin \theta$ remains finite.

In particular, for a particle charged on a Keplerian orbit of radius r₀ with former energy ${E}_{\mathrm{Kep}}$ escaping along the symmetry axis, we obtain the following relation valid in the asymptotic region:

$$\begin{eqnarray}&&E-{V}_{\mathrm{eff}}{| }_{r\gg M}={E}_{\mathrm{Kep}}-\displaystyle \frac{{aqB}}{{r}_{0}}-\sqrt{1+\displaystyle \frac{{q}^{2}{B}^{2}{r}^{2}\sin {\theta }^{2}}{4}}.\end{eqnarray} \tag{ 14 }$$

The motion is allowed only for non-negative values of $E-{V}_{\mathrm{eff}}$, and since ${E}_{\mathrm{Kep}}\lt 1$ for finite r₀ and $a\geqslant 0$ is supposed, we observe that (i) particles may only escape for ${qB}\lt 0$, (ii) escape is possible only for $a\ne 0$, and (iii) the asymptotic velocity of escaping particles is an increasing function of parameters $| {qB}|$ and a and a decreasing function of r₀.

The first condition, ${qB}\lt 0$, requires that the charge of the particle is negative for the parallel orientation of the spin and magnetic field and positive for the antiparallel orientation. The escaping orbits are not allowed in the Schwarzschild limit, a = 0, regardless of the values of q, B, and r₀. The asymptotic velocity of escaping particles rises with increasing values of a, $| q|$, and $| B|$ and will be higher for the particles charged closer to the horizon.

Thus, we conclude that escaping orbits are in principle possible, i.e., that necessary conditions for the escape may be fulfilled in a given setup. However, these conditions are not sufficient, and they give no guarantee that trajectories escaping from the equatorial plane are actually realized. To test whether ionized particles really escape, and for the discussion of their asymptotic velocity, we need to switch to the numerical approach and integrate particular trajectories.

3.2. Formation of the Escape Zone

The crucial parameters that determine the behavior of the Keplerian particle after it obtains electric charge are the black hole spin a, the ionization radius r₀, and the product qB (hereafter, magnetization parameter).¹ In the rest of the paper we scale all quantities by the black hole's rest mass M, which corresponds to setting M = 1 in the equations and leaves all the quantities dimensionless.

In the search of escaping particles we need to cover an adequate range of values of crucial parameters. Therefore, we set up a grid where the vertical axis will correspond to the spin in the range $a\in \langle 0,1\rangle$ and the horizontal axis will represent the ionization radius ${r}_{0}\in \langle 1,{r}_{\mathrm{out}}\rangle$, where ${r}_{\mathrm{out}}$ is chosen such that we do not omit any potentially interesting part of the accretion disk. The number of data points of the grid would be typically 200 × 200, and we integrate the trajectory corresponding to each point of the grid with a fixed value of qB = const. Each integration is set to run up to the value of the affine parameter $\lambda =2\times {10}^{3}$, which corresponds to $\approx 100$ azimuthal periods of the original Keplerian orbit near the ISCO. The threshold above which the particle is considered escaping is set to $r={10}^{2}$. On the other hand, the orbit is regarded as plunging if it closely approaches the horizon $r={r}_{+}+{10}^{-4}$. If the final radius lies between both limits, the respective particle is considered stable.

Results of the integration are visualized by assigning different colors to grid points corresponding to trajectories with different final states. This type of graphic representation of the dynamic system is usually denoted as a basin-boundary plot (see, e.g., Frolov & Shoom 2010; Al Zahrani et al. 2013; Al Zahrani 2014, where basin-boundary plots are applied in a similar context). We use blue color for plunging orbits, red for stable orbits, and yellow for escaping ones. We begin with the geodesic limit qB = 0, for which no escapes are allowed and the borderline between the blue plunging region and the red area of stable orbits is given by the ISCO curve. The analysis in Section 3.1 shows that only particles with ${qB}\lt 0$ may acquire sufficient energy to reach infinity. Thus, we introduce the charge and set ${qB}=-0.1$. The borderline between the regions is no longer defined by the ISCO. In general, the charge may cause the originally stable orbit to become plunging and vice versa, as some plunging neutral particles become stabilized by the ionization during their infall below the ISCO. Nevertheless, no escapes are observed for ${qB}=-0.1$.

We gradually decrease the value of qB and eventually detect escaping orbits for ${qB}=-0.5$ (top left panel of Figure 1). For this case, the escaping orbits are only allowed for spin values close to the extremal value a = 1, and they spread rather sporadically near the borderline separating the blue infall zone from the red area of stable orbits. In the top right panel of Figure 1 we present the plot for ${qB}=-1$, for which the density of escaping trajectories substantially increases and they form a narrow escape zone that extends down to spin values as low as $a\approx 0.3$. This plot also illustrates that even freely falling particles below the ISCO may still escape the attraction of the center owing to the charging process. Switching to ${qB}=-5$ (bottom left panel of Figure 1) results in shifting the escape zone further to lower spins and simultaneously to lower radial values, i.e., closer to the horizon. We also observe that for this value of qB the upper end of the escape zone detaches from the a = 1 limit and escape of particles from the vicinity of the black hole with spin close to extremal is not permitted anymore. Increasing the magnitude of the magnetization parameter to ${qB}=-20$ (bottom right panel of Figure 1) confirms both these trends and the escape zone shifts further down to lower spin (though it can never reach the static limit a = 0) and becomes generally closer to the horizon.

Zoom In Zoom Out Reset image size

Properties of the escape zone are further explored in Figures 2 and 3. In the former we observe in detail that the escape zone detaches from the extremal spin a = 1 as qB decreases. Indeed, while for ${qB}=-4.5$ there are still some escaping particles with a = 1 (left panel of Figure 2), none of them remain if qB decreases to ${qB}=-4.6$ (right panel of Figure 2). We conclude that for ${qB}\lessapprox -4.5$ the escaping trajectories are only realized for $a\lt 1$. In Figure 3 we show the effect of the magnetization parameter and spin on escaping orbits from yet another perspective. While for small values of $| {qB}|$ the escape is only possible for spin close to the extremal value (left panel of Figure 3 with ${qB}=-0.5$), for higher $| {qB}|$ the escape zone is shifted toward lower values of spin, escape is allowed even for very low a, and high spins are actually excluded (right panel of Figure 3 with ${qB}=-20$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to illustrate the dynamics of stable and escaping orbits in the escape zone, we plot the poloidal projection of typical trajectories along with the corresponding isopotential contours in Figure 4. In accordance with the results of the analysis in Section 3.1, we observe that the isopotentials for the negatively charged particles do not close and create the infinite narrow corridor parallel to the rotation axis, making the escape in principle possible. However, only some particles do actually escape (such as the one in the middle panel of Figure 4). In the left and right panels we observe that the particle may also remain oscillating on the stable orbit in or around the equatorial plane, depending on the initial radius r₀. To inspect the dynamics of stable orbits, we can employ standard techniques—e.g., Lyapunov characteristic exponents (Skokos 2010) or Poincaré sections. The latter is applied in Figure 5 (middle panel) for the case of orbits with ${r}_{0}=4.82,a=1$, and ${qB}=-1$ (the same orbit as in the left panel of Figure 4). The trajectory is integrated up to ${\lambda }_{\mathrm{fin}}={10}^{5}$, resulting in $\approx {10}^{4}$ section points that form a closed curve typical for a regular quasi-periodic orbit. Nevertheless, in the case of escaping trajectory these standard tools are not applicable, as the particle escapes too early, leaving only a few section points in the equatorial plane, which do not allow us to decide whether they form a curve or rather a surface. In addition, the computation of the Lyapunov characteristic exponents requires a much longer integration period to converge (denoted as Lyapunov time; Skokos 2010). However, we can still apply recurrence plots (Marwan et al. 2007), which are known to reflect the characteristic dynamic properties on a much shorter timescale. Recurrence plots (RPs) represent the recurrences of the trajectory to the defined ε neighborhood in the phase space. Visual survey of the RPs allows us to discriminate between the dynamic regimes of the system. In particular, regular orbits generally form diagonal patterns in the RP as shown for the stable orbit in the right panel of Figure 5. On the other hand, chaotic dynamics is characterized by disrupted diagonal patterns as observed for the escaping trajectory in the right panel of Figure 6. The given particle is integrated up to ${\lambda }_{\mathrm{fin}}=575$ (when its escape begins), which is sufficient to produce an indicative RP. For the details on recurrence analysis we refer to the review by Marwan et al. (2007), while the full description of the application of RPs in the analogous context may be found in our previous paper (Kopáček et al. 2010a). For a recent application of this method in relativistic astrophysics, see Semerák & Suková (2012), Witzany et al. (2015), Suková & Janiuk (2016), and Lukes-Gerakopoulos & Kopáček (2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the above we conclude that escaping particles fundamentally differ from the stable ones in their dynamic regime. In particular, stable orbits are regular (quasi-periodic), while the escaping particles show hallmarks of the deterministic chaos. Chaotic orbits generally tend to ergodically fill the allowed region, which, in our case, gets them sooner or later on the escaping trajectory. The defining property of the deterministic chaos is the sensitivity to the initial condition, which the escaping particles exhibit in a particularly spectacular manner—we observe that the orientation of the escaping particle (i.e., whether it escapes upward or downward from the equatorial plane) is sensitive to the initial condition, and even the tiniest shift of the parameters may switch the orientation of the outflow.

We have revealed the dynamic nature of both the stable and escaping trajectories, which allows us to clarify the complex structure of the escape zone. In Figure 7 we explore the particular escape zone (with ${qB}=-1$) in detail. The region inside the zone is densely populated by both the chaotic escaping orbits and the regular quasi-periodic orbits that remain oscillating in (or around) the equatorial plane. The density of escaping trajectories generally increases with the spin. We notice the red strip of stability close to the inner edge of the escape corridor, which corresponds to the main resonance of fundamental frequencies of the system (e.g., Contopoulos 2002). If we magnify this section of the escape zone, we find minor stability strips corresponding to the higher-order resonances. Self-similarity is observed in the distribution of the escaping orbits. Such a behavior is typical for the dynamical system with a nonintegrable perturbation (which is introduced by the electric charge and magnetic field in our case). Dynamics of the perturbed system may be understood using the general results of the chaos theory (e.g., Lichtenberg & Lieberman 1992; Tél & Gruiz 2006). In particular, the fundamental Kolmogorov–Arnold–Moser theorem (KAM theorem) and Poincaré–Birkhoff theorem show why some orbits resist the perturbation while others turn chaotic, and they explain the behavior of the perturbed system close to resonances.

Zoom In Zoom Out Reset image size

3.3. Acceleration of Escaping Particles

In the previous section we explored the formation of the escape zone and discussed its location in the ${r}_{0}\times a$ plane with respect to the magnetization parameter qB. Here we investigate the final velocity of escaping particles and, in particular, seek the most accelerated ones. In order to quantitatively discuss the acceleration of escaping particles of four-velocity ${u}^{\mu }$, we determine its linear velocity ${v}^{(i)}$ with respect to the locally nonrotating frame (Bardeen et al. 1972) with the tetrad basis ${e}_{\mu }^{(i)}$ as follows:

$$\begin{eqnarray}&&{v}^{(i)}=\displaystyle \frac{{u}^{(i)}}{{u}^{(t)}}=\displaystyle \frac{{e}_{\mu }^{(i)}{u}^{\mu }}{{e}_{\mu }^{(t)}{u}^{\mu }},\end{eqnarray} \tag{ 15 }$$

and we use it to calculate Lorentz factor $\gamma ={(1-{v}^{2})}^{-1/2}$, where $v=\sqrt{{[{v}^{(r)}]}^{2}+{[{v}^{(\theta )}]}^{2}+{[{v}^{(\varphi )}]}^{2}}$.

In Figure 8 we show the distribution of final values of γ in the escape zone using the color scale to distinguish particles with different γ. We observe that acceleration of the escaping particles systematically changes with the spin, as well as with the ionization radius, as suggested by the analysis in Section 3.1. In particular, in the left panel of Figure 8 we demonstrate the spin dependence of the final γ and confirm that the highest allowed spin (a = 1 for this particular case of ${qB}=-4$) leads to the highest possible acceleration. By choosing a smaller range of spin values, we also notice the radial dependence of final γ, which is generally less prominent than the former (right panel of Figure 8), and we confirm that lower ionization radius generally leads to higher acceleration. We conclude that for given qB the maximally accelerated particle (with ${\gamma }_{\max }$) corresponds to the highest allowed spin and the lowest possible ionization radius, i.e., in the presented basin-boundary plots it is always located in the upper left corner of the escape zone.

Zoom In Zoom Out Reset image size

In the following we discuss the dependence of ${\gamma }_{\max }$ on the value of qB. For each qB we iteratively integrate trajectories in the relevant segment of the escape zone in order to localize the maximally accelerated particle and obtain the value of ${\gamma }_{\max }$ with sufficient precision. In Figure 9 we observe that at first ${\gamma }_{\max }$ rises rapidly with $| {qB}|$; however, then the growth slows down, and eventually the value seems to saturate at ${\gamma }_{\max }\approx 6$. In the right panel of Figure 9 we present the same data using the logarithmic scale for $| {qB}|$, which highlights the sudden change of the shape of the function at $| {qB}| \approx 4.5$, reflecting the detachment of the escape zone from the extremal spin a = 1. Higher values of $| {qB}|$ require gradually lower spins to allow the escape of ionized particles as shown in Figure 1. On the other hand, higher values of spin generally lead to more accelerated escaping orbits as observed in Figure 8. These two mutually opposing tendencies cause the apparent slowdown in the increase of ${\gamma }_{\max }$ above the critical value of $| {qB}| \approx 4.5$.

Zoom In Zoom Out Reset image size

Increasing the magnetization parameter $| {qB}|$ shifts the whole escape zone closer to the horizon (as seen in Figure 1), and as a result, the initial radius of the maximally accelerated escaping trajectory ${r}_{\max }$ also gradually approaches the horizon as $| {qB}|$ increases (see the left panel of Figure 10). Especially for $| {qB}| \gtrsim 100$, the value of ${r}_{\max }$ becomes so close to ${r}_{+}$ that this tendency leads to growing computational difficulties in tracing the escape zone, as the Boyer–Lindquist coordinates diverge when the horizon is approached. The right panel of Figure 10 shows the value of spin ${a}_{\max }$ corresponding to the most accelerated particle. Above the threshold $| {qB}| \approx 4.5$ the value of ${a}_{\max }$ falls steeply and approaches the static limit a = 0, which it can never reach, though.

Zoom In Zoom Out Reset image size

We conclude that the maximal acceleration that can be achieved by charging the Keplerian particle within our model is ${\gamma }_{\max }\approx 6$, and it is realized in the case of high magnetization parameter, $| {qB}| \gt 100$, by particles charged near the horizon of a slowly spinning black hole, $a\lt 0.1$. For extremal black holes, a = 1, the maximal acceleration is ${\gamma }_{\max }\approx 2.5$, corresponding to ${qB}\approx -4.5$. In accordance with asymptotic behavior of the effective potential (see Equation (14)), we observe that ${\gamma }_{\max }$ increases with $| {qB}|$ and within the escape zone it grows with spin a and decreases with r₀. However, higher values of $| {qB}|$ only allow the formation of the escape zone for gradually decreasing a, and these opposing trends eventually compensate and let ${\gamma }_{\max }$ saturate at ${\gamma }_{\max }\approx 6$.

We have numerically evolved particles escaping from the magnetized accretion disk owing to the charging process occurring at inner radii of the disk. The magnetic field is known to act as a nonintegrable perturbation inducing deterministic chaos to the charged particle dynamics. However, in the given setup, the field itself never initiates the motion off the equatorial plane, which can be seen directly from the equations of motion, as well as from the general symmetry considerations. The actual impulse to leave the plane necessarily comes from a different source. In the real accretion disk, this minute yet essential dynamical fluctuation of the orbit in the direction perpendicular to the equatorial plane may be plausibly attributed to the interaction with the nearby matter. In the numerical model of the system the role of this fluctuation is played by the integration error, which is always present, although its value may be stabilized by the proper choice of the integration scheme (Kopáček et al. 2014). In particular, we manage to keep the relative error at the end of each integration below 10⁻¹⁰ (and the relative local truncation error in every iteration remains at the level $\approx {10}^{-14}$).

The crucial aspect that draws the distinction between escaping and stable oscillating orbits is the response of the trajectory to the fluctuation. If the particular trajectory is located in the chaotic domain of the phase space, then, due to its sensitivity to initial conditions and exponential growth of any deviation, it will escape if allowed by the topology of the effective potential. On the other hand, the regular trajectory does not further promote the infinitesimal disturbance and remains oscillating around the equatorial plane even if the effective potential allows the escape. In agreement with the general results of chaos theory (see, e.g., Lichtenberg & Lieberman 1992; Contopoulos 2002), we indeed observe that some parts of the phase space might be densely populated by both the chaotic escaping trajectories and the regular oscillating orbits (which is well seen especially in magnified views of the escape zones in Figures 2 and 3). Self-similar patterns are typically found in the vicinity of prominent resonances of fundamental frequencies of the system, and we detect them close to the inner edge of the escape zone (Figure 7). On the other hand, some portions of the phase space are dominated by regular motion ("islands of stability"), while some portions are mostly chaotic ("chaotic sea"). Nevertheless, even in the chaotic sea we may find small islands of stability (well seen in the right panel of Figure 8, where red-colored stable trajectories are spread over the interior of the escape zone).

Initial orbits of neutral particles were chosen as circular co-rotating geodesics (above the ISCO) and freely inspiraling geodesics (below the ISCO), i.e., they were confined to the equatorial plane ($\theta =\pi /2$). This choice was made for two main reasons. First, the analysis has an obvious astrophysical motivation, and the given initial setup maintains the connection of our model with the standard theory of thin accretion disks. Second, such a choice allows us to study the role that spin and magnetic field play in the collimation and acceleration of the outflow in the especially clear form. Particles on equatorial orbits have a zero-velocity component in the direction normal to the disk (${v}_{\perp }=0$), and even the charging process itself does not directly produce any (no kick is supposed here). Compared to the scenarios with the kick (i.e., with an ad hoc introduction of v_⊥) or initially nonequatorial orbits that already have nonzero v_⊥, the escape in our setup is less probable. We intentionally focused on this case to investigate whether the production of the outflow is possible even under these unfavorable conditions.

Without a preferred direction given by v_⊥, the orientation of the outflow (i.e., whether the particle escapes upward or rather downward from the equatorial plane) is random in the sense that it depends on the direction of the fluctuation. In the numerical representation of the system we confirm this by observing that the orientation of the escaping trajectory might change by switching the integration scheme or even by changing some of the parameters of the integration routine that affect the step size and the treatment of the local truncation error. In the astrophysical system of accreting black holes we expect the turbulent motion of nearby plasma to disturb the trajectory of charged particles, although the negligible change of the momentum accompanying the photoionization would also suffice to initiate the escape. Due to symmetry of the system, the orientation of a particular escaping trajectory becomes stochastic in a statistical sense, with the equal probability of escaping upward or downward. From the ensemble of charged particles we expect the symmetrical outflows of both orientations to be formed.

The employed model of the black hole magnetosphere was described by the electromagnetic test field. The test-field solutions are actually a limit of exact solutions for magnetized black holes that have been linearized in electromagnetic terms of the imposed Maxwell fields owing to external sources. Nevertheless, the exact electro-vacuum solution of magnetized Kerr–Newman black holes (Ernst & Wild 1976; García Díaz 1985) can be treated, in principle, by similar approaches to that adopted in the present paper (Karas & Vokrouhlický 1991). Nonetheless, the structure of these exact solutions has not yet been understood in the full generality, as they are often asymptotically curved and contain singularities outside the event horizon of the black hole; moreover, nonaxisymmetric exact solutions must be also nonstationary. However, for most astrophysically relevant situations, including the model discussed here, the linearized test-field description remains adequate.

The analysis was performed in dimensionless units where all quantities are geometrized and scaled by the rest mass of the black hole M. To retrieve their values in SI, we express the strength of the magnetic field as follows:

$$\begin{eqnarray}&&{B}_{\mathrm{SI}}=\displaystyle \frac{{{qBc}}_{\mathrm{SI}}}{{q}_{\mathrm{SI}}\left(\tfrac{M}{{M}_{\odot }}\right)1472\,{\rm{m}}},\end{eqnarray} \tag{ 16 }$$

where the quantities without subscript $\mathrm{SI}$ are dimensionless and the length $1472\,{\rm{m}}$ is the solar mass in geometrized units, ${M}_{\odot }=1472\ {\rm{m}}$. As a typical value of magnetization parameter qB for which the escape zone forms also for high spin values, we insert ${qB}=-1$.

To explore the upper limit of the specific charge ${q}_{\mathrm{SI}}$, we set the value for electrons, i.e., ${q}_{\mathrm{SI}}=-1.76\times {10}^{11}\,{\rm{C}}\,{\mathrm{kg}}^{-1}$, for which we find that for a stellar mass black hole of $M=10\,{M}_{\odot }$ the corresponding magnetic field reads ${B}_{\mathrm{SI}}=1.16\times {10}^{-7}\,{\rm{T}}$. Such a value is consistent with nonthermal filaments observed in the Galactic Center (LaRosa et al. 2004; Ferrière 2010).

Nevertheless, a typical specific charge obtained by a dust grain would be considerably lower. If we consider the supermassive black hole in the center of the giant elliptical galaxy M87 with mass $M={7.22}_{-0.40}^{+0.34}\times {10}^{9}\,{M}_{\odot }$ (Oldham & Auger 2016), we may use the estimates of the field strength derived from the Event Horizon Telescope data by Kino et al. (2014, 2015). The field at the base of the jet in M87 is constrained as $5\times {10}^{-3}\,{\rm{T}}\leqslant {\text{}}B\leqslant 1.24\times {10}^{-2}\,{\rm{T}}$ (Kino et al. 2015), while at a distance of the scale of 10 Schwarzschild radii from the source they get ${10}^{-4}\,{\rm{T}}\leqslant {\text{}}B\leqslant 1.5\times {10}^{-3}\,{\rm{T}}$ (Kino et al. 2014). Inserting these values into Equation (16), we obtain a specific charge of the particle in the range $-0.3\,{\rm{C}}\,{\mathrm{kg}}^{-1}\leqslant {\text{}}{q}_{\mathrm{SI}}\leqslant -2\times {10}^{-3}\,{\rm{C}}\,{\mathrm{kg}}^{-1}$, which is plausible for the charged dust grains. For the microquasar GRS 1915+105, which is powered by the stellar mass black hole with $M={12.4}_{-1.6}^{+2.0}\,{M}_{\odot }$ (Reid et al. 2014), the magnetic field strength in the wind zone was constrained by the upper limit in the range $0.1\,{\rm{T}}\leqslant {\text{}}B\leqslant 10\,{\rm{T}}$ (Miller et al. 2016). For these values Equation (16) yields the specific charge $-7.1\,\times {10}^{-4}\,{\rm{C}}\,{\mathrm{kg}}^{-1}\leqslant {\text{}}{q}_{\mathrm{SI}}\leqslant -5.3\times {10}^{-6}\,{\rm{C}}\,{\mathrm{kg}}^{-1}$, which suggests that ionization to lower charges (compared to the supermassive black hole in M87) is sufficient for the escape of particles from GRS 1915+105 system.

Although we do not reduce our analysis to describe a particular astrophysical object, it is apparent that the given scenario is astrophysically relevant and generally compatible with the conditions encountered in accreting black hole systems.

We have analyzed the effect of outflow of matter from the inner part of the magnetized accretion disk due to charging of initially neutral accreted dust. In particular, we have shown that a rather simple setup is sufficient to model an outflow even in the single particle approximation, where magnetohydrodynamic terms are neglected. Indeed, it appears that essential and sufficient ingredients to allow the matter to escape are the rotation of the central black hole, $a\ne 0$, and the presence of uniform magnetic field B. We assumed a basic accretion scenario in which the matter slowly sinks between circular Keplerian orbits (and becomes free-falling below the ISCO) in the equatorial plane, until it reaches charging radius r₀, where it obtains nonzero charge q. In general, the charge may stabilize the plunging orbit, as well as possibly turning the stable orbit into plunging. From the analysis of the effective potential we found that only negatively charged particles could possibly escape the attraction of a black hole with parallel orientation of spin and magnetic field.

A numerical survey revealed that escaping trajectories are actually realized only for ${qB}\lessapprox -0.5$. Until ${qB}\approx -4.5$, the escape zone stretches from extremal spin a = 1 to gradually lower spin values. However, for ${qB}\lessapprox -4.5$ it disconnects from the a = 1 limit and the whole escape zone is gradually compressed to lower spins and smaller radii as $| {qB}|$ increases. The maximal Lorentz factor ${\gamma }_{\max }$ that can be achieved by escaping particles was calculated. To some surprise, it seems to saturate at ${\gamma }_{\max }\approx 6$, which is attained with $| {qB}| \gtrapprox 100$ for $a\lessapprox 0.1$. Maximally spinning black holes accelerate the escaping matter up to ${\gamma }_{\max }\approx 2.5$ for ${qB}\approx -4.5$.

The presence of the magnetic field was crucial for the outflow for two complementary reasons. First, it changes the topology of the effective potential, so that it makes the asymptotic region in principle accessible ("potential valleys" connecting the equatorial plane with the asymptotic region are formed). The necessary condition for the existence of escaping trajectories is thus fulfilled. Another essential ingredient is, however, the impact of the field on the dynamic regime of charged particles. The magnetic field in the vicinity of a rotating black hole is known to trigger chaotic motion. And, as we suggest, only chaotic particles may escape from the circular orbits in the accretion disk, while regular ones remain oscillating in the equatorial plane even if the charge gives them enough energy to escape. This also explains why the escape zone may only form close to the black hole, where the system is sufficiently nonintegrable. Farther from the horizon, the spacetime becomes almost flat and magnetic field almost uniform, which leads to regular motion of charged particles.

Even though we considered a basic model, the analysis provides deep insight into the undisputed and crucial role of ordered magnetic field for the formation of the astrophysical jet. Most importantly, it shows that even in the particle approximation where magnetohydrodynamic terms are neglected, the outflow may be formed and attain relativistic velocities.

Authors acknowledge the support from the COST CZ program of the Czech Ministry of Education (project LD15061). O.K. is grateful for support from the Czech Science Foundation (GAČR-17-06962Y.). V.K. thanks the Czech Science Foundation for grant (GAČR-14-37086G) "Albert Einstein Center for Gravitation and Astrophysics." Discussions with Martin Kološ are highly appreciated.