Determination of chaotic behaviour in time series generated by charged particle motion around magnetized Schwarzschild black holes

TL;DR

This paper investigates the chaotic behavior of charged particles in a magnetized Schwarzschild black hole environment, using various chaos detection methods including machine learning, to understand how magnetic and geometric parameters influence particle dynamics.

Contribution

It introduces a comprehensive analysis of chaos in charged particle motion around black holes, employing machine learning techniques to identify chaotic regimes more efficiently.

Findings

Chaotic behavior increases with the inclination angle.

Machine learning outperforms traditional methods in detecting chaos.

High inclination angles lead to more chaotic particle trajectories.

Abstract

We study behaviour of ionized region of a Keplerian disk orbiting a Schwarzschild black hole immersed in an asymptotically uniform magnetic field. In dependence on the magnetic parameter ${\cal B}$, and inclination angle $\theta$ of the disk plane with respect to the magnetic field direction, the charged particles of the ionized disk can enter three regimes: a) regular oscillatory motion, b) destruction due to capture by the magnetized black hole, c) chaotic regime of the motion. In order to study transition between the regular and chaotic type of the charged particle motion, we generate time series of the solution of equations of motion under various conditions, and study them by non-linear (box counting, correlation dimension, Lyapunov exponent, recurrence analysis, machine learning) methods of chaos determination. We demonstrate that the machine learning method appears to be the most efficient in determining the chaotic region of the $\theta-r$ space. We show that the chaotic character of the ionized particle motion increases with the inclination angle. For the inclination angles $\theta \sim 0$ whole the ionized internal part of the Keplerian disk is captured by the black hole.

Tables

Table 1. Table 1: Magnitudes of magnetic fields B 𝐵 B in Gauss units for magnetic parameter ℬ = 0.01 ℬ 0.01 {\cal{B}}=0.01 and stellar M = 10 ​ M ⊙ 𝑀 10 subscript 𝑀 direct-product M={10\leavevmode\nobreak\ {M}_{\odot}} and supermassive M = 10 8 ​ M ⊙ 𝑀 superscript 10 8 subscript 𝑀 direct-product M={10^{8}\leavevmode\nobreak\ {M}_{\odot}} black holes. Test particles with different specific charges q / m 𝑞 𝑚 q/m ratios are assumed - from electron and proton to charged dust grain (one electron lost, m = 2 × 10 − 18 𝑚 2 superscript 10 18 m=2\times 10^{-18}  kg).

$ℬ=0.01$	electron	proton	charged dust
$M=10M_{\odot }$	$2\cdot {10}^{-5}G$	$0.04G$	$2\cdot {10}^{8}G$
$M={10}^8{M}_{\odot }$	$2\cdot {10}^{-12}G$	$4\cdot {10}^{-9}G$	$20G$

Table 2. Table 2: Comparison of values of different nonlinear methods applied on regular, chaotic and pseudo-random generated series of the length 10 000.

Met./ Type	Regular	Weak ch.	Strong ch.	P.-rand.
Box-count	0.0969	0.8343	1	1
Corr. dim.	$5\ast {10}^{-15}$	0.8943	0.8735	0.9978
Lyap. exp.	$-7\ast {10}^{-16}$	0.1889	0.3289	0.0842
RQA-RR	0.4999	0.0183	0.0116	0.0057
RQA-DET	1	0.6959	0.6114	0.0113
RQA-ENTR	8.5152	0.8951	0.8202	0.0296
RQA-LAM	5001	5.6299	2.8507	2.0006
ML-Rand. For.	0.5915	0.6917	0.7869	0.3788

Table 3. Table 3: Comparison of time in seconds required for different nonlinear methods applied on time series generated by logistic map. For parameter r 𝑟 r varying from 2 2 2 to 4 4 4 with the step of 0.01 0.01 0.01 leads to 200 200 200 time series, with initial value x 0 = 0.1 subscript 𝑥 0 0.1 x_{0}=0.1 we did set up the iterations to length 100 , 1000 100 1000 100,1000 and 10000 10000 10000 and compared the time needed for the calculation for given methods.

Method / Length	100	1000	10000
Box-count	0.354	0.416	3.519
Correlation dim.	4.155	26.54	2353
Lyapunov exp.	3.294	119.5	12574
RQA- RR	0.164	6.332	507.8
RQA-DET	0.209	6.069	916.9
RQA-ENTR	0.197	7.247	1128
RQA-LAM	0.197	7.248	1128
ML-Rand. For.	0.677	1.212	7.479

Equations

$$\mathrm{d}{s}^2=-f(r)\mathrm{d}{t}^2+f^{-1}(r)\mathrm{d}{r}^2+r^2(\mathrm{d}{\theta }^2+{\sin }^2\theta \mathrm{d}{\varphi }^2),$$

$$f(r)=1-\frac{2M}{r}.$$

$$A_{\varphi }=\frac{B}{2}{g}_{\varphi \varphi }=\frac{B}{2}{r}^2{\sin }^2\theta .$$

$$H_{\mathrm{p}}=\frac{1}{2}{g}^{\alpha \beta }({\pi }_{\alpha }-q{A}_{\alpha })({\pi }_{\beta }-q{A}_{\beta })+\frac{1}{2}{m}^2,$$

$${\pi }^{\mu }=p^{\mu }+q{A}^{\mu },$$

$$\frac{\mathrm{d}{x}^{\mu }}{\mathrm{d}\zeta }\equiv p^{\mu }=\frac{\partial H}{\partial {\pi }_{\mu }},\frac{\mathrm{d}{\pi }_{\mu }}{\mathrm{d}\zeta }=-\frac{\partial H}{\partial x^{\mu }}.$$

$$E$$

$$L$$

$$x=r\cos (\varphi )\sin (\theta ),y=r\sin (\varphi )\sin (\theta ),z=r\cos (\theta ).$$

$$\mathcal{E}=\frac{E}{m},\mathcal{L}=\frac{L}{m},\mathcal{B}=\frac{qB}{2m},$$

$$H=\frac{1}{2}f(r)p_r^2+\frac{1}{2r^2}{p}_{\theta }^2+\frac{1}{2}\frac{m^2}{f(r)}\left[V_{\mathrm{eff}}(r,\theta )-{\mathcal{E}}^2\right],$$

$$V_{\mathrm{eff}}(r,\theta )\equiv f(r)\left[1+{\left(\frac{\mathcal{L}}{r\sin \theta }-\mathcal{B}r\sin \theta \right)}^2\right].$$

$${\mathcal{E}}^2=V_{\mathrm{eff}}(r,\theta ;\mathcal{L},\mathcal{B}).$$

$$0=q_2+q_3,$$

$${\pi }_{\alpha (1)}={\pi }_{\alpha (2)}+{\pi }_{\alpha (3)}$$

$$p_{\alpha (1)}=p_{\alpha (2)}+p_{\alpha (3)}.$$

$$p_{\alpha (1)}\approx p_{\alpha (2)}\gg p_{\alpha (3)}.$$

$$m_{(\mathrm{I})}=m_{(\mathrm{II})},p_{(\mathrm{I})}^{\mu }=p_{(\mathrm{II})}^{\mu }.$$

$${\pi }_{\alpha (\mathrm{II})}=p_{\alpha (\mathrm{I})}+q{A}_{\alpha }.$$

$$x^{\alpha }$$

$$u_{\alpha }$$

$${\mathcal{L}}_{(\mathrm{I})}=\frac{r_0\sin {\theta }_0}{\sqrt{r_0-3}},{\mathcal{E}}_{(\mathrm{I})}=\frac{r_0-2}{\sqrt{r_0^2-3r_0}}.$$

$${\mathcal{L}}_{(\mathrm{II})}={\mathcal{L}}_{(\mathrm{I})}+\mathcal{B}{r}_0^2{\sin }^2{\theta }_0,{\mathcal{E}}_{(\mathrm{II})}={\mathcal{E}}_{(\mathrm{I})}.$$

$$B=\frac{2m{c}^4\mathcal{B}}{qGM}\text{leavevmode}[\mathrm{G}],$$

$$D_0=\underset{ϵ\to 0}{\lim }\frac{\ln N(ϵ)}{\ln \frac{1}{ϵ}},$$

$$D_2=\underset{ϵ\to 0}{\lim }\frac{\ln C(ϵ)}{\ln ϵ},$$

$$D_q=\underset{ϵ\to 0}{\lim }\frac{1}{q-1}\frac{\ln {\sum }_k{p}_k^q}{\ln ϵ}-\infty \le q\le \infty ,$$

$$\hat{C}(ϵ)=\underset{N\to \infty }{\lim }\frac{2}{N(N-1)}\sum _{i<j}H(ϵ-|x_i-x_j|),$$

$$x_i=(x_{i-m+1},x_{i-m+2},...,x_i)\in {\mathbb{R}}^m.$$

$${\lambda }_m=\underset{N\to \infty }{\lim }\frac{1}{N}{\sum _{n=0}}^{N-1}\ln [|f^{\prime }(f^n(x_0))|].$$

$$R_{i,j}=H(ϵ-\Vert x_i-x_j\Vert \mid )i,j=1,...,N,$$

Full text

black holes, magnetic fields, charged particle motion, time series, chaos and regularity, Lyapunov exponent, correlation dimension, machine learning

1. I Introduction

2. II Charged particle motion around magnetized Schwarzschild black holes

3. III Ionization of the Keplerian disks

4. III.1 Initial conditions of the motion of particles forming the ionized disk

5. III.1.1 Ionization scenarios

6. III.1.2 Initial conditions of the charged particle motion

7. IV Motion of charged particles forming ionized disk

8. IV.1 Typical trajectories of charged particles: chaos and regularity

9. IV.2 Non-linear methods determining chaos-regularity transition

10. IV.2.1 Description of the numerical methods

11. IV.2.2 Evolution of ionized Keplerian accretion disk

12. IV.3 Astrophysical relevance

13. V Conclusions

14. A Methods of determination of chaotic behaviour in time series

15. A.1 Methods detecting chaos

16. A.1.1 Box-counting

17. A.1.2 Correlation dimension

18. A.1.3 Lyapunov exponent

19. A.1.4 Recurrence quantification analysis

20. A.1.5 Machine learning

21. A.2 Test of the non-linear methods detecting chaos

I Introduction

It is well known that the test particle motion in the field of Kerr black holes is fully regular (Carter, 1968). The same is true for both charged and uncharged particles moving in the field of Kerr-Newman black holes with an internal electromagnetic field related to their electric or tidal charge (Misner et al., 1973; Bičák et al., 1989; Blaschke and Stuchlík, 2016). The regularity holds also for Kerr-Newman-de Sitter black holes (Carter, 1973), even if they are dyonic, i.e., carrying both electric and magnetic charges (Stuchlík, 1983). The situation changes radically for black holes immersed in an external electromagnetic field since the test charged particle motion takes generally chaotic character Karas and Vokrouhlický (1992); Stuchlík and Kološ (2016).

Most of the observed black hole candidates have an accretion disk constituted from conducting plasma which dynamics can generate a magnetic field external to the black hole. Another possibility is represented by an external galactic magnetic field that can be amplified by the black hole strong gravity. Such magnetic fields satisfy the test field approximation condition, having thus negligible effect on the spacetime structure and the motion of neutral particles, however, for particle with large specific charge the electromagnetic Lorentz force is relevant and leads generally to the chaotic motion of charged particles Frolov and Shoom (2010); Kovář et al. (2008); Kološ et al. (2017).

The exact shape and structure of the magnetic fields around compact objects is still under examination, but the uniform magnetic field assumption introduced by Wald Wald (1974) can be used as first simple approximation to more complex fields. The charged test particle motion in such an asymptotically uniform configuration has been already studied in a large variety of papers that give significant insight into the astrophysical processes in vicinity of magnetized black holes Prasanna and Vishveshwara (1978); Prasanna (1980); Aliev and Galtsov (1981); Frolov and Shoom (2010); Kovář et al. (2010); Kopáček et al. (2010); Abdujabbarov et al. (2013a); Al Zahrani et al. (2013); Abdujabbarov et al. (2013b); Kopáček and Karas (2014); Shiose et al. (2014); Kološ et al. (2015); Tursunov et al. (2016); Stuchlík and Kološ (2016); Kopáček and Karas (2018); Benavides-Gallego et al. (2019). In the present paper, we examine chaotic charged test particle dynamics around a Schwarzschild black hole immersed in an external uniform magnetic field, resulting under special initial condition of ionized Keplerian accretion disk.

The matter forming an electrically neutral (Keplerian) accretion disk orbiting such a magnetized black hole can get ionized (e.g. by an irradiation), and start to feel the external magnetic field. Under the influence of the magnetic field, the original purely circular motion of the electrically neutral matter has to be transformed into one of the following regimes of the motion of created charged particles (Stuchlík and Kološ, 2019): a) regular oscillatory motion possibly reflecting the high-frequency X-ray quasiperiodic oscillations observed in microquasars Kološ et al. (2017), b) destruction of the ionized region of the disk due to the radial infall into the black hole, c) chaotic motion governing transformation of the Keplerian disk into thick toroidal structure, in combination with creation of winds (or relativistic jets for the case of rotating black holes Stuchlík and Kološ (2016)).

We explore charged particle destiny for a variety of initial conditions, namely the electromagnetic interaction intensity parameter reflecting intensity of the magnetic field and the specific charge of the particle ${\cal{B}}$, initial position of the particle given by the radius $r$, and the inclination angle of the disk to the magnetic field lines $\theta$. We concentrate our attention on the transition between the regular and chaotic character of the motion of the resulting charged particle dynamics in dependence on the initial conditions, namely for fixed values of the electromagnetic interaction parameter ${\cal{B}}$, we determine the distribution of the regular and chaotic states in the plane $\theta-r$ of the initial conditions; the measure of the chaos is also reflected.

The equations of charged particle motion around magnetized black holes are characterized generally as a system demonstrating deterministic chaos – its determination is not trivial and demands application of non-linear methods as the box-counting on one side, or the machine learning on the other.

Observations/measurements of a (physical or general) quantity in time produce sequences of numbers, and hence time sequence (series) analysis is important tool in almost every corner of current science. Non-linear systems can produce time series with very complex behaviour that are on the first sight similar to the series of random numbers. Non-linear deterministic systems can demonstrate chaotic behaviour that is only apparently random behaviour. The deterministic chaos is hard to distinguish from random noise in observed data, especially when the degree of freedom of the generating non-linear system is high.

Very well established linear tools, like Fourier analysis, can easily distinguish between regular and random sequences of numbers, but they fail to distinguish between a deterministic-chaos sequence and a random sequence, giving flat (constant) power spectral density in both cases. An observer using the Fourier analysis for analysis of chaotic data could make wrong conclusion and claim the data are completely random with no information inside. If we would like to extract any information from chaotic time series, the non-linear tools of detection should be used. In the present paper, we test behaviour of real-number sequences using non-linear tools – namely, the box-counting, correlation dimension, Lyapunov exponent, recurrence quantification analysis (RQA) and machine learning. These non-linear tools are shortly presented in Appendix, and tested using the simple case of toy model based on the simple logistic map for regular/chaotic sequence generation. All codes of the methods presented in this article are written in the Wolfram Mathematica 11, since the Mathematica software provides the high-level programming language with already implemented subroutines (functions) for machine learning. The sequences of the tested regular/chaotic data are generated by the solution of the motion equations of the charged particle motion with various initial conditions.

Combination of deterministic chaos and random noise in real measured data, i.e. small signal to noise ratio, can bring new problems to the task of distinguishing between chaos/noise. In our study such a problem is avoided, as we test the non-linear methods of detecting the deterministic chaos on time series of regular/chaotic data generated theoretically by the solution of the equations of the charged particle motion.

The present paper can be considered as a preliminary study where we test the non-linear chaos detecting methods on system which can be controlled and where the dynamics is known. In future we plan to use such methods to study the timing signals (photon number time dependence) related to the quasiperiodic oscillations observed in the microquasar sources Remillard and McClintock (2006). The non-linear test of the observed microquasar timing data has been already applied in different context Karak et al. (2010); Mannattil et al. (2016); Suková et al. (2016); Huppenkothen et al. (2017).

Throughout the paper, we use the spacetime signature $(-,+,+,+)$, and the system of geometric units in which $G=1=c$. Greek indices are taken to run from 0 to 3.

II Charged particle motion around magnetized Schwarzschild black holes

The line element of the Schwarzschild black hole spacetime with mass $M$ reads

[TABLE]

where the so called lapse function $f(r)$ takes the form

[TABLE]

Hereafter, we put for simplicity $M=1$, i.e. we use dimensionless radial coordinate $r$ (and time coordinate $t$).

We assume an asymptotically uniform magnetic field having strength $B$ at the spatial infinity (i.e. at large distances from the black hole). The magnetic field lines are oriented perpendicularly to the equatorial plane of the black hole spacetime. The only non-zero covariant component of the potential of the electromagnetic field takes the form Wald (1984)

[TABLE]

The Hamiltonian governing the test (non-radiating) charged particle motion can be written in the form Wald (1984)

[TABLE]

where the kinetic four-momentum $p^{\mu}=mu^{\mu}$ is related to the generalized (canonical) four-momentum $\pi^{\mu}$ by the relation

[TABLE]

that satisfy the Hamilton equations presented in the form

[TABLE]

The affine parameter $\zeta$ is related to the proper time $\tau$ of the particle by the relation $\zeta=\tau/m$.

Due to the symmetries of the Schwarzschild spacetime (1) and the asymptotically uniform magnetic field (3), one can easily find the conserved quantities of the particle motion – the energy and the axial angular momentum that can be expressed as

[TABLE]

The dynamical equations for the charged particle motion in the Cartesian coordinates can be found using the coordinate transformations

[TABLE]

Introducing for convenience the specific constants of the motion, energy ${\cal{E}}$, axial angular momentum ${\cal{L}}$, and the magnetic parameter ${\cal{B}}$ governing the intensity of the electromagnetic interaction, by the relations Frolov and Shoom (2010); Kološ et al. (2015)

[TABLE]

one can rewrite the Hamiltonian (4) in the form

[TABLE]

where $V_{\rm eff}(r,\theta;{\cal{L}},{\cal{B}})$ denotes the effective potential governing the turning points of the radial and latitudinal motion, given by the relation

[TABLE]

The terms in the parentheses corresponds to the central force potential given by the specific angular momentum ${\cal{L}}$, and electromagnetic potential energy specified by the magnetic parameter ${\cal{B}}$.

The effective potential (12) demonstrates clear symmetry $({\cal{L}},{\cal{B}})\leftrightarrow(-{\cal{L}},-{\cal{B}})$ that allows to distinguish in the following only two situations corresponding to the minus configuration (${\cal{L}}>0,{\cal{B}}<0$) where the magnetic parameter and the axial angular momentum have opposite signs and the Lorentz force is attracting the charged particle to the $z$-axis towards the black hole, and to the plus configuration (${\cal{L}}>0,{\cal{B}}>0$) where magnetic parameter and the axial angular momentum have the same signs and the Lorentz force is repulsive, acting outward the black hole. The positive angular momentum of a particle, ${\cal{L}}>0$, means that the particle is revolving in the counter-clockwise motion around the black hole. If charge of the particle is considered to be positive, $q>0$, the minus configuration, ${\cal{B}}<0$, corresponds to the vector of the magnetic field $\vec{B}$ pointing downwards, while plus configuration, ${\cal{B}}>0$, corresponds to the vector of the magnetic field $\vec{B}$ pointing upwards the $z$-axis Kološ et al. (2015).

The charged particle motion is limited by the energetic boundaries given by

[TABLE]

The axial symmetry of the background of the combined gravitational and magnetic fields implies independence of the effective potential $V_{\rm eff}$ from the coordinate $\phi$ which allows us to examine $V_{\rm eff}(r,\theta)$ as a 2D function of the spherical coordinates $r,\theta$, or Cartesian $x,z$ coordinates (9). The effective potential is positive outside the black hole horizon, and diverges at the horizon $r=2$. The region within the horizon is excluded from our investigation.

The effective potential (12) enables us to demonstrate general properties of the charged particle dynamics and has already been explored in detail in Kološ et al. (2015) – here we are directly applying these results.

III Ionization of the Keplerian disks

For the magnetized Schwarzschild black holes the spherical symmetry of the spacetime is combined with the axial symmetry of the uniform magnetic field. If we assume the lines of the uniform magnetic field oriented in the direction of the $z$-axis (vertical direction) Wald (1974), the equatorial plane of the spacetime has to be perpendicular, i.e., the $x-y$ plane – see Fig. 2.

III.1 Initial conditions of the motion of particles forming the ionized disk

In the initial state we consider a Keplerian accretion disk orbiting the magnetized Schwarzschild black hole that consists from electrically neutral test particles following circular geodesics. Due to the spherical symmetry of the Schwarzschild spacetime the Keplerian disk can be located in any central plane of the spacetime, being thus inclined to the equatorial plane of the magnetized Schwarzschild black hole by a latitudinal angle $\theta_{0}$, giving one of the initial conditions for the motion of particles of the ionized Keplerian disk, see Fig. 2.

III.1.1 Ionization scenarios

In order to demonstrate the influence of the magnetic field on an ionized Keplerian accretion disk, some realistic ionization scenarios for the neutral particles of the disk have to be considered. As ionization model of a neutral particle, one can consider the Magnetic Penrose Process (MPP) Parthasarathy et al. (1986); Dadhich et al. (2018), where the original 1st neutral particle splits into two charged particles - 2nd and 3rd. Conservation of total charge of particles entering the ionization process takes the form

[TABLE]

and the law of conservation of the canonical momenta (5)

[TABLE]

takes due to the charge conservation the form

[TABLE]

In many realistic scenarios, like neutron $\beta$ decay or neutral atom ionization, one of the created charged particles is much more massive then the other one, $m_{2}/m_{3}\gg\leavevmode\nobreak\ 1$ – as follows from the proton/electron or ion/electron mass ratios. The more massive charged particle (proton or ion) takes almost all the initial momentum of the original neutral particle, and the dynamical influence of the lighter charged particle (electron) can be neglected

[TABLE]

In another realistic scenario one can consider the Keplerian accretion disk created by a plasma that can be considered as a quasi-neutral soup of charged particles - electrons and ions - orbiting around the central black hole along circular geodesic orbits. If the accretion disk is dense enough, the main free path of the charged particles will be very short in comparison to the length of the orbit around the central hole. This means that the influence of the magnetic field on the charged particle motion is effectively suppressed, and they move as a composite neutral body. However, when the density of plasma of such accretion disks significantly decreases, the charged particles start to feel effectively the influence of the magnetic field that modifies substantially their trajectories. Of course, such kind of effective ionization of the Keplerian disks is restricted to the regions located near their inner edge corresponding to the marginally stable circular geodesic (ISCO).

Both previously mentioned scenarios enable a simplified ionization model, where the neutral particle just obtains an electric charge while its mechanical momentum remains conserved. Such an ionization model (17) has been already studied, but in the field of rotating Kerr black hole Stuchlík and Kološ (2016); Kopáček and Karas (2018), where the charged particle escape velocities and structure of escape zones were explored. Here we focus our attention on the complementary study of measure of chaos of the ionized particle motion in dependence on the particle initial conditions related to the inclination of the Keplerian disk, and the particle initial position (radius). We test which initial conditions lead to destruction of the Keplerian disk, or to transitions between the regular and chaotic motion.

III.1.2 Initial conditions of the charged particle motion

We thus study the simplified particle ionization model (17) that considers only the heavy particles in the Magnetic Penrose Process, and can be characterized by simple conditions of the mass conservation and kinetic momentum conservation

[TABLE]

The test particle kinetic momenta before (I) and after (II) ionization are thus equal and the canonic momentum of the ionized particle is given by the relation

[TABLE]

At the moment of ionization, the neutral or charged test particle is considered to be located on inclined circular orbit with initial position $x^{\alpha}$ and initial four-velocity $u_{\alpha}$

[TABLE]

The specific angular momentum ${\cal{L}}$ and the specific energy ${\cal{E}}$ of the neutral test particles on the inclined circular orbits of the Keplerian accretion disk are given by the relations Wald (1984)

[TABLE]

The simplified ionization condition (19) and use of the definitions of the energy (7) and the angular momentum (8) enable to write the formulas of the ionized test particle specific angular momentum ${\cal{L}}$ and specific energy ${\cal{E}}$ in the form

[TABLE]

The magnetic field has in the Schwarzschild metric only one non-vanishing component, $A_{\phi}$ (3), hence only the (canonical) specific angular momentum ${\cal{L}}$ is changed due to the ionization, while the particle specific energy ${\cal{E}}$ remains constant. The energy of the neutral particles on the circular geodesic orbits, given by (22), is always ${\cal{E}}\leq\leavevmode\nobreak\ 1$, but energy ${\cal{E}}>1$ is needed for the charged particle to escape to infinity along the $z$-axis. Therefore, no escape to infinity of ionized matter from the Keplerian disks is possible in the Schwarzschild spacetime Stuchlík and Kološ (2016).

IV Motion of charged particles forming ionized disk

Since in the Schwarzschild metric the ionized particle following originally a circular geodesic cannot escape to infinity, the capture by the black hole, or bound motion in vicinity of the original circular orbit are the only options. If the charged particle is not captured, its motion remains bounded in some closed region around the black hole – the motion of such a bounded charged particle is in general chaotic. The character of the bounded motion depends mainly on the inclination of the original Keplerian disk, and on the magnitude of the electromagnetic interaction. As we shall see, for large inclination angles of the Keplerian disks to the magnetic field lines, $\theta_{0}\sim\pi/2$, the bounded motion could be regular.

Our study of the simplified ionization of the Keplerian disks is related to the heavier charged particles resulting from the considered scenarios of the MPP process, i.e., to the motion of protons and ions, not to the electron motion in general. However, we have to note that in some special situations our results could be applied with reasonable precision also for the electron motion – these situations have to correspond to the case when the electron has initially the specific energy and specific angular momentum close to those of the circular geodesic motion. Such a condition can be clearly satisfied for some ionized atoms, or while the electron following the circular geodesics starts to feel the influence of the magnetic field due to decreasing density of the disk.

IV.1 Typical trajectories of charged particles: chaos and regularity

Due to the spherical symmetry of the Schwarzschild spacetimes, the Keplerian disk can be located in any central plane, contrary to the case of the Kerr black hole spacetimes where the inner parts of the Keplerian disks have to be located in the equatorial plane of the spacetime due to the so called Bardeen-Peterson effect reflecting the interplay of the disk viscous stresses and the frame dragging of the spacetime Bardeen and Petterson (1975). We thus first discuss properties of the ionized disks orbiting the magnetized Schwarzschild black hole with large enough initial inclination angle, $\theta_{0}\sim\pi/2$, in dependence on the magnetic parameter characterizing the intensity of the electromagnetic interaction. As a second case we discuss the role of the inclination angle for the case of the magnetic interaction parameter giving the strongest chaotic behaviour for the near-equatorial disks. Finally, we study the combined role of the magnetic parameter ${\cal{B}}$ and the inclination angle for charged particles located initially at the ISCO.

The influence of the magnitude of the magnetic field (parameter) on the character of the ionized near-equatorial accretion disk is demonstrated in Fig. 3. When the magnetic field is missing (${\cal{B}}=0$ case), all the orbits of the ionized Keplerian disk keep their circular shape and we see just inclined razor thin disk. When a slight magnetic field is switched-on (${\cal{B}}=\pm 0.001$ cases), the charged particles forming the ionized disk become to be unsettled, and the circular orbits are perturbed, but the ionized particle orbits are still bounded in vicinity of the equatorial plane, following oscillatory motion in the radial and latitudinal coordinates, and the accretion disk becomes slightly thick. Notice that in the case ${\cal{B}}=-0.001$ the particles starting very close to the ISCO are captured by the black hole, while such trapping is not possible for the case ${\cal{B}}=0.001$. When a larger magnetic field is switched-on (${\cal{B}}=\pm 0.01$ cases), the charged particle trajectories become chaotic and a lot of trajectories finishes inside the black hole. The complete destruction of the ionized Keplerian disk can be seen for the ${\cal{B}}=-0.01$ case, when all the orbits are captured by black hole. For even larger magnetic field (${\cal{B}}=\pm 0.1$ cases), the Lorentz force becomes more influential and charged particle trajectories start to be strongly chaotic creating a toroidal structure, and demonstrate tendency to wind up and down along the magnetic field lines, especially in the case of ${\cal{B}}=-0.1$ when the structure of orbiting particles resembles a sphere. When the magnetic field parameter is large (${\cal{B}}=\pm 1$ cases), the Lorentz force becomes to be the dominant force for the particle motion. The charged particles wind up and down along magnetic field lines, demonstrating simultaneously slow shifting around the black hole in the clockwise (${\cal{B}}>0$) or counter-clockwise (${\cal{B}}<0$) direction. Once again, the thin Keplerian disk has been transformed into some thick structure spread around the equatorial plane of the magnetized black hole. For larger values of the magnetic parameter, $|{\cal{B}}|>1$, similar behaviour is demonstrated by the ionized Keplerian disk.

The influence of the initial inclination angle $\theta_{0}$ of the Keplerian accretion disk to the magnetic field lines on the character of the motion of particles of the ionized disk is demonstrated in Fig. 4 for the magnetic parameter fixed to the value of ${\cal{B}}=0.1$ corresponding to the case when the ionized near-equatorial Keplerian disks demonstrate the most chaotic behaviour. Our goal is to demonstrate clearly how the modifications of the inclination angle influence the chaotic character of the particle motion of the ionized disk.

Clearly, for large inclinations angle of a near-equatorial Keplerian disk, $\theta_{0}=1.51$ case, we can see that the charged particle trajectories remain close to the equatorial plane, demonstrating epicyclic oscillations in the radial and vertical directions – this kind of regular motion around magnetized black holes has been discussed in detail for both Schwarzschild Kološ et al. (2015) and Kerr Tursunov et al. (2016) black holes and can be related to the high-frequency quasiperiodic oscillations observed in microquasars Stuchlík et al. (2013); Kološ et al. (2017). For Keplerian disks that are slightly off-equatorial, $\theta_{0}=1.21$ case, we observe after ionization a transition to the chaotic particle motion demonstrating slight tendency to wind up and down along the magnetic field lines – the disk seems to be deformed into a toroidal structure. For Keplerian disks with middle inclination, $\theta_{0}=0.91$ case, the particles of the ionized disk demonstrate strongly chaotic motion with extension comparable in both radial and vertical dimensions, causing that the transformed disk resembles a quasi-spherical structure of orbiting charged particles. If the initial inclination of the Keplerian disk is slightly lowered, $\theta_{0}=0.61$ case, we observe a clear tendency of the charged particle motion to follow the direction of the magnetic field lines, but the transformed disk still resembles a quasi-spherical structure. On the other hand, for small and very small initial inclination angles of the Keplerian disk, $\theta_{0}=0.31,0.01$ cases, we observe fall of the charged particles into the black hole (direct fall in the case of very small inclination angle) – in such situations the internal region of the Keplerian disk becomes unstable after ionization, being completely destroyed and captured by the black hole.

IV.2 Non-linear methods determining chaos-regularity transition

Using various non-linear methods for determination of chaos and its measure, we realize detailed study of the stability of the Keplerian disks orbiting a magnetized Schwarzschild black hole, and the measure of chaos in the motion of particles of the ionized disk, in dependence on the initial inclination angle between the disk plane and the magnetic field lines direction $\theta_{0}$, and the position of the particles at the Keplerian disk $r_{0}$. We shall study the cases presented above, for the characteristic values of the magnetic parameter determining the intensity of the electromagnetic interaction between the ionized particles and the external magnetic field ${\cal{B}}$. We first study the case of the magnetic parameter ${\cal{B}}=0.1$, testing the efficiency of a large variety of the non-linear methods. The results are presented in Figs 5, 6. We test the following non-linear methods: Box-Counting, Correlation Dimension, Lyapunov Exponent, Reccurence Quantification Analysis with variants RR (recurrence rate), DET (determinism – predictability of the system), LL (diagonal line length), ENTR (Shannon entropy – probability distribution of the line lengths), and Machine Learning with variants Random Forest, Neural Network, Linear Regression, Nearest Neighbours, Gradient Boosted Trees. (All of these non-linear methods of detecting chaos are presented and tested by a simple logistic map function in the Appendix.) Finally, we study in detail the cases of ${\cal{B}}=\pm 0.001,\pm 0.01,\pm 0.1,\pm 1$, comparing the application of the Correlation Dimension method and the Machine Learning method with the Random Forest algorithm.

IV.2.1 Description of the numerical methods

For numerical trajectory calculations, the Wolfram Mathematica software was used with fixed step size fourth order explicit Runge-Kutta method. The size of the step has been set to $10^{-2}$ and the time of integration to $10^{4}$, then every 100-th point was taken; therefore, all trajectories which have not fallen into the black hole are of the length of $10^{4}$ points. The fixed step method has been used to ensure that the non-linear methods for calculation of the chaos-measure of every pair of the initial conditions take the input of the same length. All the trajectories that fall into the black hole were not taken as input for calculation of the chaos-measure, and are denoted by white colour in the final figures 5, 6, 7, and 8. Initial conditions for the charged particles from the ionized Keplerian disk radius $r_{0}$ have been taken from interval $r_{0}\in[6,12]$, while inclination between the disk plane and the magnetic field lines $\theta_{0}$ has been taken from interval $\theta_{0}\in[0,\pi/2]$. Both these intervals have been evenly distributed among $10^{2}$ points; all their combinations lead to $10^{4}$ couples of conditions for particle trajectories, plotted for every figure from Figs 5, 6, 7, and 8. All the computations have been realized on DELL 7577 with processor Intel® Core i7-7700HQ and 16GB 2400MHz DDR4 SDRAM.

As can be seen from Fig. 5, the ionized particle trajectories are mostly regular when the Keplerian disk is near the equatorial plane, i.e., the inclination angle is of the value $\theta\sim\pi/2$ and the magnetic field is almost perpendicular to the disk plane. As the inclination angle is gradually decreased, we see gradual increase of the chaotic character of the charged particle motion. Near the inclination angle $\theta\sim\leavevmode\nobreak\ 1.1$, we enter a fractal like zone, where the strongly chaotic trajectories are located in close vicinity to those not so much chaotic (or regular). Finally, when the inclination angle is small enough (say $\theta\sim\leavevmode\nobreak\ 0.7$ for $r=7$) all the trajectories will eventually end being captured by the black hole. This ”capture inclination angle” depends on the radial position inside the disk - trajectories with initial position closer to the black hole get captured for larger inclination angle. When the magnetic field will be parallel (almost parallel) to the ionized Keplerian disk ($\theta\sim\leavevmode\nobreak\ 0$), all the trajectories will be captured by BH. But one can argue, that such configuration is not so realistic, since the currents inside the ionized accretion disk always produce some internal magnetic field with component parallel to the disk plane.

Except the fate of the captured parts of the ionized Keplerian disks, we have to follow the fate of the rest of such disks whose charged particles remain bound in vicinity of the black hole, following generally chaotic trajectories with possible exceptions of special regular motions. The measure of the chaotic character of the motion, and existence of possible ”islands of regularity” can be treated by the non-linear methods of determination of deterministic chaos.

All the non-linear methods we have used are revealing similar structures at the Figs. 6, 7, and 8. The colour scale (colour palette) is different for each method, because each of the non-linear methods denotes trajectories by different scale of the chaos-measure estimation – see section A. Therefore, the sensitivity of different non-linear methods to the chaotic behaviour is different, so it seems to be quite complex task to use one common colour scale for all methods. Moreover, when plotting various chaos estimations in one figure, the presence of outreach values is inconvenient – if there is some value far away from the others, then it is problematic to cover the differences between closer values by visible colour differences. At Figs. 5, 6, 7, and 8, we should focus on structures produced by different colours - such structures will indicate transition between regular and chaotic motion.

Box-Counting method A.1.1 (see Appendix) is clearly the fastest method and that is a big advantage while studying large datasets. The reason is the simplicity of the algorithm which does not need to make many numerical operations, and therefore does not need much computing power. We also have used only one-dimensional input in the sense of non-embedding dimension (for the Box-Counting, but for the Correlation Dimension we did use it), which non-linearly enlarges the computing, while its application did not seem to improve our final figures. The Correlation Dimension method A.1.2 is, along with the Lyapunov Exponent A.1.3 method, the most time dependent, when we consider the fact that from Recurrence Quantification Analysis (RQA) method A.1.4 computation we get four figures. This property varies also with the density of points in time series. For the Lyapunov Exponent method, together with ML methods, the colour range seems to be highly continuous, while the Correlation Dimension method shows some structure details, similarly to the RQA methods. RQA demonstrates similar structure for all the tools (RR - recurrence rate, DET - determinism, LL - line length, ENTR - entropy). Time needed for the RR variant could be shorter than for the DET variant, which also needs lower computing time than the ENTR and LL variants. This is caused by the necessary numerical operations needed for computing each of them – while more complicated numerical operations for the ENTR and LL variants are already done, the RR and DET variants can be derived from them, but not vice versa. As already mentioned, sometimes in the results of these non-linear methods a big gap between the most and the least chaotic trajectory occurs, implying that the trajectories in some interval of the chaos-measure are not distinguished by distinct colours, and the chaotic structure is not observable with sufficient precision. This case can be handled by adjusting of the outreach values, or by defining special rules for plotting. The implementation of the ML models A.1.5 is also very fast, as the time needed to build up a model has not been counted in Figs 6, 7, and 8. When considering time to build up such a model, several factors have to be taken into account, like choice of the programming language and its library, concrete algorithm and its parameters, data preprocessing, model validation, etc. We can observe that different ML algorithms are producing structures which are not so similar, as in the case of the previous methods (with the same training set input). This is affected by the fact that behind each ML algorithm is different principle of achieving the goal (data classification or regression). It has been proved that some of them are relevant mostly for some specific tasks – by this not only the specific application for particular tasks like, e.g., image recognition is meant, but also the usage according to available data, the model generalization (bias - variance trade-off), memory consumption, or the possibility of additional customization (neural network).

IV.2.2 Evolution of ionized Keplerian accretion disk

We have assumed simple ionization scenario for particles from a neutral Keplerian disk, where the neutral particles will split into two oppositely charged particles, and the newly charged particles start to feel influence of an inclined external uniform magnetic field. In the most common process of such ionization, the produced charged particles have quite different masses - proton for example is by three orders more massive than electron, and hence takes almost whole initial momentum of the former neutral particle, see (17). On the other hand, the electron, with specific charge $q/m$ by three orders higher than that of the proton, experiences the influence of the magnetic field due to the Lorentz force by three orders of magnitude stronger than the proton, i.e., for electrons the magnetic parameter ${\cal{B}}$ will be higher by three orders in comparison with those related to protons – see eq. (10). In order to consider all possible inclinations of the disk to the magnetic field lines, we restrict our attention to the magnetized non-rotating Schwarzschild black holes when particle energy after ionization process remains constant. In the present paper we have studied the motion of the heavier particles (protons, ions) after the ionization, such particles at bounded orbits around black hole will mix with neutral particles of the Keplerian disk and form new deformed lightly charged accretion disk, or the inner region of the disk could be ionized completely. The electrons are assumed to escape from the black hole neighbourhood in the form of wind, or form a corona around the accretion disk.

In Fig. 7 all possible inclinations of magnetic field vector to the accretion disk plane $\theta\in(0,\pi/2)$ are considered for various characteristic values of the interaction with the magnetic field ${\cal{B}}$. As shown in Fig. 3, the accretion disk destiny strongly depends on the magnetic field strength.

For very small magnetic field magnitudes $|{\cal{B}}|=0.001$ the motion of ionized particles is almost regular for any inclination and surprisingly it is more regular when the inclination of the magnetic field is small. Also almost all particle trajectories will remain in disk, only the inner edges of the disk in ${\cal{B}}=-0.001$ case are captured by the Schwarzschild black hole.

For larger magnetic field magnitudes, the situation is completely different and more complicated - the ionized accretion disk has completely different evolution in the ${\cal{B}}>1$ and ${\cal{B}}<1$ cases. For ${\cal{B}}=-0.01$ (and ${\cal{B}}=-0.1$) all (and almost all) ionized trajectories are captured by the black hole, and the disk (inner parts of the disk) will be destroyed. For ${\cal{B}}=0.01$, only the disks with small and large inclinations are preserved, for large inclinations some thick disk of ionized particles will be formed, for small inclinations the ionized trajectories are forming some sphere around the central black hole. The case ${\cal{B}}=0.1$ has been already shown in Figs. 5 and 6 - disks with small inclination are destroyed, mildly inclined disks are formed by particles following strongly chaotic trajectories creating a (quasi-)spherical configuration, and charged particles from disks with almost perpendicular inclination of magnetic field will remain near the disk plane.

If the magnetic field parameter is really large, $|{\cal{B}}|\geq 1$, as it could be for elementary particles with large specific charge $q/m$, the Lorentz force is leading force in the system. The charged particles are winding up and down in vertical direction along the magnetic field lines, while slowly moving along the central black hole. Both cases ${\cal{B}}=-1$ and ${\cal{B}}=1$ are giving the same resulting trajectories, differing only in the direction of slow motion around the central black hole – see Fig. 3.

Our numerical study thus determines the relation between regions of regular and chaotic motion of the charged particles of the ionized Keplerian disk. Moreover, also the critical ”capture inclination angle” is determined in dependence on the initial position in the Keplerian disk for fixed values of magnetic field parameter ${\cal{B}}$.

In order to illustrate in detail the role of the combined effect of the magnetic field parameter ${\cal{B}}$ and the disk inclination angle, we study distribution of the capture by the black hole, and the relation of the regular and chaotic motion, for charged particles located initially at the ISCO. The results are given in Fig. 8 for ${\cal{B}}\in[-1,1]$ and inclination $\theta_{0}\in[0,\pi/2]$. All the non-linear methods of chaos determination give again the same results clearly reflecting the asymmetry of the capturing process in relation to the sign of the magnetic parameter.

IV.3 Astrophysical relevance

We have studied the influence of an external uniform magnetic field showing that even small ${\cal{B}}$ parameter has significant influence on the thin Keplerian accretion disk, if the particles of the disk become ionized, as shown in Fig. 3. Such qualitative description should be followed by qualitative analysis to clear up if all these effect are astrophysically relevant.

In this article we assume the magnetic field to be uniform, but real magnetic fields around microquasar or supermassive black holes, and their accretion disks, are far away from being completely regular and uniform. The Wald uniform magnetic field solution is used as a useful approximation describing properly at least the magnetic field magnitude. Moreover, the magnetic parameter ${\cal{B}}$ contains, together with the field strength, also the specific charge of the ionized test particles, see Eq. (10). This implies that in order to make proper estimation of the magnetic field magnitude, one needs to identify first the type of matter inside accretion disk.

In our approach, the ”charged particle” can represent matter ranging from electron to some charged inhomogeneity orbiting in the innermost region of the accretion disk. The specific charges $q/m$ for any of such structures will then range from the electron maximum to zero. Recalling the physical constants in the dimensionless magnetic parameter as ${\cal{B}}=|q|BGM/(2mc^{4})$, we get from eq. (10) the magnetic field strength in Gauss

[TABLE]

where the quantities are given in CGS units, see Tab. 1.

The range of the magnetic field magnitude for ${\cal{B}}=0.01$ magnetic parameter and stellar $M={10\leavevmode\nobreak\ {M}_{\odot}}$, or supermassive $M={10^{8}\leavevmode\nobreak\ {M}_{\odot}}$, BH and various specific charges of ionized material from accretion disk is quite wide and hence it should not be complicated to find astrophysically relevant situation for presented model of ionized disk. Even weak magnetic field, like for example Galactic magnetic fields $10^{-5}\leavevmode\nobreak\ {\rm{G}}$, could have strong impact on ionized accretion disks. Another point to address is the difference of Lorentz force magnitude action on electron and proton. For the same magnitude of magnetic field $B$ in Gauss, the electron magnetic field parameter ${\cal{B}}_{\rm e}$ will be 1836 time bigger then the proton magnetic field parameter ${\cal{B}}_{\rm p}$, hence we can assume the protons should stay in accretion disk plane while electrons will escape disk plane along magnetic field lines in properly tuned magnetic field.

V Conclusions

When the inclination of the magnetic field to the neutral Keplerian disk is almost perpendicular, the ionized particle motion can be almost completely regular. As the magnetic field inclination to the field decreases, the charged particle motion becomes chaotic with chaos-measure increasing with decreasing inclination to middle angles, until some limiting value of the inclination after which the charged particle is captured by the black hole. The most chaotic behaviour of the ionized particles, and the originally Keplerian disk destruction can be expected for magnetic field strength $|{\cal{B}}|\in(0.01,1)$, weaker magnetic fields, with ${\cal{B}}\in(0,0.01)$, are not strong enough to destroy the disk. The very strong magnetic fields, with $|{\cal{B}}|\geq 1$, which is from the astrophysics points of view also the most relevant case for the elementary particles, will take complete control over the particle dynamics due to the dominating Lorentz force, and the charged particles mainly follow the magnetic field lines.

All presented non-linear methods have been used for chaos determination in the charged particle motion - all methods are giving similar results, but the CPU time consumption differs for different methods. This could be due to algorithm programming, but it is dominantly caused by the fact that different numerical methods have different computational complexity requiring different computational times. One of the most promising methods is the Machine Learning, which can determine distinction between regular/chaotic trajectory very fast. All the considered methods are giving similar results while measuring the chaos in the charged particle motion – the structures presented in Fig. 6 are similar. They have different scales for the chaos-measure and thus they also differ in the colour palette reflecting different ranges of the chaos measure given by different non-linear methods.

Variety of phenomena related to the combined gravo-magnetic effects on charged test particle motion around rotating Kerr black holes is much richer than for the Schwarzschild black holes, especially, the black hole rotation allows the ionized particles to escape to infinity along the direction of the external magnetic field Stuchlík and Kološ (2016). The study of charged particles ’kicked’ from the innermost stable circular orbit (ISCO) in the equatorial plane, and hence escaping to infinity, has been treated in Al Zahrani et al. (2013); Shiose et al. (2014); Al Zahrani (2014). Actually, for some large enough magnetic field parameter $|{\cal{B}}|\geq 1$, the charged particles must escape from the equatorial plane even with zero kick Stuchlík and Kološ (2016). The reason for such an effect is that the particle orbit is stable in the radial $r/x$ direction, being unstable in the vertical $\theta/z$ direction, and the particle is forced to escape in the vertical direction – see the second trajectory in Fig. 11 from Tursunov et al. (2016). The exploration of MPP effect in such configuration, and the possibility of the black hole rotational energy extraction due to the discharge of the Wald black hole charge is an interesting new phenomenon governing acceleration of the ultra-relativistic jets Tursunov et al. (2018a).

It will be also interesting to examine how the chaos of the charged particle motion will be modified, and how the ionized Keplerian disk will be influenced by the external magnetic field, when the charged particle radiation reaction Tursunov et al. (2018b) will be taken into account. As one can expect, the radiation reaction force will act as dumping in harmonic oscillator, forcing the chaotic particle trajectory to become more regular. As related to the observational data, we have to recall the argument that the signal from the accretion disks orbiting neutron stars with magnetic field having much larger magnitude than those related to black holes tends to be more chaotic in comparison with the signal from the accretion disks orbiting magnetized black holes Karak et al. (2010). We plan to develop both detailed theoretical models of observational effects related to the ionized accretion disks and more complex structures, and further testing of non-linear methods of measuring the level of chaotic motion applied in the more detailed treatment of theoretical data, and the related real observational data.

Acknowledgments

The authors acknowledge the Silesian University in Opava Grant No. SGS/14/2016. M.K. acknowledges the Czech Science Foundation Grant No. 16-03564Y, Z.S. acknowledges the Albert Einstein Center for Gravitation and Astrophysics supported by the Czech Science Foundation Grant No. 14-37086G.

Appendix A Methods of determination of chaotic behaviour in time series

The motion of charged particles in the field of a magnetized black hole demonstrates a mixture of regularity and deterministic chaos. Deterministic chaos is a hardly predictable and apparently random behaviour which can appear in dynamical systems, being detectable by non-linear methods only. We first present a short survey of the methods detecting the deterministic chaos, and transitions between the regular and chaotic behaviour of dynamical systems, and then we test these methods in the case of so called logistic map representing one of the simplest dynamical systems.

A.1 Methods detecting chaos

Dynamical systems, represented by time series of data, can demonstrate regular behaviour, completely random behaviour, or the so called deterministic chaos. The standard linear methods as the Fourier analysis are not relevant in distinguishing these three types of data, therefore, non-linear methods of treating such data samples are necessary. Such methods are able to give even a measure of chaos in the detected data. Here we give a survey of these methods.

A.1.1 Box-counting

Box-counting ($D_{0}$) or box dimension is one of the most widely used methods of estimations of fractal dimension. The calculation and empirical application of this method is simple as compared to other methods. We present here the general idea behind this method, for detailed description see Schroeder (1991). For a set $S$ in a Euclidean space $\mathbb{R}^{n}$, we define the Box-counting measure as

[TABLE]

where $N(\epsilon)$ is the number of boxes of side length $\epsilon$ required to cover the set. The dimension $D_{0}$ of the set $S$ is estimated by seeing how the logarithmic rate of $N(\epsilon)$ increases as ${\epsilon\to 0}$, or in other words, as we make the grid finer.

A.1.2 Correlation dimension

Very popular tool for detecting chaos in experimental data is calculation of the Correlation dimension ($D_{2}$). The general idea behind computing correlation dimension is to find out for some small $\epsilon$ the number of points $C(\epsilon)$ (correlation sum), which Euclidean distance is smaller than $\epsilon$. The definition of the Correlation dimension takes the form

[TABLE]

and one can compute the number $C$ for various values of $\epsilon$ and $D_{2}$ can be approximated again by fitting the logarithmic values.

It is worth to mention that $D_{2}$ and $D_{0}$ are subclasses of the general definition of $D_{q}$ determining the family of fractal dimensions Schroeder (1991). The fractal dimensions are defined by the relation

[TABLE]

where, $p_{k}$ denotes relative frequency with which the fractal points are falling inside the k-th cell.

For $q=0$ we obtain the Box-counting dimension, for $q\to 1$ we obtain the so called information dimension, the numerator of which is denoted as Shannon’s entropy, and for $q=2$ we obtain the Correlation dimension.

There is several approaches to the calculation procedure of the correlation dimension – let us mention the approximation of $C(\epsilon)$ by the algorithm of Grassberger and Procaccia Grassberger and Procaccia (1983)

[TABLE]

where $H$ is the Heaviside step function.

When counting correlation dimension (fractal dimension), one should not omit importance of the embedding dimension, and also the close connection to Takens’s theorem about reconstruction of state space from sequence of observations. The embedding dimension is created from series of length $N+m-1$ for some given $m$ series of $N$ vectors, where i-th component takes the form

[TABLE]

One of the purposes of the embedding dimension is to distinguish between chaotic and random time series. In the chaotic series, for increasing $m\in\mathbb{N}$, the fractal dimension stabilizes at some value $D<m$, while in the random series, the fractal dimension goes along with $m$ to infinity.

A.1.3 Lyapunov exponent

In short, a Lyapunov exponent is a number giving the measure of separation of trajectories that are initially infinitesimally close. Near trajectories in chaotic systems diverge exponentially, what corresponds to positive Lyapunov exponents. The amount of separation can differ for various directions of the initial separation vector. Because of this fact, there is a spectrum of the Lyapunov exponents, which corresponds to the phase space dimension.

For example, if we consider the logistic map $f(x)=rx(1-x),r\in[0,4]$ as a typical example of simple chaotic system, which represents one-dimensional non-linear difference equation, the maximal Lyapunov exponent can be calculated directly from the expression of the $f$ function from the formula

[TABLE]

However, we assume here handling with time series inputs only, forbidding application of such a formula. This approach leads us to use method, which is determining the Lyapunov exponents from time series, being described in Kantz (1994). The maximal Lyapunov exponent characterizes the spectrum, and therefore denotes amount of predictability for the considered dynamical system 111Note that the Lyapunov coefficients have a crucial role also in other branches of physics, e.g., they characterize instability of quasinormal modes of perturbation fields around black holes (Cardoso et al., 2009; Konoplya and Stuchlík, 2017; Toshmatov et al., 2018)..

The Lyapunov exponents can be calculated without knowledge of a model behind the dynamical system, being based only on the observed time series. We use the method based on the statistical properties of the divergence of the neighbouring trajectories, introduced by Kantz (1994) Kantz (1994). We denote our numerical estimation by $\lambda$, which is actually sophisticated way of estimating true Lyapunov exponent as the average of maximal effective exponents $\lambda_{\tau}$ along the trajectory. The algorithm is well explained in Kantz’s article, it is quite fast and powerful method. The result of this numerical approach applied to logistic map is very similar to the result of 30 as one can compare in Fig. 9. The observable differences appears in this case by non-chaotic or even constant time series, which are negligible for the purposes of our study.

A.1.4 Recurrence quantification analysis

The recurrence quantification analysis (RQA) is a widely used tool for investigation of the state space trajectories – it determines the number and duration of recurrences of a dynamical system. RQA was developed since 1992 by Zbilut, Webber Jr. and Marwan Zbilut and Webber (1992); Marwan (2008). Recurrence plots (RP) provide a graphical tool for observing periodicity of phase space trajectories and was introduced by Eckmann et al.in 1987 Eckmann et al. (1987). This observing is possible through visualization of a symmetrical square matrix, in which the elements correspond to times at which a state of a dynamical system recurs.

We can define a RP which measures recurrences of a trajectory $x_{i}\in R^{d}$ in phase space

[TABLE]

where $N$ is the number of measured points $x_{i}$, $\epsilon$ is a threshold distance, and $\|\cdot\|$ is a norm. From this equation we obtain the already mentioned symmetrical square matrix of zeroes and ones. When we represent the two repeating elements with different colours in a plot, we obtain the discussed RP. Threshold value parameter determines density of the RP plot.

For investigating of the chaotic trajectories we use the following RQA variant tools:

RR - The recurrence rate is simplest tool, measuring the density of the recurrence points in the RP – in another words, it divides the number of the recurrence points in the RP by the number of all elements in the matrix. RR tool reflects the chance that some state of the system will recur

[TABLE] 2. 2.

DET - Determinism is the rate of the recurrence points which build the diagonal lines. DET determines how predictable the system is being formally defined as

[TABLE]

where $P(l)$ denotes the frequency distribution of the lengths $l$ of the diagonal lines. 3. 3.

LL - Line length tool represents average diagonal line length, which is in relation with the time of predictability of the system. It reflects the average time for which any two parts of the considered trajectories are close – this time can be denoted as mean prediction time

[TABLE] 4. 4.

ENTR - Entropy (or the Shannon entropy) of the probability distribution of the diagonal line lengths p(l), which are reflecting complexity of the system deterministic structure. The definition reads

[TABLE]

where $p(l)$ is the probability that a diagonal line is exactly of the length $l$ – it can be estimated from the frequency distribution $P(l)$ with

$p(l)=\frac{P(l)}{\sum_{l=l_{min}}^{N}P(l)}$.

A.1.5 Machine learning

Machine learning is a very powerful tool, which finds application in many quite diverse fields – language translating algorithms, computer vision, beating best Go player in the world Silver et al. (2017a), or Chess programs with different architecture Silver et al. (2017b). The purpose of this paper is to briefly introduce its basic principles and application to the study the chaotic motion.

Roughly speaking, machine learning (ML) is a field of computer science, which is strongly connected with another fields of mathematics as optimization, statistics, linear algebra, etc. Its beginning goes to 1950’s and as many inventions in computer science, or in science generally, the ML was not invented by a single person, but we mention at least A. Samuel, who used for the first time the term ”Machine learning”. The ML is using algorithms on data samples to discover known or unknown patterns in data. This dividing of patterns leads to basic divisions of the ML, namely the supervised, semi-supervised and unsupervised learning. The wide range of applications announces the good ability of handling non-linear dynamical systems data.

Our intention of using the ML method is to decide whether a trajectory of a charged particle is chaotic or not. For this purpose we use the supervised ML, where we train various ML algorithms with examples from the non-linear methods we have already described. The training set consists overall of 4725 examples. Chaos is very complex phenomenon which is reflected by various another phenomena like sensitivity to initial conditions, fractal dimension of coordinates time series, attractors, or the predictability of the dynamical system. We include into our training data set the results from all the already mentioned methods (Box-counting, Correlation dimension, Lyapunov exponent, RQA - RR, DET, LL, ENTR ) with effort to use all the different properties of them. The appropriate number of trajectories has been assigned to evenly weighted sum of previously calculated chaos-measure estimations based on previously mentioned methods. Various types of the ML algorithms has been trained with these data, the most of them are presented in Fig. 6.

The ML algorithms considered in our study are the following; the description is exactly taken from the cited sources:

Decision Tree Wolfram—Alpha (2018a) – A decision tree is a flow chart–like structure in which each internal node represents a ”test” on a feature, each branch represents the outcome of the test, and each leaf represents a class or value distribution. 2. 2.

Neural Network Wolfram—Alpha (2018b) – A neural network consists of stacked layers, each performing a simple computation. Information is processed layer by layer from the input layer to the output layer. The neural network is trained to minimize a loss function on the training set using gradient descent. 3. 3.

Linear Regression Wolfram—Alpha (2018c) – The linear regression predicts the numerical output $y$ using a linear combination of numerical features. The conditional probability $P(y|x)$ is modeled according to

[TABLE]

with $f(\theta,x)=x\cdot\theta$.

The estimation of the parameter vector $\theta$ is done by minimizing the loss function

[TABLE]

where $m$ is the number of examples and $n$ is the number of numerical features. 4. 4.

Random Forest Wolfram—Alpha (2018d) – Random forest is an ensemble learning method for classification and regression that operates by constructing a multitude of decision trees. The forest prediction is obtained by taking the most common class or the mean-value tree predictions. Each decision tree is trained on a random subset of the training set and only uses a random subset of the features (bootstrap aggregating algorithm). 5. 5.

Gaussian Process Wolfram—Alpha (2018e) – The ”GaussianProcess” method assumes that the function to be modeled has been generated from a Gaussian process. The Gaussian process is defined by its covariance function (also called kernel). In the training phase, the method will estimate the parameters of this covariance function. The Gaussian process is then conditioned on the training data and used to infer the value of a new example using a Bayesian inference. 6. 6.

Nearest Neighbours Wolfram—Alpha (2018f) – Nearest neighbors is a type of instance-based learning. In its simplest form, it picks the commonest class or averages the values among the k nearest neighbors. 7. 7.

Gradient Boosted Trees Wolfram—Alpha (2018g) – Gradient boosting is a machine learning technique for regression and classification problems that produces a prediction model in the form of an ensemble of trees. Trees are trained sequentially with the goal of compensating the weaknesses of previous trees. The current implementation uses the LightGBM framework in the back end.

In the final stages of our study, represented by Fig. 7, we decided to compare the results of the Correlation dimension method with the efficient ML - Random Forest method, applied to the sample of $10^{4}$ pairs of initial conditions for various magnetic field parameters values.

A.2 Test of the non-linear methods detecting chaos

Classical example of a simple non-linear dynamical system is the so called logistic map, given by the quadratic recurrence equation

[TABLE]

If the initial value of the variable is given, $x_{0}\in(0,1)$, the logistic map (38) generates sequences of real numbers $x_{n}\in(0,1)$ in dependence on the map parameter $r$. The behaviour of such sequence $x_{n}$ strongly depends on logistic map parameter that is considered in the interval $r\in[2,4]$. Roughly speaking, the behaviour on the interval $r\in(2,r_{0})$ is regular (rather predictable) and on the interval $r\in[r_{0},4]$ it is chaotic (hardly predictable or rather unpredictable), with some occasional ”islands of regularity”. The transition between regular and chaotic behaviour happens for logistic map parameter $r=r_{0}\approx\leavevmode\nobreak\ 3.56995$. Bifurcation diagram for the logistic map, with asymptotically approached values of the sequence, is shown in Fig. 9 by gray background points.

For chaos determination and chaotic behaviour description in sequences of numbers, the following methods have been tested: Box-Counting (section A.1.1), Correlation Dimension (A.1.2), Lyapunov Exponent (A.1.3), Recurrence Quantification Analysis (RQA) with variants (A.1.4), and Machine Learning with various algorithms (A.1.5). All of the algorithms we use in this article are capable to be applied for one dimensional sequences of real numbers (time series) taken as an input. From the point of observation, it is important to distinguish between chaotic sequences ruled by some (unknown) laws, and random sequences, where the numbers are chosen randomly. The theoretical boundary of distinguishing between chaos and random sequences is sequence length Ott (1993). For the non-linear dynamical system having many degrees of freedom, short sequences are not sufficient for any decision.

To clear up how the non-linear methods work, we use four representative sequences of length $10^{4}$

[TABLE]

where $x_{n}^{(1)}$ is a regular sequence, $x_{n}^{(2)}$ is a weakly chaotic sequence generated by the logistic map with $r=3.6$, $x_{n}^{(3)}$ is a strongly chaotic sequence generated by the logistic map with $r=4$, and $x_{n}^{(4)}$ is a sequence of pseudo-random distribution of numbers. The results are presented in Tab. 2.

We have tested in detail all the introduced non-linear methods of chaos determination on sequences of numbers generated by the logistic map for various values of the parameter $r$ – for results see Fig. 9. The important common sign of all these methods is that they are able to detect higher chaos-measure when divergence of the trajectories in the bifurcation point of the diagram at $r=3$ occurs, i.e., when there is more than one fixed point. However, this is not undeniable true for all the cases, what leads to different estimations of the chaos-measure for different non-linear methods. The Box-counting and Lyapunov exponent methods show approximately linear behaviour when moving to the point, where the divergence of the trajectories begins – while moving to this point, the chaos-measure is increasing, and it is decreasing while leaving it. The Correlation dimension method demonstrates little bit different behaviour, where non-chaotic region is denoted by values close to zero, and the level of chaos starts to rapidly grow only in small distance from the region of divergence of trajectories in the bifurcation point – this type of the behaviour is also demonstrated by the RQA tools, however, in a reversed way. The determination of chaotic behaviour for $r\geq 3.56995$, when all the non-linear methods are showing the highest values of the chaos-measure, can be considered as highly satisfactory. However, as in this region also the ”islands of regularity” are present, not all of the chaos-measure estimations should be high there.

As expected, different ML algorithms produce different results. The Random Forest algorithm, used also in (Huppenkothen et al., 2017) seems to produce quite reliable results – our estimation of the chaos-measure is getting higher as the $r$ parameter grows, with exceptions that could be explained by the fact there are also ”islands of regularity” in the chaotic region. This is the reason why we decided to include the ML Random Forest algorithm to the figure 7; but we have to stress that the others ML algorithms seem to produce interesting and relevant results as well.

The time spent for the calculations of every single figure from Fig. 9 is presented in Tab. 3. From here we can see that the Box-counting method is the fastest, but when looking on the results, the other non-linear methods seem to be more precise. The time spent on the calculations also strongly depends on the sequence length - such dependence is non-linear in the case of almost all the methods, but we should point out that Box-counting and ML seem to behave generally in rather linear way while working with larger data sets (however, this general conclusion is not fully confirmed by the results presented in Tab 3).

For most of the non-linear chaos detecting methods we have programmed and tested several variants according to the algorithm architecture, and settings of the adjustable parameters, for the code see Pánis (2019). The goal was to obtain relevant results and acceptable computing time.