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

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.

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 [9][10][11].
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 [12] 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 [8,9,[13][14][15][16][17][18][19][20][21][22][23][24][25][26]. 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: (1) regular oscillatory motion possibly reflecting the high-frequency X-ray quasiperiodic oscillations observed in microquasars [11], (2) destruction of the ionized region of the disk due to the radial infall into the black hole, (3) 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 [8]), see Fig. 1.
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 B, initial position of the particle given by the radius r , and the inclination angle of the disk to the magnetic field lines θ . 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 B, we determine the distribution of the regular and chaotic states in the plane θ − 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 nonlinear 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 [27]. The non-linear test of the observed microquasar timing data has been already applied in different context [28][29][30][31].
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.

Charged particle motion around magnetized Schwarzschild black holes
The line element of the Schwarzschild black hole spacetime with mass M reads (1) where the so called lapse function f (r ) takes the form 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 [32] The Hamiltonian governing the test (non-radiating) charged particle motion can be written in the form [32] where the kinetic four-momentum p μ = mu μ is related to the generalized (canonical) four-momentum π μ by the relation that satisfy the Hamilton equations presented in the form Eur. Phys. J. C (2019) 79:479

Fig. 2
Schematic figure of accretion disk around black hole immersed into uniform magnetic field with magnetic field vector B aligned with z axis. Neutral accretion disk plane is inclined to the z axis by θ 0 ∈ (0, π/2) angle and disk inclination to equatorial plane is π/2 − θ 0 . The accretion disk is created by neutral particles orbiting the black hole along circular orbits with various radii r 0 . For spherically symmetric Schwarzschild BHs, the presented set up of aligned magnetic field with z axis and disk inclined to equatorial plane is equivalent to the disk in the equatorial plane with the magnetic field being inclined to the z axis The affine parameter ζ is related to the proper time τ of the particle by the relation ζ = τ/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 The dynamical equations for the charged particle motion in the Cartesian coordinates can be found using the coordinate transformations x = r cos(φ) sin(θ ), y = r sin(φ) sin(θ ), z = r cos(θ ). (9) Introducing for convenience the specific constants of the motion, energy E, axial angular momentum L, and the magnetic parameter B governing the intensity of the electromagnetic interaction, by the relations [9,23] one can rewrite the Hamiltonian (4) in the form where V eff (r, θ; L, B) denotes the effective potential governing the turning points of the radial and latitudinal motion, given by the relation The terms in the parentheses corresponds to the central force potential given by the specific angular momentum L, and electromagnetic potential energy specified by the magnetic parameter B.
The effective potential (12) demonstrates clear symmetry (L, B) ↔ (−L, −B) that allows to distinguish in the following only two situations corresponding to the minus configuration (L > 0, 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 (L > 0, 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, 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, B < 0, corresponds to the vector of the magnetic field B pointing downwards, while plus configuration, B > 0, corresponds to the vector of the magnetic field B pointing upwards the z-axis [23].
The charged particle motion is limited by the energetic boundaries given by The axial symmetry of the background of the combined gravitational and magnetic fields implies independence of the effective potential V eff from the coordinate φ which allows us to examine V eff (r, θ) as a 2D function of the spherical coordinates r, θ, 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 [23] -here we are directly applying these results.

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) [12], the equatorial plane of the Thin Keplerian accretion disk, created by neutral test particles following circular geodesics, and its evolution (thickening/destruction) when the influence of the magnetic field is switched on. The accretion disk particles following initially circular orbits is located in the central plane with inclination θ 0 . = 1.37 to the equatorial plane, while the magnetic field lines are everywhere aligned with the z axis (vertical direction). The uncharged test particles following the innermost stable circular geodesic (ISCO -depicted by the dashed circle) of the Schwarzschild spacetime represent the inner edge of the Keplerian accretion disk. If the disk will remain neutral or the magnetic field is missing (B = 0 case, middle figure) all the orbits will remain in their circular shape and we see just inclined razor thin disk. If a slightly strong electromagnetic interaction is switched-on (B = ± 0.001 cases, middle row), the charged particles forming the disk that is originally almost perpendicu-lar to the magnetic field lines start to follow epicyclic oscillations around the circular orbit in both radial and latitudinal directions; the accretion disk becomes to be slightly thick. If larger magnetic field is switched-on (B = ± 0.01, ± 0.1 cases), the charged particle motion becomes quite chaotic and the accretion disk is destroyed or transformed into thick toroidal structure. The complete destruction of the Keplerian disk can be seen in the B = − 0.01 case, when all the particles are captured by the black hole. If the magnetic parameter of the field that is switched-on is large (|B| ≥ 1 cases), the Lorentz force dominates the particle motion. The charged particles are spiralling up and down along the magnetic field lines, while slowly moving around the black hole in the clockwise (B > 0) or the counter-clockwise (B < 0) direction. The thin Keplerian disk has been destroyed and transformed into some special thick toroidal structure spacetime has to be perpendicular, i.e., the x − y plane -see In the initial state we consider a Keplerian accretion disk orbiting the magnetized Schwarzschild black hole that con-sists 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 θ 0 , giving one of the initial conditions for the motion of particles of the ionized Keplerian disk, see Fig. 2.
Eur. Phys. J. C (2019) 79:479 Fig. 4 The evolution of thin Keplerian accretion disk with various initial inclinations to the magnetic field lines is demonstrated for the magnetic parameter B = 0.1. The accretion disk coexists of many circular orbits initially located in central plane given by the inclination angle θ 0 related to the z axis. Neutral particle ISCO represents the inner edge of the Keplerian disk and is depicted by the dashed circle. For small inclinations of magnetic field lines to the Keplerian disk, the capture by the black hole is the only option for charged particles after ioniza-tion; the accretion disk is completely destroyed in this case. For middle inclinations, the ionized particles enter quite chaotic behaviour, moving up and down along the magnetic field lines -in such case the inner region of the thin disk is transformed into a quasi-spherical structure. For Keplerian disk almost perpendicular to the magnetic field lines, the ionized particles follow regular epicyclic trajectories in vicinity of the equatorial plane, increasing slightly the disk thickness

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) [33,34], 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 and the law of conservation of the canonical momenta (5) takes due to the charge conservation the form In many realistic scenarios, like neutron β decay or neutral atom ionization, one of the created charged particles is much more massive then the other one, m 2 /m 3 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 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. Page 7 of 21 479 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 [8,25], 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.

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 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 At the moment of ionization, the neutral or charged test particle is considered to be located on inclined circular orbit with initial position x α and initial four-velocity u α The specific angular momentum L and the specific energy E of the neutral test particles on the inclined circular orbits of the Keplerian accretion disk are given by the relations [32] 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 L and specific energy E in the form The magnetic field has in the Schwarzschild metric only one non-vanishing component, A φ (3), hence only the (canonical) specific angular momentum L is changed due to the ionization, while the particle specific energy E remains constant. The energy of the neutral particles on the circular geodesic orbits, given by (22), is always E ≤ 1, but energy 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 [8].

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, θ 0 ∼ π/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.

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 [35]. We thus first discuss properties of the ionized disks orbiting the magnetized Schwarzschild black hole with large enough initial inclination angle, θ 0 ∼ π/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 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 (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 (B = ± 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 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 B = 0.001. When a larger magnetic field is switched-on (B = ± 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 B = − 0.01 case, when all the orbits are captured by black hole. For even larger magnetic field (B = ± 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 B = − 0.1 when the structure of orbiting particles resembles a sphere. When the magnetic field parameter is large (B = ± 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 (B > 0) or counter-clockwise (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, |B| > 1, similar behaviour is demonstrated by the ionized Keplerian disk.
The influence of the initial inclination angle θ 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 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, θ 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 [23] and Kerr [24] black holes and can be related to the high-frequency quasiperiodic oscillations observed in microquasars [11,36]. For Keplerian disks that are slightly off-equatorial, θ 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, θ 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, θ 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, θ 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.

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 θ 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 B. We first study the case of the magnetic parameter B = 0.  Determination of chaotic behaviour of ionized particle motion. The charged particle motion is represented by time series created from the time dependence of its radial component r (τ ). The equi-distance sequence of points r (τ i ) is used to determine the measure of the chaos in the ionized particle motion by applying the non-linear method of Correlation dimension. Every point in this figure represents one charged particle trajectory for which the chaos-measure is determined. The ionized particles are represented by the initial disk radius r 0 , given on the horizontal axis, and the initial disk inclination angle θ 0 , given on the vertical axis; the magnetic parameter has been chosen to be B = 0.1. We consider the range of disk radii r 0 ∈ [6,12] corresponding to the inner region, and whole the range of the inclination angles θ 0 ∈ [0, π/2]. The chaos measure of the particle trajectory is represented by the colour of the point, white area represents the particle trajectories captured by the black hole. The colour scale related to the quantity giving the measure of chaos is shown on the bar -increasing number means increasing chaos-measure. The gray colour thus means the upper degree of chaos.
Since the chaos-measure quantities are defined in different ways for different non-linear methods, and their ranges are also different, in the following figures the colour mode is chosen in a way that tends to normalize the results obtained by different methods. In the present and all the following figures an initial Keplerian disk with given inclination angle θ 0 is determined by the horizontal slice with constant inclination θ = θ 0 ods of detecting chaos are presented and tested by a simple logistic map function in the Appendix.) Finally, we study in detail the cases of B = ± 0.001, ± 0.01, ± 0.1, ± 1, comparing the application of the Correlation Dimension method and the Machine Learning method with the Random Forest algorithm.

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 100th 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 Figs. 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 ∈ [6,12], while inclination between the disk plane and the magnetic field lines θ 0 has been taken from interval θ 0 ∈ [0, π/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 θ ∼ π/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 θ ∼ 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 θ ∼ 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 disktrajectories 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 (θ ∼ 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 Fig. 6 Comparison of the considered non-linear methods of chaos determination and measurement, as applied on the charged particle trajectories from initially Keplerian accretion disk as described in Fig. 5. Every figure corresponds to the disk (initial position of ionized particle) radii r 0 ∈ [6, 12] and inclination of magnetic field lines to the Keplerian disk plane θ 0 ∈[0, π/2]. Name of the used non-linear method, and total time for construction of the whole figure in format (hh/mm/ss), are also presented. The range of the quantities giving the chaos-measure in various methods are different, the colour scale is chosen accordingly, to reflect the maximum of the chaos-measure by a specific colour. The RQA methods have been reversed by transformation. We subtracted the calculated data from the corresponding maximal data value, for RQA -LL we also have lowered one outreach maximal value to the second maximum in order to make the structure visible 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 A1a (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 A1b is, along with the Lyapunov Exponent A1c method, the most time dependent, when we consider the fact that from Recurrence Quantification Analysis (RQA) method A1d 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 numeri-cal 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 chaosmeasure 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 A1e 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).

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 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 par-ticles 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 θ ∈ (0, π/2) are considered for various characteristic values of the interaction with the magnetic field B. As shown in Fig. 3, the accretion disk destiny strongly depends on the magnetic field strength.
For very small magnetic field magnitudes |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 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 B > 1 and B < 1 cases. For B = − 0.01 (and 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 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 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, |B| ≥ 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 B = − 1 and B = 1 are giving the same resulting trajectories, differing only in the direction of slow motion around the central black holesee 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 B.
In order to illustrate in detail the role of the combined effect of the magnetic field parameter B and the disk inclination angle, we study distribution of the capture by the black hole, and the relation of the Fig. 7 Evolution of an ionized Keplerian accretion disk in an external magnetic field characterized by typical magnetic parameters B = ± 0.001, ± 0.01, ± 0.1, ± 1. As in Figs. 5 and 6, the Keplerian disk (initial position of the ionized particle) radius r 0 ∈ [6,12] is given on the horizontal axis, while magnetic field inclination to disk θ 0 ∈ [0, π/2] is given on the vertical axis. Two non-linear methods of the Correlation Dimension (upper row), and the Machine Learning with Random Forest algorithm (lower row) has been applied and their results are compared. The presented results are similar for both methods, the Machine Learning method takes considerably less amount of computer time for the calculation regular and chaotic motion, for charged particles located initially at the ISCO. The results are given in Fig. 8 for B ∈ [−1, 1] and inclination θ 0 ∈ [0, π/2]. All the nonlinear 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.

Astrophysical relevance
We have studied the influence of an external uniform magnetic field showing that even small 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 qualita- Fig. 8 Similarly as in Fig. 6, comparison of the considered non-linear methods of chaos determination and measurement is applied on the charged particle trajectories from initially Keplerian accretion disk as described in Fig. 5, but instead of the radial coordinate, the magnetic field parameter B is used. Therefore every figure corresponds to the disk magnetic field parameter B ∈ [−1, 1] and inclination of magnetic field lines to the Keplerian disk plane θ 0 ∈ [0, π/2], while the initial position of ionized particle is fixed to the ISCO value r 0 = 6. For both parameters B and θ we used 101 equidistantly distributed values. Name of the used non-linear method, and total time for construction of the whole figure in format (hh/mm/ss), are also presented. The range of the quantities giving the chaos-measure in various methods are different, the colour scale is chosen accordingly, to reflect the maximum of the chaos-measure by a specific colour. As in Fig. 6, the RQA methods have been reversed by transformation. We subtracted the calculated data from the corresponding maximal data value, for RQA -RR, ENTR, LL we also have lowered some outreach maximal values in order to make the structure visible tive 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 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 B = |q|BG M/(2mc 4 ), we get from Eq. (10) the magnetic field strength in Gauss where the quantities are given in CGS units, see Table 1.
The range of the magnetic field magnitude for B = 0.01 magnetic parameter and stellar M = 10 M , or supermassive M = 10 8 M , 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 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 B e will be 1836 time bigger then the proton magnetic field parameter B 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.

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 |B| ∈ (0.01, 1), weaker magnetic fields, with B ∈ (0, 0.01), are not strong enough to destroy the disk. The very strong magnetic fields, with |B| ≥ 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 dif-ferent 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 motionthe 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 gravomagnetic 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 [8]. 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 [19,22,37]. Actually, for some large enough magnetic field parameter |B| ≥ 1, the charged particles must escape from the equatorial plane even with zero kick [8]. The reason for such an effect is that the particle orbit is stable in the radial r/x direction, being unstable in the vertical θ/z direction, and the particle is forced to escape in the vertical direction -see the second trajectory in Fig. 11 from [24]. 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.
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 [38] 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 [28]. 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. Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: It is not necessary to upload any data because they can be generated by using the codes in the citation [53].] Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

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.
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, nonlinear 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 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 [39]. For a set S in a Euclidean space R n , we define the Box-counting measure as where N ( ) is the number of boxes of side length required to cover the set. The dimension D 0 of the set S is estimated by seeing how the logarithmic rate of N ( ) increases as → 0, or in other words, as we make the grid finer.

b 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 the number of points C( ) (correlation sum), which Euclidean distance is smaller than . The definition of the Correlation dimension takes the form and one can compute the number C for various values of 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 [39]. The fractal dimensions are defined by the relation 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 → 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( ) by the algorithm of Grassberger and Procaccia [40] 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 One of the purposes of the embedding dimension is to distinguish between chaotic and random time series. In the chaotic series, for increasing m ∈ N, the fractal dimension stabilizes at some value D < m, while in the random series, the fractal dimension goes along with m to infinity.

c Lyapunov exponent
In short, a Lyapunov exponent is a number giving the measure of separation of trajectories that are initially infinitesimally 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 [41]. The maximal Lyapunov exponent characterizes the spectrum, and therefore denotes amount of predictability for the considered dynamical system. 1 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) [41]. We denote our numerical estimation by λ, which is actually sophisticated way of estimating true Lyapunov exponent as the average of maximal effective exponents λ τ 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 A6 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.

d Recurrence quantification analysis
The recurrence quantification analysis (RQA) is a widely used tool for investigation of the state space trajectoriesit determines the number and duration of recurrences of a dynamical system. RQA was developed since 1992 by Zbilut and Webber [42] and Marwan [43]. Recurrence plots (RP) provide a graphical tool for observing periodicity of phase space trajectories and was introduced by Eckmann et al. in 1987 [44]. 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 ∈ R d in phase space where N is the number of measured points x i , is a threshold distance, and · 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: 1. 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 1 Note 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 [56][57][58].
tool reflects the chance that some state of the system will recur 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 where P(l) denotes the frequency distribution of the lengths l of the diagonal lines. .

e 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 [45], or Chess programs with different architecture [46]. 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 nonlinear 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: 1. Decision Tree [47] -A decision tree is a flow chartlike 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. Neural Network [48] -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. Linear Regression [49] -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 The estimation of the parameter vector θ is done by minimizing the loss function where m is the number of examples and n is the number of numerical features. 4. Random Forest [50] -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. Gaussian Process [51] -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. Nearest Neighbours [52] -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. Gradient Boosted Trees [53] -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.
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 If the initial value of the variable is given, x 0 ∈ (0, 1), the logistic map (A14) generates sequences of real numbers x n ∈ (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 ∈ [2,4]. Roughly speaking, the behaviour on the interval r ∈ (2, r 0 ) is regular (rather predictable) and on the interval r ∈ [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  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 A1a), Correlation Dimension (A1b), Lyapunov Exponent (A1c), Recurrence Quantification Analysis (RQA) with variants (A1d), and Machine Learning with various algorithms (A1e). 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 [54]. 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 x ( where x (1) n is a regular sequence, x (2) n is a weakly chaotic sequence generated by the logistic map with r = 3.6, x (3) n is a strongly chaotic sequence generated by the logistic map with r = 4, and x (4) n is a sequence of pseudo-random distribution of numbers. The results are presented in Table 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 Table 3 Comparison of time in seconds required for different nonlinear methods applied on time series generated by logistic map. For parameter r varying from 2 to 4 with the step of 0.01 leads to 200 time series, with initial value x 0 = 0.1 we did set up the iterations to length 100, 1000 and 10,000 and compared the time needed for the calculation for given methods r -for results see Fig. 9. The important common sign of all these methods is that they are able to detect higher chaosmeasure 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 Boxcounting 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 ≥ 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 [31] 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 Fig. 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 Table 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 Table 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 [55]. The goal was to obtain relevant results and acceptable computing time.