The Lorenz system: hidden boundary of practical stability and the Lyapunov dimension

On the example of the famous Lorenz system, the difficulties and opportunities of reliable numerical analysis of chaotic dynamical systems are discussed in this article. For the Lorenz system, the boundaries of global stability are estimated and the difficulties of numerically studying the birth of self-excited and hidden attractors, caused by the loss of global stability, are discussed. The problem of reliable numerical computation of the finite-time Lyapunov dimension along the trajectories over large time intervals is discussed. Estimating the Lyapunov dimension of attractors via the Pyragas time-delayed feedback control technique and the Leonov method is demonstrated. Taking into account the problems of reliable numerical experiments in the context of the shadowing and hyperbolicity theories, experiments are carried out on small time intervals and for trajectories on a grid of initial points in the attractor’s basin of attraction.


Introduction
In 1963, meteorologist Edward Lorenz suggested an approximate mathematical model (the Lorenz system) for the Rayleigh-Bénard convection and discovered numerically a chaotic attractor in this model [76]. This discovery stimulated rapid development of the chaos theory, numerical methods for attractor investigation, and still attracts much attention of scientists from different fields 1 (see, e.g., [26,89,96,102,104]). The aim of this work is to discuss the estimates of the global stability and the difficulties of numerically studying the birth of self-excited and hidden attractors, caused by the loss of global stability.
Consider the classical Lorenz system ⎧ ⎪ ⎨ ⎪ ⎩ẋ with physically sound parameters σ, r > 0, b ∈ (0, 4] [76]. An important task in both qualitative and quantitative investigation of dynamical systems is the study of limiting behavior of a system after a transient process, i.e., the problem of revealing and analysis of all possible limiting oscillations. In a theoretical perspective, this problem relates to a generalization [44,63] of the second part of the celebrated Hilbert's 16th problem [34] on the number and mutual disposition of attractors and repellers in the chaotic multidimensional dynamical systems, and, in particular, their dependence on the degree of polynomials in the model (see also corresponding discussion, e.g., in [99,113]). Remark that 16th Hilbert problem is far from complete solution even for two-dimensional polynomial quadratic systems admitting only periodic attractors (see, e.g., [42,62]).
A limiting oscillation (or a set of oscillations) is called an attractor, if the initial data from its open neighborhood in the phase space (with the exception of a minor set of points) lead to a long-term behavior that approaches it, and its attracting set is called the basin of attraction (i.e., a set of initial data for which the trajectories tend to the attractor).
The examination of all limiting oscillations could be facilitated if for a system it is possible to prove the dissipativeness in the sense of Levinson (or D-property) ( [74], see also [67])-the existence of a bounded convex globally absorbing set containing all such oscillations. The Lorenz system (1) is dissipative in the sense of Levinson, and one can consider, e.g., the following absorbing set [60,61]: This implies that all solutions of (1) exist for t ∈ [0, +∞) and, thus, system (1) generates a dynamical system. For considered assumptions on parameters, if r < 1, then system (1) has a unique equilibrium S 0 = (0, 0, 0), and if r > 1, then system (1) has three equilibria: S 0 = (0, 0, 0) and two symmetric equilibria: S ± = (± (r − 1)b, ± (r − 1)b, r − 1).
For a system that is dissipative in the sense of Levinson, further study of existence of oscillations is usually connected with the stability of a stationary set. If any trajectory of a system tends to the stationary set, then the system is called globally stable. 2 In this definition, the stationary set can contain both stable (trivial oscillations) and unstable equilibrium states, i.e., the local stability of all equilibria is not required 3 . Various analytical methods, based on the Lyapunov ideas, have been developed (see survey, e.g., [47]) for the study of global stability and estimates of its boundaries in the space of parameters.
Farther, the birth of nontrivial oscillations in a system can be observed when crossing the global stability boundary in the space of parameters. In the simplest case, this is due to the loss of stability of the stationary set. Within this framework, it is naturally to classify oscillations in systems as self-excited or hidden [39,62,67,68]. Basin of attraction of a hidden oscillation in the phase space does not intersect with small neighborhoods of any equilibria, whereas a self-excited oscillation is excited from an unstable equilibrium. A self-excited oscillation is a nontrivial one if it does not approach the stationary states (i.e., ω-limit set of corresponding trajectory does not consist of equilibria).
In the general case, when the boundaries of global stability are violated, the birth of oscillations can occur due to local bifurcations in the vicinity of the stationary set (trivial boundary of the global stability) 4 or nonlocal bifurcations (hidden boundary of the global stability). If an attractor is born via such a nonlocal bifurcation of the loss of global stability, then the attractor is a hidden one since the basin of its attraction is separated from the locally attractive stationary set.
In practice, only problem of the birth of new nontrivial attractors is of interest because there are roundoff and truncation errors in numerical experiments and there is noise in physical experiments. That is, if the 2 We use the term "global stability" for simplicity of further discussion, while in the literature there are used different terms like "globally asymptotically stable" [108, p. 137], [31, p. 144], "gradient-like" [72, p. 2], [111, p. 56], "quasi-gradient-like" [72, p. 2], [111, p. 56] and others, reflecting the features of the stationary set and the convergence of trajectories to it. 3 This can be demonstrated on the example suggested by Vinograd [109] (see also [32, p. 191] Here, a planar system with a unique globally attractive equilibrium of the saddle type has a family of homoclinic orbits. Also, global stability in the presence of locally unstable equilibria is a typical case for systems describing pendulums, PLLs [11,21,46,69], and electric machines [64]. model has a global bounded convex absorbing set, then over time, the state of the system, observed experimentally, will be attracted to the local attractor contained in the absorbing set. This problem can be considered as a practical interpretation of the problem of determining the boundary of global stability.

Inner estimation: the global stability and trivial attractors
For the Lorenz system (1), using several approaches based on the construction of Lyapunov functions it is possible to obtain the following result.
Theorem 1 (criterion for the absence of self-excited and hidden oscillations) If for parameters of system (1) one of the following cases holds: then there are no nontrivial self-excited and hidden oscillations in the phase space of system (1), and any of its solution (x(t), y(t), z(t)) tends to the stationary set as t → +∞.
Proof Let us present the sketch of the proof (see [58,63,65,73]) and consider the following three cases.
(1) For r ≤ 1, system (1) has the unique equilibrium S 0 . Thus, the absence of self-excited oscillations follows from the Routh-Hurwitz criterion on local stability for the equilibrium S 0 . The absence of hidden oscillations can be obtained by Barbashin-Krasovskii theorem (see, e.g., [7,47]) and the Lyapunov function V (x, y, z) = 1 2 x 2 σ + y 2 + z 2 . For r ≤ 1, this function satisfies the following relation on the derivative with respect to system (1): thus, all conditions of Barbashin-Krasovskii theorem are satisfied implying global stability of the unique equilibrium S 0 . For r > 1, system (1) has three equilibria: the saddle equilibrium S 0 and the equilibria S ± , whose stability depends on the parameters r , σ , b. Thus, here the Barbashin-Krasovskii theorem is not applicable. In this case, system (1) may possess nontrivial self-excited oscillations with respect to unstable equilibria and also hidden oscillations. (2) For r > 1 and 2σ < b, the absence of nontrivial oscillations (and, thus, the global stability of the stationary set {S 0 , S ± }) can be demonstrated (see, e.g., [59]) by the LaSalle principal [54]. For that by the time and coordinate transformations t → τ, ψ : 0, and consider the following Lyapunov function: Note that the LaSalle principle requires the compactness of the set, where the Lyapunov function V is defined to show its boundedness from below. In our case, one can show that the inequality V (χ, ϑ, υ) > − 1 2 is valid for any (χ, ϑ, υ) ∈ R 3 . Since β < 0, function (5) satisfies the following relation on the derivative with respect to system (4): From the relationV (χ, ϑ, υ) = 0, it follows that the largest invariant set M ⊂ {(χ, ϑ, υ) ∈ R 3 V (χ, ϑ, υ) = 0} consists of the equilibrium points of system (4). Then, according to LaSalle principle any solution of system (4) (and thus system (1)) tends to an equilibrium state as τ → +∞. The level lines of Lyapunov function (5) surrounding the stationary set are presented in Fig. 1. For 2σ = b and β = 0, the third equation in (4) yields υ(t) = υ 0 exp(−αt) and thus any trajectory approaches an equilibrium as τ → +∞.
If 2σ ≥ b, then it is known [60,61] the following relation for the solutions of system (1): Using this relation, one can obtain the following estimation for all z ≥ x 2 2σ : The matrix Q has the following eigenvalues: for which relation Then, condition (6) takes the following form: , then using (12), one gets the following estimation: The later inequality If one considers the following coefficients (see, e.g., [58]): which are satisfied for b ∈ [2, 4] and r ≥ r 0 . Thus, for these values condition (6) takes the following form: (3). This completes the proof. The case corresponding to conditions (3) of Theorem 1, when all trajectories tend to the stationary set, however, not all equilibria of the stationary set are locally stable, is illustrated in Fig. 2.
Beyond the estimate (3) in Theorem 1, the analysis of global stability and the birth of nontrivial attractors can be performed numerically. It is further known that the separatrix of saddle S 0 can form a homoclinic loop from which unstable cycles can arise and violate the global stability (however, a set of measure zero does not affect the global attraction on a stationary set from a practical point of view). Using the Fishing principle [57,67,71] for the Lorenz system (1), it is possible to prove the following: For instance, for the classical values of parameters σ = 10, b = 8/3 of system (1), it is possible to find numerically the approximate value of such homoclinic bifurcation r h ≈ 13.926, when two symmetric homoclinic orbits appear forming a homoclinic butterfly (Fig. 3). A further increase in the parameter r leads to the birth of two saddle periodic orbits from each homoclinic orbit [93].

Outer estimation: the absence of trivial attractors
For systems with a global absorbing set and an unstable stationary set, the existence of self-excited attractors is obvious. From a computational point of view, this allows one to use a standard computational procedure, in which after a transient process, a trajectory, starting from a point of unstable manifold in a neighborhood of equilibrium, reaches a state of oscillation; therefore, one can easily identify it. System (1) possesses the absorbing set B (defined by Eq. (2)) and for σ > b + 1, r > r cr = σ ( σ +b+3 σ −b−1 ) all equilibria are unstable. Thus, in this case, system (1) has a nontrivial self-excited attractor. If we consider classical values of parameters σ = 10, b = 8/3, then for r > r cr , e.g., for r = 28, it is possible to observe the self-excited chaotic attractor with respect to all three equilibria S 0 , S ± (Fig. 4). This gives an outer estimation of the practical global stability.
In this manner, Lorenz had discovered numerically in the phase space of system (1) the existence of the celebrated chaotic Lorenz attractor from a vicinity of unstable zero equilibrium.

The boundary of practical stability and absence of nontrivial attractors
The presence of an absorbing set (Fig. 5) implies the existence of a globally attractor A glob , which contains all local self-excited and hidden attractors, and a stationary set. Thus, inside the set B it is possible to study numerically the presence of nontrivial self-excited and hidden attractors for parameters r , σ , b not satisfying conditions (3) of global stability, i.e., by fixing σ and b and by decreasing r from r cr . For σ = 10, b = 8/3, this gives us the following region r ∈ (r gs , r cr ), where r gs ≈ 4.64, r cr ≈ 24.74.
A nontrivial self-exited attractor can be observed numerically for 24.06 r < r cr ≈ 24.74 (see, e.g., [96]). In this case of nontrivial multistability, system (1) possesses a local chaotic attractor A which is selfexcited with respect to equilibrium S 0 and coexists with the trivial attractors S ± (Fig. 6).

Hidden attractor or hidden transient chaotic sets?
For the Lorenz system (1), it is still an open question [39, p. 14], whether for some parameters there exists a hidden chaotic attractor, i.e., whether it is possible by changing parameters to disconnect the basin of attraction from equilibria S 0 , S ± (e.g., for the parameters σ = 10, b = 8 3 : if r = 28, then attractor is connected with S 0 , S ± ; if r = 24.5, then attractor is connected with only S 0 ). The current results on the existence of the hidden attractors in the Lorenz system are the following. Recently reported hidden attractors in the Lorenz system with r < r cr and locally stable equilibria S ± turn out to be a transient chaotic set (a set in the phase space, which can persist for a long time, but after all collapses), but not a sustained hidden chaotic attractor [79,112].
Since in a numerical computation of a trajectory over a finite-time interval it is difficult to distinguish  (10)) containing self-excited chaotic attractor in the Lorenz system (1) with r = 28, σ = 10, b = 8 3 , when the stationary set {S 0 , S ± } is unstable a sustained chaos from a transient chaos [30,52], it is reasonable to give a similar classification for transient chaotic sets [14,18]. A transient chaotic set is called a hidden transient chaotic set if it does not involve and attract trajectories from a small neighborhood of equilibria; otherwise, it is called self-excited.
For an arbitrary system possessing a transient chaotic set, the time of transient process depends strongly on the choice of initial data in the phase space and also on the parameters of numerical solvers to integrate a trajectory (e.g., order of the method, step of integration, relative and absolute tolerances). This complicates the task of distinguishing a transient chaotic set from a sustained chaotic set (attractor) in numerical experiments. For the Lorenz system (1), suppose that σ = 10, b = 8 3 are fixed and r varies. Near the point r ≈ 24.06, it is possible to observe a long living transient chaotic set, which is hidden since its basin of attraction does not intersect with the small vicinities of equilibrium S 0 . For example, for r = 24 a hidden transient chaotic set can be visualized [112] from the initial point (2, 2, 2) (Fig. 7). In [79], hidden transient chaotic set was obtained in system (1) with r = 29, σ = 4, b = 2.
In our experiments, consider system (1) with r = 24. For a trajectory with a certain initial point, which is computed by a certain solver with specific parameters, we estimate the moment of the end of transient behavior as the moment of time when the trajectory falls into a small vicinity of one of the stable equilibria S ± . Using MATLAB's standard procedure ode45 with default parameters (relative tolerance 10 −3 , absolute tolerance 10 −6 ) for system (1) 10 7 ]. Remark that, if we consider the same initial points, but use MATLAB's procedure ode45 with relative tolerance 10 −6 , for all these initial points the chaotic transient behavior will last over a time interval of more than [0, 10 6 ], and corresponding transient chaotic sets will not collapse. As a conclusion from this numerical study of transient chaotic behavior, we may suggest to specify precisely numerical solver, its 3 by a trajectory that starts in the vicinity of the unstable equilibrium S 0 . This attractor coexists with stable equilibria S ± (trivial attractors) In [79,112], the study of existence of the hidden transient chaotic sets and hidden attractors in the Lorenz system was done by the numerical study of basins of attraction on some 2d cross sections. Next, we demonstrate an approach based on the extension of the parameter space and the numerical continuation method (see "Appendix A"). Consider the following generalization of the classical Lorenz system [76]: with parameters σ, r, b > 0 and parameter a = 0. Consider For the following linear transformation (see, e.g., [60]): system (1) is transformed to the Glukhovsky-Dolzhansky system [27]: where The Glukhovsky-Dolzhansky system describes the convective fluid motion inside a rotating ellipsoidal cavity.
If we set then after the linear transformation (see, e.g., [60]): with positive ν 1 , ν 2 , h, we obtain the Rabinovich system [83,87], describing the interaction of three resonantly coupled waves, two of which being parametrically excited: where Following [14], let us briefly describe the experiment, which allows to organize the transition from the hidden attractors in Glukhovsky-Dolzhansky and Rabinovich systems to the hidden transient chaotic set in the Lorenz system.
In this experiment for system (15), we consider three sets of parameters: P GD r = 346, a = 0.01, σ = 4, b = 1 (for the Glukhovsky-Dolzhansky system -GD), P L r = 24, a = 0, σ = 10, b = 8/3 (for the Lorenz system -L), and P R r = 24, a = −1/r − 0.01, σ = −ar, b = b cr + 0.14 (for the Rabinovich system -R). Here, in contrast to the case of classical Lorenz system considered in previous section, we change the parameters in such a way that hidden Glukhovsky-Dolzhansky and Rabinovich attractors are located not too close to the unstable zero equilibrium so as to avoid a situation that numerically integrated trajectory persists for a long time and then falls on an unstable manifold of the unstable zero equilibrium, then leaves the transient chaotic set, and finally tends to one of the stable equilibria.
These sets of parameters define three points, P GD , P L , and P R , in the 4D parameter space (r, a, σ, b). Consider two line segments, P GD → P L and P L → P R , defining two parts of the path in the continuation procedure. Choose the partition of the line segments into N st = 10 parts and define intermediate values of parameters as follows: . . , N st . Initial points for trajectories of system (15) that define hidden chaotic sets are presented in Table 1. At each iteration of the procedure, a chaotic set defined by the trajectory in the phase space of system (15) is computed using MATLAB's procedure ode45 with relative and absolute tolerances equal to 10 −8 . The last computed point of the trajectory at the previous step is used as the initial point for computation at the next step.
Thus, using numerical continuation procedure and starting from the hidden attractor in the Glukhovsky-Dolzhansky system (18), it is possible to organize a transition and localize numerically the hidden transient chaotic sets in both Lorenz (1) and Rabinovich (21) systems. One can ensure that in all three cases the trajectories (including separatrices), starting in small neighborhoods of the unstable equilibrium S 0 = (0, 0, 0), are not attracted by the computed chaotic set, but tend to the symmetric stable equilibria S ± .

Exact and finite-time Lyapunov dimension & Lyapunov exponents: global attractor, local attractors, and transient chaotic sets
In this section, we will demonstrate the difficulties in the reliable numerical computation of the finite-time Lyapunov exponents (FTLE) and finite-time Lyapunov dimension (FTLD) 5 , as well as distinguishing with their help an attractor from a transition set. Before proceeding to the description and discussion of these numerical experiments, let us remark that investigations of the numerical reliability of computergenerated chaotic trajectories have been extensively Step 2 : r = 24, a = 0, σ = 10, b = 8/3, Lorenz system,separatrices of S 0 → S ± ; Step 1 : r = 346, a = 0.01, σ = 4, b = 1, Glukhovsky-Dolzhansky system, separatrices of S 0 → S ± ; Step 3 : r = 24, a = 1/r − 0.01, σ = −ar, b = b cr + 0.14, Rabinovich system, separatrices of S 0 → S ± ; performed in the context of the shadowing theory and hyperbolicity breakdown (see, e.g., [29,33,81,84,92]). A computer-generated chaotic trajectory shadows a "true" chaotic trajectory if the former stays uniformly close to the latter and vice versa. Under this assumption, the shadowing of numerical trajectories for a reasonable time span would be a minimum requirement for a meaningful computer simulation of a dynamical system. The obstructions to shadowability vary from mild to severe. For example, hyperbolic chaotic systems have been proved by Anosov to be shadowable for an infinite period of time [3]. However, if nonhyperbolic behavior is caused by manifold tangencies, it has been proved that computer-generated chaotic trajectories shadow "true" trajectories for a finite period of time [91]. Finally, strongly non-hyperbolic systems, like those with unstable dimension variability [38], can be shown to be practically nonshadowable since the shadowing time may be very small. In particular, the existence of fluctuations of the FTLE is a numerical fingerprint of shadowability breakdown via unstable dimension variability (see, e.g., [90,[105][106][107]).
Since the classical Lorenz system (1) has a property of singular hyperbolicity (or quasi-hyperbolicity), only finite shadowability is possible for it (see, e.g., [4,16]). Taking into account these problems, the experiments are carried out on small time intervals and for trajectories on a grid of initial points in the attractor's basin of attraction.  10 6 ].

Fig. 9
The trajectory forms a chaotic set, which looks like an "attractor" (navy blue) and then tends to S + (cyan). (Color figure online) eters (relative tolerance 10 −3 , absolute tolerance 10 −6 ). We numerically approximate the finite-time Lyapunov exponents and finite-time Lyapunov dimension, using the algorithm from [44]. Integration with t > T 1 ≈ 18082 leads to the collapse of the "attractor," i.e., the "attractor" turns out to be a transient chaotic set. However, on the time interval t ∈ [0, T 3 ≈ 6.9 × 10 5 ] we have LE 1 (t, u 0 ) > 0 and, thus, may conclude that the behavior is chaotic, and for the time interval t ∈ [0, T 2 ≈ 3.443 × 10 5 ] we have dim L (t, u 0 ) > 2. Thus, the estimation of duration of the transient behavior by analyzing the sign of the largest Lyapunov exponent could lead to the wrong conclusion. This numerical phenomenon can be explained due to the fact that the finite-time Lyapunov exponents and Lyapunov dimension are averaged during computation over the considered time interval. Since the "lifetime" of a transient chaotic process can be extremely long and in view of the limitations of reliable integration of chaotic ODEs (which we will discuss in the next subsection), even longtime computation of the finitetime Lyapunov exponents and the finite-time Lyapunov dimension does not guarantee a relevant approximation of the Lyapunov exponents and the Lyapunov dimension.
Remark that, if we choose the same initial point u 0 = (20,20,20), as in the previous experiment, but use MATLAB's procedure ode45 with relative tolerance 10 −6 , the chaotic transient process will last over a time interval of more than [0, 10 6 ], and the transient chaotic set will not collapse. If now we choose a new initial point u 0 = (10, 10, 10), it would be possible to observe a chaotic transient process over the time interval [0, T 1 ≈ 6.85 × 10 5 ], while LE 1 (t, u 0 ) > 0 and dim L (t, u 0 ) > 2 would be satisfied on the time intervals t ∈ [0, T 2,3 ], where T 2,3 10 6 .

5.2
The problem of choosing the time interval and initial data: attractor and embedded unstable periodic orbits Consider the challenges of the finite-time Lyapunov dimension computation along the trajectories over large time intervals (see, e.g., [48][49][50][51]), which is connected with the existence of unstable periodic orbits (UPOs) embedded in chaotic attractor. The "skeleton" of a chaotic attractor comprises embedded unstable periodic orbits (UPOs) (see, e.g., [1,6,17]), and one of the effective methods among others for the computation of UPOs is the delayed feedback control (DFC) approach, suggested by K. Pyragas [85] (see also discussions in [15,45,55]). This approach allows Pyragas and his progeny to stabilize and study UPOs in various chaotic dynamical systems. Nevertheless, some general analytical results have been obtained [35], showing that DFC has a certain limitation, called the odd number limitation (ONL), which is connected with an odd number of real Floquet multipliers larger than unity.      In order to overcome ONL, later Pyragas suggested a modification of the classical DFC technique, which was called the unstable delayed feedback control (UDFC) [86].
Rewrite system (1) in a general forṁ Let u upo (t, u ). To compute the UPO and overcome ONL, we add the UDFC in the following form:u (24) where 0 ≤ R < 1 is an extended DFC parameter, N = 1, 2, . . . , ∞ defines the number of previous states involved in delayed feedback function F N (t), λ 0 c > 0, and λ ∞ c < 0 are additional unstable degree of freedom parameters, B, C are vectors and K > 0 is a feedback gain. For initial condition u upo 1 0 and T = τ , we have and, thus, the solution of system (24) coincides with the periodic solution of initial system (23). For the Lorenz system (1) with parameters r = 28, σ = 10, b = 8/3 using (24) with B * = (0, 1, 0), one can stabilize a period-1 UPO u upo 1 (t, u 0 ) with period τ 1 = 1.5586... from the initial point u 0 = (1, 1, 1), w 0 = 0 (Fig. 12). Results of this experiment could be repeated using various other numerical approaches (see, e.g., [13,82,110]) and are in agreement with similar results on the existence of UPOs embedded in the Lorenz attractor [8,25]. However, the Pyragas procedure, in general, is more convenient for UPOs numerical visualization and finds widespread application in a rich variety of chaos control problems (see, e.g., [12,23,94]).
For the initial point u ), we numerically compute the trajectory of system (24) without the stabilization (i.e., with K = 0) on the time interval [0, T = 100] (Fig. 12b). We denote it bỹ u(t, u upo 1 0 ) to distinguish this pseudo-trajectory from the periodic orbit u(t, u ) and visualizes a local chaotic attractor A.
We use the adaptive algorithm for the computation of the finite-time Lyapunov dimension and exponents for trajectories on the local attractor A [44]. In order to distinguish the corresponding values for the stabilized UPO u(t, u ) coincide up to the 4th decimal place inclusive on the interval [0, t 1 m ≈ 7τ 1 ]. After t > t 1 m , the difference in values becomes significant and the corresponding graphics diverge in such a way that the part of the graph corresponding to the unstabilized trajectory is lower than the part of the graph corresponding to the UPO (see Figs. 13b, 14).
The Jacobi matrix at the saddle-foci equilibria S ± has simple eigenvalues, which give the following: dim L S ± = 2.0136. The UPO u upo 1 with period τ 1 = 1.5586 has the following Floquet multipliers: ρ 1 = 4.7127, ρ 2 = 1, ρ 3 = −1.19 × 10 −10 and corresponding Lyapunov exponents:   ). Also remark that system (1) has the analytical solution u(t) = (0, 0, z 0 exp(−bt)) which tends to the equilibrium S 0 = (0, 0, 0) from any initial point (0, 0, z 0 ) ∈ R 3 . In general, the existence of such solutions in the phase space complicates the procedure of visualization of a chaotic attractor (pseudo-attractor) by one pseudo-trajectory with arbitrary initial data computed for a sufficiently large time interval (see, e.g., [79,97]). In particular, the numerical computation of finite-time local Lyapunov exponents along this trajectory during for any time interval does not lead to aver-aging of these values across the attractor, but to tending of these values to the finite-time local Lyapunov exponents of S 0 .

Rigorous analytical computation for the global attractor
Using an effective analytical technique, proposed by Leonov [40,56], it is possible to obtain [65,70] the exact formula of the Lyapunov dimension for the global attractor A glob of the Lorenz system (1): for the case, when conditions (3) of the Theorem 1 are not satisfied. For the Lorenz system (1) with classical values of parameters r = 28, σ = 10, b = 8/3, we have the following relations: Here, since the global Lorenz attractor contains a period-1 UPO: A glob ⊃ u upo 1 , we have the following lower-bound estimate for the Lyapunov dimension: dim L A glob ≥ 2.0678 = dim L u upo 1 . Similar experiment and results for the Rössler system [88] are presented in [48,50].
Concerning the time of integration, remark that while the time series obtained from a physical experiment are assumed to be reliable on the whole considered time interval, the time series produced by the integration of mathematical dynamical model can be reliable on a limited time interval only due to computational errors (caused by finite precision arithmetic and numerical integration of ODE). Thus, in general, the closeness of the real trajectory u(t, u 0 ) and the corresponding pseudo-trajectoryũ(t, u 0 ) calculated numerically can be guaranteed on a limited short time interval only.
In our experiment, if we continue computation over a long time interval [0, 10000] of FTLD along the stabilized UPO and the pseudo-trajectory obtained without Pyragas stabilization, as a result, completely different values will be obtained (Fig. 14). Evolution of dim L (u(t, ·), u upo 1 0 ) along the stabilized UPO will tend to the analytical value dim L u upo 1 = 2.0678, computed via Floquet multipliers, while evolution of dim L (ũ(t, ·), u upo 1 0 ) along the pseudo-trajectory will converge to the value 2.0622 6 . These results are in good agreement with the rigorous analysis of the time interval choices for reliable numerical computation of trajectories for the Lorenz system. The time interval for reliable computation with 16 significant digits and error 10 −4 is estimated as [0, 36], with error 10 −8 is estimated as [0, 26] (see [36,37]), and reliable computation for a longer time interval, e.g., [0, 10000] in [75], is a challenging task that requires a significant increase in the precision of the floating-point representation and the use of supercomputers.

Conclusion
In this work, for the classical Lorenz system, the boundaries of practical global stability are estimated and the difficulties of numerically studying the birth of selfexcited and hidden attractors, caused by the loss of global stability, are discussed. Difficulties in localization and numerical analysis of hidden transient chaotic sets in the Lorenz system are emphasized in the experiments on computation of the largest finite-time Lyapunov exponent and finite-time Lyapunov dimension along trajectories with different initial data and large time intervals. The classical self-excited Lorenz attractor is considered, and the applications of the Pyragas time-delayed feedback control technique and Leonov analytical method are demonstrated for the Lyapunov dimension estimation, as well as for the verification of the famous Eden's conjecture on the maximum value of this characteristic that could be achieved on the attractor. The problem of reliable numerical computation of the finite-time Lyapunov dimension along the trajectories over large time intervals is discussed.
Acknowledgements Open access funding provided by University of Jyväskylä (JYU). The authors declare that they have no conflict of interest. We dedicate this paper to the memory of Gennady A. Leonov (1947Leonov ( -2018, with whom we began this work in 2018. This study was partially funded by the Russian Science Foundation (Project 19-41-02002).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Appendix A: Hidden attractors localization via numerical continuation method
One of the effective methods for numerical localization of hidden attractors in multidimensional dynamical systems is based on the homotopy and numerical continuation method (NCM). The idea is to construct a sequence of similar systems such that for the first (starting) system the initial point for numerical computation of oscillating solution (starting oscillation) can be obtained analytically. Thus, it is often possible to consider the starting system with self-excited starting oscillation; then, the transformation of this starting oscillation in the phase space is tracked numerically while passing from one system to another; the last system corresponds to the system in which a hidden attractor is searched.
For studying the scenario of transition to chaos, we consider system (23) with f (u) = f (u, λ), where λ ∈ ⊂ R d is a vector of parameters, whose variation in the parameter space determines the scenario. Let λ end ∈ define a point corresponding to the system, where a hidden attractor is searched. Choose a point λ begin ∈ such that we can analytically or numerically localize a certain nontrivial (oscillating) attractor A 1 in system (23) with λ = λ begin (e.g., one can consider an initial self-excited attractor defined by a trajectory u 1 (t) numerically integrated on a sufficiently large time interval t ∈ [0, T ] with initial point u 1 (0) in the vicinity of an unstable equilibrium). Consider a path 7 in the parameter space , i.e., a continuous function γ : [0, 1] → , for which γ (0) = λ begin and γ (1) = λ end , and a sequence of points {λ j } k j=1 on the path, where λ 1 = λ begin , λ k = λ end , such that the distance between λ j and λ j+1 is sufficiently small. On each next step of the procedure, the initial point for a trajectory to be integrated is chosen as the last point of the trajectory integrated on the previous step: u j+1 (0) = u j (T ). Following this procedure and sequentially increasing j, two alternatives are possible: the point of A j are in the basin of attraction of attractor A j+1 , or while passing from system (23) with λ = λ j system with respect to a set A is defined as: The Douady-Oesterlé theorem [20] implies that for any fixed t > 0 the finite-time Lyapunov dimension on a compact invariant set A, defined by (28), is an upper estimate of the Hausdorff dimension: dim H A ≤ dim L (t, A). The best estimation is called the Lyapunov dimension [40] dim L A= inf