Near-horizon Structure of Escape Zones of Electrically Charged Particles around Weakly Magnetized Rotating Black Hole. II. Acceleration and Escape in the Oblique Magnetosphere

We assume a thin-disk accretion geometry (Novikov & Thorne 1973) in which the electrically neutral matter gradually descends between corotating (prograde, equatorial) Keplerian orbits with energy E_Kep(r) and angular momentum L_Kep(r), given as follows (Bardeen et al. 1972):

$$\begin{eqnarray}\begin{array}{rcl}{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{array}\end{eqnarray} \tag{ 15 }$$

Circular geodesics are not allowed below the radius of the innermost stable circular orbit (ISCO), whose position for corotating particles is given as

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

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{\tfrac{3{a}^{2}}{{M}^{2}}+{Z}_{1}^{2}}$.

Below the ISCO, the geodesics are supposed to turn into freely falling inspirals maintaining the energy and angular momentum corresponding to the ISCO radius; i.e., the particles keep $E={E}_{\mathrm{Kep}}({r}_{\mathrm{ISCO}})$ and $L={L}_{\mathrm{Kep}}({r}_{\mathrm{ISCO}})$ during their infall.

However, initially neutral elements may undergo a sudden charging process (by photoionization, for example) at a given radius r₀ and obtain a nonzero specific charge q. While the change of the rest mass m may be neglected, the particle dynamics changes due to the presence of the electromagnetic field A_μ. In particular, the conserved value of energy is modified as

$$\begin{eqnarray}&&E={E}_{\mathrm{Kep}}-{{qA}}_{t},\end{eqnarray} \tag{ 17 }$$

while the values of the spatial components of canonical momentum are changed as

$$\begin{eqnarray}&&{\pi }_{r}={\pi }_{r}^{0}+{{qA}}_{r},\ \ \ {\pi }_{\theta }={{qA}}_{\theta },\ \ \ {\pi }_{\varphi }={L}_{\mathrm{Kep}}+{{qA}}_{\varphi },\end{eqnarray} \tag{ 18 }$$

where π_r⁰ is zero for particles ionized above/at the ISCO, and for infalling particles with ionization radius ${r}_{0}\lt {r}_{\mathrm{ISCO}}$, the value is calculated using the assumption that the particle is confined to the equatorial plane (${\pi }_{\theta }^{0}=0$) before it obtains the electric charge.

Once the charge is introduced, the value of the effective potential (Equation (13)) changes accordingly. In order to study the trajectories of escaping particles, we examine the behavior of the potential for r ≫ M. In particular, for the initially neutral particle with energy E_Kep ionized in the equatorial plane at r₀, we obtain the following relation valid in the asymptotic region:

$$\begin{eqnarray}&&E-{V}_{\mathrm{eff}}{| }_{r\gg M}={E}_{\mathrm{Kep}}-1-\displaystyle \frac{{{qB}}_{z}a}{{r}_{0}}+{ \mathcal O }\left({r}^{-1}\right).\end{eqnarray} \tag{ 19 }$$

Motion is possible only for E ≥ V_eff. Since E_Kep < 1 with finite r₀ and a ≥ 0 is considered, we observe that (i) particles may only escape for qB_z < 0, (ii) the escape is possible only for $a\ne 0$, (iii) the asymptotic velocity of escaping particles is an increasing function of parameters $| {{qB}}_{z}|$ and a and a decreasing function of r₀, and (iv) the escape is not allowed for the perpendicular inclination, $\alpha \equiv \arctan \left({B}_{x}/{B}_{z}\right)=\pi /2$.

Conditions (i)–(iii) are analogous to those obtained in Paper I for an aligned case with an asymptotic magnetic field B = B_z. Indeed, the perpendicular component B_x does not contribute to the energy of the charged particle in the given setup, since the corresponding terms in the time component of the vector potential (Equation (3)) vanish in the equatorial plane and the B_x terms in V_eff decrease as ${ \mathcal O }\left({r}^{-1}\right)$ or faster. Therefore, for a given asymptotic magnitude of the field $B\equiv \sqrt{{B}_{x}^{2}+{B}_{z}^{2}}$, the energy and thus the final escape velocity should decrease with increasing inclination α. For the case of the field perpendicular to the rotation axis (α = π/2), the escape is not allowed at all according to condition (iv). The effective potential for the axisymmetric system defined by Equation (6) of Paper I was more specific and also predicted the escape direction along the axis of symmetry. Although the directionality of trajectories cannot be inferred from Equation (19), we expect the particles at r ≫ M to move along the magnetic field, which is asymptotically uniform in asymptotically flat spacetime.

Analysis of the effective potential shows that escaping orbits are, in principle, possible in the given setup. The first condition, qB_z < 0, requires that the charge of the particle is negative for the parallel orientation of the spin and the magnetic field component B_z, and it is positive for the antiparallel orientation. The escaping orbits are not allowed in the Schwarzschild limit, a = 0, regardless of the values of q, B_z, B_x, and r₀. The asymptotic velocity of escaping particles rises with increasing values of a, $| q|$, and $| {B}_{z}|$ and will be higher for the particles charged closer to the horizon. The escape of particles is not allowed in the perpendicular configuration of the field (B_z = 0) regardless of the values of B_x and other parameters.

We conclude that all necessary conditions for the escape may be fulfilled in the given setup. However, as we have demonstrated in Paper I, for the special case of an axisymmetric system with B_x = 0, these conditions are not sufficient; to test whether particles really escape and discuss their asymptotic velocity, we need to employ a numerical approach and integrate the equations of motion (Equation (9)). In the rest of the paper, we switch to dimensionless units and set M = 1.

The escape zones describe regions where particles get accelerated to attain the escape velocity from the system. The emergence of escape zones and acceleration of escaping particles in the special case of an aligned magnetic field (B = B_z, B_x = 0) was discussed in Paper I. A numerical survey revealed that escaping trajectories are realized only in a limited region of parameter space. In particular, escaping trajectories were found only for qB ⪅ −0.5. In the interval −4.5 ⪅qB ⪅ −0.5, we could observe in the r₀ × a plane (ionization radius versus spin) a narrow escape zone stretching from extremal spin a = 1 to lower spin values. However, for qB ⪅ −4.5, it disconnects from the a = 1 limit, and the whole escape zone is gradually compressed to lower spin and smaller radii as $| {qB}|$ increases. For high values of spin, the escaping trajectories were allowed only for a small range of the magnetization parameter qB. The resulting acceleration of escaping particles was then limited. One of the key questions we address in the current paper is whether these limitations are also present in the oblique magnetospheres, and, in particular, whether the nonaxisymetric system could support the acceleration to ultrarelativistic velocities.

The inclination of the field breaks the axial symmetry of the system and cancels the conservation of the axial component of the angular momentum. The discussion of particle trajectories is thus enriched by new parameters compared to the axisymmetric case. The magnetic field is now specified by two parameters, namely, components B_z and B_x or, equivalently, the field magnitude $B\equiv \sqrt{{B}_{z}^{2}+{B}_{x}^{2}}$ and inclination angle $\alpha \equiv \arctan \left({B}_{x}/{B}_{z}\right)$. The initial position of the particle in the equatorial plane (θ₀ = π/2) is defined by the ionization radius r₀ and azimuthal angle φ₀. In the considered scenario, it follows from Equations (5) and (18) that ${\pi }_{\theta }^{0}(\varphi +\pi )\,=-{\pi }_{\theta }^{0}(\varphi )$, and trajectories with initial azimuthal angles φ₀ and φ₀ + π are thus equivalent.

For a given value of spin a, inclination α, and magnetization parameter qB, we numerically integrate the system of Hamiltonian equations of motion (Equation (9)) with the set of initial positions in the equatorial plane (i.e., with particles differing in initial (ionization) radius r₀ and ${\varphi }_{0}$). If the particle does not plunge into the horizon, we integrate its trajectory up to λ_fin = 1000 and consider it escaping if r_fin ≥ 200. In the plots, we use the following color-coding based on the final states: blue for the plunging orbits, red for the stable trajectories (with r_fin < 200), and yellow for the escaping particles. The typical resolution of the plots is 200 trajectories per diameter of the investigated portion of the equatorial plane. Integration of nonlinear Equation (9) must be performed with caution, as an inappropriate choice of the integration scheme might easily lead to unreliable results. To control the global integration error, we evaluate the relative error of energy ${E}_{{rr}}\equiv \left|{E}_{n}/E-1\right|$, where E_n is an actual value of energy expressed from the normalization condition ${p}^{\mu }{p}_{\mu }=-1$, which is affected by the integration errors of the coordinates and momenta, while E denotes its initial value, which is conserved by a real trajectory, since the relevant equation of motion reads $\tfrac{d{\pi }_{t}}{d\lambda }=0$. A numerical error of energy induces the artificial excitation or damping of the system, and since we deal with the nonintegrable system that is sensitive to initial conditions, the value of E_rr must be kept within reasonable bounds. Details of the employed integration routine are provided in the Appendix.

The emergence and evolution of the primary escape zone with respect to the inclination angle α for the particular choice of parameters (a = 0.98 and qB = −5) are shown in Figure 1. We intentionally choose such values of the spin and magnetization parameters for which the escape zone does not form in the aligned field. Indeed, for α = 0°, we observe no escaping orbits in the first panel of Figure 1. However, a slight perturbation of the axial symmetry with the field inclination α = 1° appears sufficient to trigger the formation of the escape zone. With increasing inclination, the escape zone becomes more pronounced and extends to higher radii. The azimuthal asymmetry of the zone becomes apparent for higher inclinations, as shown in the bottom panels of Figure 1.

Zoom In Zoom Out Reset image size

For the case of equal values of the asymptotic field components (B_z = B_x, i.e., α = 45°), the primary escape zone becomes substantially deformed. If the inclination further increases, secondary and tertiary escape zones emerge and rise in size (Figure 2). The maximum extent of the escape zones is reached for α ≈ 60°, and with higher inclinations, the zones merge and gradually decline (top panels of Figure 3). For α ≳ 80°, the escaping orbits become rare until they vanish completely for α = 84° (bottom panels of Figure 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We may conclude that the inclination of the magnetic field strongly supports the formation of the escape zones of charged particles. Analyzing the case for which there are no escaping orbits in the aligned setup, we have seen that even a very small inclination (α = 1°) is sufficient to induce the escaping trajectories; higher inclinations lead to the formation of large escape zones. However, as predicted by the analysis of the effective potential in Section 3.1, the role of the aligned component B_z is crucial, and the escape zones vanish as the inclination approaches α = π/2.

In the previous section, we discussed the formation and evolution of equatorial escape zones with respect to the inclination angle α for a fixed value of magnetization qB and spin a. In order to investigate the final velocity of escaping particles of four-velocity u^μ and seek the most accelerated ones, we determine the linear velocity v⁽ⁱ⁾ 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{ 20 }$$

We use these velocity components to express the Lorentz factor $\gamma ={(1-{v}^{2})}^{-1/2}$, where $v=\sqrt{{[{v}^{(r)}]}^{2}+{[{v}^{(\theta )}]}^{2}+{[{v}^{(\varphi )}]}^{2}}$.

The final value of γ of trajectories in the escape zone is computed for the two examples of highly inclined magnetospheres analyzed in the previous section. We use the color scale to show the distribution of values of final γ within the escape zones for α = 45° (left panel of Figure 4) and α = 54° (right panel). As predicted by the analysis of effective potential in Section 3.1, the maximal γ is found for the trajectories with minimal r₀. We observe that minimal r₀ is achieved by particles with φ₀ ≈ π/2, where the angle φ is measured counterclockwise from the positive direction of the x-axis (which is the direction of the field inclination). Comparing both panels of Figure 4, we also confirm that with a higher inclination of the field (of fixed magnitude qB), we obtain lower values of γ. While the perpendicular component B_x is crucial as a perturbation that also allows the escape in cases where it was prohibited in the parallel field and supports the formation of extended escape zones, the acceleration of escaping particles is actually controlled by the parallel component B_z (see Equation (19)).

Zoom In Zoom Out Reset image size

The maximal Lorentz factor γ_max that can be achieved by escaping particles in the aligned field was investigated in Paper I. In particular, we found that the value of γ_max saturates at γ_max ≈ 6, which is attained with $| {qB}| \gtrapprox 100$ for a ⪅ 0.1, while the maximally spinning black hole accelerates the escaping matter only up to γ_max ≈ 2.5 for qB ≈ −4.5. However, as demonstrated in the previous section, the inclination of the field also supports the formation of escape zones around rapidly spinning black holes. On the other hand, the analysis of the effective potential (Equation (13)) suggests that γ increases with B_z and fixed B; it thus decreases as α grows. Maximally accelerated particles are therefore expected from a slightly inclined magnetosphere, where the inclination α is sufficient to support the escape zone; however, the parallel component B_z, which is the actual source of energy for the acceleration, remains dominant.

For the above reasons, we restrict our search for maximally accelerated escaping particles to the small inclination of the field (we set B_x/B_z = 0.1, i.e., α ≈ 6°), extremal spin a = 1, and φ₀ = π/2 in order to support the escaping trajectories of highly accelerated particles from the small radii. The expected value of the final Lorentz factor γ may be derived from Equation (19) as

$$\begin{eqnarray}&&\gamma (a,{r}_{0},| {{qB}}_{z}| )={E}_{\mathrm{Kep}}(a,{r}_{0})+\displaystyle \frac{| {{qB}}_{z}| a}{{r}_{0}}.\end{eqnarray} \tag{ 21 }$$

However, the values of r₀ for which the escaping orbits are actually realized are not known a priori, and we have to numerically investigate the trajectories in the relevant range of initial radii to localize the inner edge of the escape zone, i.e., to identify the escaping trajectory with minimal value ${r}_{0}^{\min }$ corresponding to γ_max. Numerical analysis confirms that ${r}_{0}^{\min }$ generally decreases as $| {qB}|$ grows. In particular, for the given parameters (B_x/B_z = 0.1, ${\varphi }_{0}=\pi /2$), we observe that while $| {qB}| \approx 1$ leads to ${r}_{0}^{\min }\approx 5$, it gradually decreases to ${r}_{0}^{\min }\approx 2$ for $| {qB}| \approx 10$ and eventually stabilizes at ${r}_{0}^{\min }\approx 1.5$ for $| {qB}| \gtrapprox 100$. In Figure 5, we show the corresponding values of γ_max as a function of $| {qB}|$ and compare them with the analogous data for the aligned magnetic field from Paper I. They both coincide for small values of $| {qB}| \lessapprox 4.5$; however, then the curves split as the stronger aligned field excludes the escape from the vicinity of rapidly spinning black holes and the acceleration saturates at γ_max ≈ 6. The inclined field seems to allow the escape from the maximally spinning hole for any value of $| {qB}|$, and γ_max grows steadily with the linear trend suggested by Equation (21). The actual values of γ_max are compared with the values predicted for a given numerically obtained r₀, and they agree very well (dotted line in Figure 5). The dashed–dotted line shows the theoretical maximum of γ_max obtained by substituting a = 1 and r₀ = r₊ = 1 in Equation (21).

Zoom In Zoom Out Reset image size

Breaking the axial symmetry appears to have profound consequences regarding the acceleration of escaping particles. While the axisymmetric setup allows the acceleration up to γ ⪅ 6 (and even more strictly limited for rapidly rotating black holes), here, with the oblique magnetosphere of sufficient $| {qB}|$, we may obtain ultrarelativistic particles with γ ≫ 1. In particular, we numerically confirmed the value of γ = 680 for $| {qB}| ={10}^{3}$. Despite numerical difficulties in tracing it in strong fields, γ is supposed to further grow with increasing $| {qB}|$ and does not seem to have any fixed limit, since radiation losses (Tursunov et al. 2018) are not considered in the present analysis. The small inclination of the field thus proves sufficient to perturb the dynamics of stable orbits, and, unlike the aligned field, it admits the acceleration of particles to ultrarelativistic velocities. The efficiency of the process grows with the spin, reaching its maximum for extremal black holes; however, large γ may also be obtained for moderately rotating black holes.

The escape zones in the axisymmetric setup studied in Paper I revealed a complicated inner structure; in particular, we found regions with self-similar patterns at different magnification (Figure 7 in Paper I). The escape zone was typically composed of escaping orbits intermixed with stable orbits in a manner characteristic of objects with fractal geometry described by a noninteger dimension. This observation, together with other indications, suggests that deterministic chaos plays an important role in the escape zone, and transient chaos is a key ingredient in the dynamics of escaping particles. We expect that, compared to the previous axisymmetric case, the amount of chaos increases when the symmetry breaks, since we have previously performed the analysis of bound orbits in this system (Kopáček & Karas 2014) clearly showing the growth of chaos due to increasing inclination.

Nevertheless, chaotic episodes of unbound escaping orbits are rather short, as the particle quickly escapes far enough from the vicinity of the black hole to the region where the spacetime is almost flat and magnetic field almost uniform; thus, the nonintegrable perturbation that induces chaotic behavior diminishes there. For this reason, it becomes problematic to quantify these trajectories by the standard chaotic indicators given by the Lyapunov characteristic exponents (applied for bound orbits in Kopáček & Karas 2014), since these typically require a very long time to converge (Skokos 2010). Also, the visualization of the trajectories in the Poincaré surfaces of a section is not effective here, since the unbound motions of escaping trajectories do not provide enough intersection points. Moreover, the nonaxisymmetric system has three degrees of freedom, which further reduces the applicability of this method. Therefore, we need to employ alternative tools, e.g., the technique of recurrence analysis (Marwan et al. 2007), which is able to indicate chaotic behavior on a substantially shorter timescale regardless of the dimensionality of the system.

In Figure 6, we study the details of a particular escape zone with qB = −4.1, inclination α = 14°, and initial azimuthal angle φ₀ = 0. The escape zone is plotted in the r₀ × a plane (ionization radius versus spin) with the same color-coding of trajectories as in Figures 1–3 (blue for plunging, yellow for escaping, and red for oscillating).³ We observe the complex structure of the primary escape zone, which is characteristic of objects with fractal geometry. In particular, we notice the regions where all three types of trajectories intermix, while in the axisymmetric case, the blue region of plunging orbits is always connected and has a well-defined boundary. The inclination of the field thus substantially increases the level of complexity of the primary escape zone, which indicates the presence of strong chaos.

Zoom In Zoom Out Reset image size

In Figure 7, we explore the structure of the escape zones for the highly inclined magnetosphere (α = 54°, qB = −5), where all three classes of escape zones are present (as already shown in the last panel of Figure 2). An appropriate choice of the initial azimuthal angle (φ₀ = 60°) allows one to depict them in a single r₀ × a plot (panel (a) of Figure 7). The tertiary escape zone (panel (b)) reveals the well-defined boundary without traces of nonlinear effects. On the other hand, the edge of the secondary escape zone (panels (c) and (d)) shows a narrow fuzzy layer where plunging and stable orbits intermix, while the region of escaping orbits remains connected. Nevertheless, in the primary escape zone (panels (e) and (f)), all three types of orbits intermix, and their domains are interwoven in a complex way with a fractal structure characteristic of chaotic systems.

Zoom In Zoom Out Reset image size

Visual survey of the escape zones (Figures 6 and 7 and ) suggests that chaos plays an important role mostly in the primary escape zone (the only zone encountered with small to moderate inclinations), while there is a minor indication of chaotic dynamics on the boundary of the secondary escape zone and no traces of chaos in the tertiary zone. In the following, we employ several quantitative indicators of chaos in order to further inspect this conjecture. Since we are dealing with transient chaos of escaping trajectories, standard tools like Lyapunov exponents or Poincaré sections become ineffective. Therefore, we employ the technique of recurrence analysis (Marwan et al. 2007), which allows one to detect chaotic behavior on a short timescale. This method is based on the statistical analysis of the recurrences of the trajectory to the neighborhood of its previous states in the phase space. In particular, the binary-valued recurrence matrix  R_ij is constructed as

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{{ij}}(\varepsilon )={\rm{\Theta }}(\varepsilon -| | {\boldsymbol{x}}(i)-{\boldsymbol{x}}(j)| | )\ \ \ i,j=1,\,\ldots ,\,N,\end{eqnarray} \tag{ 22 }$$

where ε is a predefined threshold parameter, Θ is the Heaviside step function, and the vector time series ${\boldsymbol{x}}(t)$ of length N represents the analyzed segment of the trajectory in the phase space. The standard Euclidean norm L² is applied to evaluate distances after the normalization of individual time series of vector components to zero mean and unit variance, so that they contribute proportionally to the distance.

The recurrence matrix may be directly visualized in recurrence plots (RPs), which are useful for a qualitative inspection of dynamics and intuitive detection of chaos. Previously, RPs were applied in the context of relativistic astrophysics for simulated data (see, e.g., Kopáček et al. 2010; Semerák & Suková 2012; Kovář et al. 2013) and also in X-ray astronomy for the analysis of the light curves of black hole binaries (Suková & Janiuk 2016; Suková et al. 2016), while applications in gravitational-wave astronomy have been discussed recently (Lukes-Gerakopoulos & Kopáček 2018; Zelenka et al. 2020).

Various statistical measures of RPs are employed for recurrence quantification analysis (RQA), which allows a systematic survey of dynamics. A prominent feature of RP to analyze is the presence of diagonal lines. In particular, the distribution of length l of diagonal lines in an RP is given by the histogram $P(\varepsilon ,l)$ as

$$\begin{eqnarray}\begin{array}{rcl}P(\varepsilon ,l) & = & \displaystyle \sum _{i,j=1}^{N}(1-{{\boldsymbol{R}}}_{i-1,j-1}(\varepsilon ))\\ & & \times \,(1-{{\boldsymbol{R}}}_{i+l,j+l}(\varepsilon ))\displaystyle \prod _{k=0}^{l-1}{{\boldsymbol{R}}}_{i+k,j+k}(\varepsilon ).\end{array}\end{eqnarray} \tag{ 23 }$$

The average length (AVL) of the diagonal lines is evaluated as

$$\begin{eqnarray}&&\mathrm{AVL}\equiv \displaystyle \frac{{\displaystyle \sum }_{l={l}_{\min }}^{{l}_{\max }}{lP}(\varepsilon ,l)}{{\displaystyle \sum }_{l={l}_{\min }}^{{}_{\max }}P(\varepsilon ,l)},\end{eqnarray} \tag{ 24 }$$

where only diagonal lines of length at least l_min count, and l_max is the length of the longest line (except the line of identity on the main diagonal).

The recurrence indicator ENTR is defined as the Shannon entropy of the probability $p(\varepsilon ,l)=P(\varepsilon ,l)/{N}_{l}$ of finding a diagonal line of length l in an RP,

$$\begin{eqnarray}&&\mathrm{ENTR}\equiv -\displaystyle \sum _{l={l}_{\min }}^{{l}_{\max }}p(\varepsilon ,l)\mathrm{ln}p(\varepsilon ,l),\end{eqnarray} \tag{ 25 }$$

where N_l is the total number of diagonal lines: ${N}_{l}(\varepsilon )\,={\sum }_{l\geqslant {l}_{\min }}P(\varepsilon ,l)$.

Other RQA measures based on the presence of diagonal and vertical lines in an RP may be defined; however, here we only employ the indicators AVL and ENTR. Definitions and properties of other RQA indicators may be found in a review by Marwan et al. (2007).

Besides the recurrence analysis, we employ a standard method for characterizing fractal sets in n-dimensional metric space, and we evaluate the Minkowski–Bouligand (box-counting) dimension DIM, defined as

$$\begin{eqnarray}&&\mathrm{DIM}=\mathop{\mathrm{lim}}\limits_{\delta \to 0}\displaystyle \frac{\mathrm{log}P(\delta )}{\mathrm{log}(1/\delta )},\end{eqnarray} \tag{ 26 }$$

where P(δ) is the number of n-dimensional boxes with side δ required to cover the set. For a complete definition of the box-counting dimension, its properties, and its relation to other fractal dimensions, see, e.g., Falconer (2003). The box-counting dimension of the fractal escape boundaries shown in Figures 6 and 7 could be evaluated; however, we calculate the indicator of the trajectories directly in order to compare their values within the escape zones. A similar application of the box-counting dimension and other chaotic indicators was recently presented by Pánis et al. (2019).

We compute the chaotic indicators AVL, ENTR, and DIM for the same set of orbits as presented in Figure 4 (i.e., for parameters qB = −5 and a = 0.98 with inclinations α = 45° and 54°, respectively); however, we use a shorter segment of the trajectory (λ_fin = 300 instead of 1000), and the trajectories are now integrated with the fixed time step as required to perform recurrence analysis correctly. Only a radial coordinate of the trajectory is employed for the analysis; a single time series r(λ) with 1000 data points is used as a input. The box-counting dimension is computed with the BoxCountfracDim function,⁴ and RQA indicators are evaluated with the CRP Toolbox (Marwan et al. 2007). The threshold parameter is set to the fixed value ε = 0.2, the minimal length of the diagonal line l_min = 2, and embedding is not used. All calculations are performed in Matlab R2014a.

In Figure 8, we present color-coded values of indicators of trajectories within the escape zones. Regarding the case of inclination α = 45° (upper panels of Figure 8), we observe that the escape zone is clearly distinguished by the values of the indicators. Moreover, the values of the box-counting dimension DIM (top right panel) to some degree correspond with the values of the final Lorentz factor (left panel of Figure 4). Highly accelerated trajectories tend to have higher DIM, and vice versa. The values of the RQA indicators AVL and ENTR do not entirely follow this trend; however, they still clearly distinguish the escaping orbits within the escape zone (especially in the case of AVL). The case of α = 54°, where all three classes of escape zone are present, is analyzed in the bottom panels of Figure 8. Comparison with the right panel of Figure 4 shows that while the shape of the primary escape zone is clearly resolved, only the inner parts of the secondary escape zone are recognized by the increased values of the indicators, and the shape of the tertiary escape zone cannot be distinguished.

Zoom In Zoom Out Reset image size

The behavior of the chaotic indicators of the trajectories in the escape zones thus corresponds with the observations of the detailed structure of the zones made in Figures 6 and 7. In particular, the comparison with Figure 8 confirms the conjecture that transient chaos plays a dominant role in the dynamics of the primary escape zone. The fractal geometry of this zone directly corresponds with the increased values of chaotic indicators. In the case of the secondary escape zone, only a narrow fuzzy layer of intermixed trajectories was observed on its edge (panels (c) and (d) of Figure 7), corresponding to increased values of the chaotic indicators in the inner part of the zone. On the other hand, regular dynamics in the tertiary escape zone with well-defined boundaries and no intermixing of trajectories is also confirmed by chaotic indicators, as their values do not increase within this zone.