{'fraction_non_alphanumeric': 0.03489810534540009, 'fraction_numerical': 0.006940556025782928, 'mean_word_length': 4.6671832337096895, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Jacobian deformation ellipsoid and Lyapunov stability analysis revisited

The stability analysis introduced by Lyapunov and extended by Oseledec is an excellent tool to describe the character of nonlinear n-dimensional flows by n global exponents if these flows are stable in time. However, there are two main shortcomings: (a) The local exponents fail to indicate the origin of instability where trajectories start to diverge. Instead, their time evolution contains a much stronger chaos than the trajectories, which is only eliminated by integrating over a long time. Therefore, shorter time intervals cannot be characterized correctly, which would be essential to analyse changes of chaotic character as in transients. (b) Moreover, although Oseledec uses an n dimensional sphere around a point x to be transformed into an n dimensional ellipse in first order, this local ellipse has yet not been evaluated. The aim of this contribution is to eliminate these two shortcomings. Problem (a) disappears if the Oseledec method is replaced by a frame with a 'constraint' as performed by Rateitschak and Klages (RK) [Phys. Rev. E 65 036209 (2002)]. The reasons why this method is better will be illustrated by comparing different systems. In order to analyze shorter time intervals, integrals between consecutive Poincare points will be evaluated. The local problems (b) will be solved analytically by introducing the symmetric 'Jacobian deformation ellipsoid' and its orthogonal submatrix, which enable to search in the full phase space for extreme local separation exponents. These are close to the RK exponents but need no time integration of the RK frame. Finally, four sets of local exponents are compared: Oseledec frame, RK frame, Jacobian deformation ellipsoid and its orthogonal submatrix.


Lyapunov method revisited
Lyapunov exponents [1] are a well-established tool [2][3][4][5][6][7][8][9] to analyse the type of chaos displayed by a trajectory x(t) as a solution of a nonlinear ndimensional equation in the phase space R{x}, . Their evaluation is based on proofs by Oseledec [2] that a vector u to any arbitrarily chosen neighbouring point of x (except the instantaneous direction along dx as discussed later) will rotate to the direction of integrated extreme expansion. After an initial period, this local direction obtained through integrated rotation is an inherent feature of each point on the trajectory. Moreover, this local uniqueness applies to the set of all n directions of a frame orthogonalised subsequently after each time step. The integration for infinite time t → ∞ of the corresponding local exponents results in a set of n Lyapunov exponents assessing the character of chaos in terms of dynamical instability of the given dynamical system.
The fact that an initial frame rotates after a short time to a locally unique orientation has fascinated many authors. This fascination includes also the corresponding local exponents which has, probably, prevented one to ask: Is the local exponent leading to the largest Lyapunov exponent already locally assessing the strongest expansion, i.e. the strongest divergence of neighbouring trajectories, thus identifying the places where trajectories 'break away', the origin of instability and chaos? Even a quick look at this exponent for the Lorenz model reveals that this is not the case resulting in misleading results for shorter time intervals.
The origin of this discrepancy is inherent to the idea of Lyapunov and Oseledec: They considered a sphere of neighbouring points. Each vector u not orthogonal to the direction dx of the trajectory contains a component parallel to dx. This component measures only the expansion along the trajectory, hence the acceleration, which is not connected to instability. Only the component orthogonal to dx analyzes the change of the distance to neighbouring trajectories and can indicate where trajectories start to diverge. This feature is obscured by the admixing of the acceleration. Since in many cases the average acceleration is zero, the integration for long times cancels this contribution, thus only the value for divergence remains. These shortcomings can easily be avoided by using a method introduced by Rateitschak and Klages (RK) [10]: Select as an initial direction u the direction of dx, thus u is parallel to the trajectory. Certainly, this direction will remain parallel to the trajectory for all times. The integral of the corresponding exponent will only describe the average acceleration which for many cases is zero, i.e., if the flow remains bounded and does neither shrink nor blow up. This exponent is not related to the problem of instability. However, all remaining directions of the initial frame, subsequently orthogonalised and their rotation integrated, are describing the divergence of neighbouring trajectories. Their corresponding exponents are the important numbers for analysing the nature of chaos. After an initial period much shorter than for the Oseledec method these directions are unique to the point x(t). The first in the orthogonalizing procedure of the corresponding local exponents is the largest one quantifying where trajectories diverge. Moreover, integrals for shorter time intervals analyse the chaotic character of the trajectories thus discriminating between intervals of strong and weak or no chaos.
This problem of finding coordinate systems in which trivial eigendirections are eliminated has already been addressed by Eckhardt and Wintgen in 1991 [11]. However, they focussed on periodic orbits in conservative two degree of freedom Hamiltonian systems for which they eliminated the two trivial neutral directions along the orbit and perpendicular to it on the energy shell by also requiring certain smoothness properties. Their Hamiltonian method was reviewed and further amended by Gaspard for calculating local Lyapunov exponents and local stretching rates [12]. Other simple numerical methods for computing local Lyapunov exponents and stretching rates have been proposed and tested by Dellago and Hoover [13] and by Rateitschak and Klages (RK) [10].
Both Oseledec and RK methods will be illustrated for the Lorenz [14] and the Rössler [15] model by comparing them with each other and by also applying them to shorter time intervals for transient chaos. It is one of the main points of this contribution to show that the method using a set orthogonal to the flow is more adequate to describe the chaotic behaviour than the Oseledec method using general directions, thus explaining why Rateitschak and Klages [10] found much improved results introducing this method for complex chaotic behaviour.

Jacobian deformation ellipsoid
In order to understand the definition and use of what later on we call 'Jacobian deformation ellipsoid' we first briefly review the origins of Lyapunov instability analysis: (i) The basic idea of Lyapunov exponents is to follow the evolution of points close to the points x(t) on a trajectory. First, these neighbouring points are chosen on an n-dimensional sphere around the starting point x(to). Then this sphere is continuously deformed during the time evolution. After a sufficiently long integration time τ this deformed object is analysed. From the largest expansion direction the largest Lyapunov exponent is evaluated. Starting with the direction of this largest expansion, further orthogonal directions are used to measure their expansions, which complete the set of n Lyapunov exponents.
(ii) In principle, the evolution of the deformation should be numerically computed for an infinite number of points on the initial sphere, a tremendous task even for large computers.
(iii) However, the works of Lyapunov [1] and Oseledec [2] propose a well-established, much less elaborate method: Arbitrarily defined at the start, choose an orthogonal frame of only n directions, follow their evolution -orthogonalised always in the same order by the Gram-Schmidt orthonormalisation method before each integration step -and measure their expansion, thus disregarding their rotation. The Lyapunov exponents are then obtained as the averages of the logarithms of these expansions.
It would be interesting to test to which extent (iii) corresponds to (i) in case (ii) could be treated with less numerical effort. It is the aim of this section to propose a simpler method.
Before using the computer for doing (ii), let's go back to the 19 th century of Jacobi. The main idea was then: Think of a reduction of the problem and retain only the essentials. Jacobi [16] used only the first derivative to define his matrix J at a point x, J(x) = dX / dx (2) . He realised that the sum of the diagonal elements of J measures the local rate of change of the volume of a sphere around x in the limit of infinitesimally small radius. However, nothing has been said then about the deformation of this sphere. Although finally global quantities should be evaluated, we will first consider 'instantaneous local' quantities defined by the deformation between a fixed time t and t+dt, starting as a sphere at t.
Applying the 19 th century method to the deformed sphere, what is the first approximation? This has been shown already by Farmer et al. [4] with a picture of an ellipse in two dimensions. Green and Kim [7] describe a general ellipsoid for n dimensions, which is continuously deformed in time but remains always an ellipsoid. The problem (i) is then solved by using the n principal axes and their corresponding expansions for the final ellipsoid. Obviously, these principal axes are orthogonal in accord with (iii), since the ellipsoid has inversion symmetry.
An n-dimensional ellipsoid is described by a symmetric n x n matrix E, with orthogonal eigenvectors as principal axes and real eigenvalues, and is determined by n(n+1)/2 parameters. At this point it seems worthwhile to note that E yields the radius for the infinite number of possible directions as diagonal elements (EU)ii of EU, after transforming to a new coordinate system with the unitary matrix U and its transpose U T as EU=U T E U. The matrix E has the same size as J. However, J is in principle not symmetric causing complex eigenvalues and non-orthogonal complex eigenvectors. Furthermore, its diagonal elements are connected with the rate of change of the radius of a sphere, not with the radius of the ellipsoid described by E. Therefore, why not try to symmetrise J while keeping the diagonal elements by defining a symmetric matrix S = ½ (J + J T ) with orthogonal eigenvectors and real eigenvalues in order to describe the rates of change of the radii when the sphere is transformed in first order into an ellipsoid? This paper aims to convince the reader that S is exactly describing in first order the 'instantaneous local rates of exponential stretching ratio' -in short 'local exponents' -for all possible directions, although only n 2 numbers are involved in J.
This symmetric result has almost been found by Greene and Kim [7] in their eqs. (26)-(28). They showed that their instantaneous local expansion exponents λU along any orthogonal set of directions U are given by the diagonal elements Kii of K=U T J U (note that in general J and K are not symmetric). Simply reduce this equation for only one direction u to the scalar product λu= (u,Ju). This form can be derived in a direct geometrical way providing probably one of the simplest approaches to define local Lyapunov exponents. The derivation is so short that it will be sketched here in the introduction, as follows: The local exponent λ for a dynamical system given by eq.(1) can be defined as λ= (1/Δt) ln (r), with r = (d+Δd)/d = 1 + Δd/d being the stretching ratio for the length d of a vector u pointing from a point x(t) towards a neighbouring point. Using the Jacobian matrix J, the new vector after time Δt is found in first order to be u(t+Δt) = u(t)+J u(t) Δt. This new vector has rotated and changed its length to d+Δd. It is now important to eliminate the rotation by projecting the change Δu = J u Δt onto u by using the scalar product Δd = (u,Ju) Δt / │u│ 2 . valid for Δt→0. Putting this into r = 1 + Δd/d and expanding ln(1+α) ≈ α for│α│<<1 results in λ= (1/Δt) ln (r) ≈ (1/Δt) (u,Ju) Δt / │u│ 2 . The simple result λ=(u,Ju) is then found if the vector u has been normalised to unity. Note that the `logarithmisation' is no longer visible in the form for λ after the above approximation has been applied. Furthermore, note that this result has the correct dimension of reciprocal time as needed for a Lyapunov exponent defined above as 'the rate of exponential stretching ratio'. Interestingly, the form λu= (u,Ju) shows that u and -u give the same value implying inversion symmetry. In more detail, the elements Jik of the matrix J appear in pairs (Jik + Jki) in the scalar product (u,Ju). Hence J can be replaced by the symmetric matrix S = ½ (J + J T ), with J T denoting the transposed matrix. This gives the same value for λ u= (u,Su) as λ u= (u,Ju)=(u,J T u).
Going now back to the above form K=U T J U of Green and Kim [7] and replacing there J by S, the diagonal elements Tii of the new matrix T=U T S U are equal to Kii. In contrast, this new T is a symmetric matrix corresponding to the symmetrised K with T = ½ (K + K T ). Furthermore, the transformation of J with K=U T J U could be interpreted as transforming J into the new reference system U and the result might be denoted by JU. Similarly, T =U T S U transforms S into SU= U T S U with its diagonal elements (SU)ii as instantaneous local exponents in the U directions. Thus S is the generating form containing all deformations of the ellipsoid and will be called 'Jacobian deformation ellipsoid' (not to be confused with 'deformed ellipsoid E' which describes directly the radii, not the rates of their changes). Note that using the symmetry argument is a consequence of the reduction of the deformation to the lowest order. In reality, the deformations are more complicated than being captured in lowest order. As a by-product, the largest eigenvalue of the main principal axis of S gives the extreme rate of expansion of the sphere, its eigenvectors rarely being parallel to any of the frames in Section 1. 1 The above first step defines the instantaneous local deformation. Before doing the second step, let us recall here three relevant features of the 'Jacobian deformation ellipsoid' S: First, the matrix S = ½ (J + J T ) is obviously symmetric by definition. Second, the `logarithmisation' is not visible explicitly since the approximation ln(1+α) ≈ α for small α has been used; nevertheless, S describes the exponential rates of instantaneous local stretching. Third, S yields this rate for all directions around the point x, although only n 2 numbers of the elements of the local Jacobian matrix J(x) are needed to define the matrix S.
The second step now consists of evaluating global quantities both for schemes (i) and (iii) above in order to test them in the long-time limit. In both cases, global exponents are defined as final deformations after integrating all deformations during a time interval τ=tfinal -tinitial for τ→∞.
Before going into further detail, two important questions have to be answered: Have `deformations of previous deformations' to be evaluated? And have rotations of the ellipsoids to be incorporated?
The answers to both questions are illustrated in figure 1 of Benettin et al. [3]: After each time step the new deformation does not account for the previous deformation, that is, the deformation of the previous deformation is of higher order and can be neglected. Thus instantaneous local quantities can be averaged if integration is replaced by small but finite time steps as in numerical work. As far as rotation is concerned, each new start in Benettin's figure 1 is from a rotated direction. Since S does not contain any information about the rotation, it has to be incorporated separately, thus at each time step the transformed SU(t) = U T (t) S (t) U(t) has to be evaluated. Note that for U(t) any of the two orthonormalised frames of section 1.1 could be used. Furthermore, note that SU(t) has usually nonzero off-diagonal elements implying that the diagonal elements are no eigenvalues. Now, on the basis of the new Jacobian deformation ellipsoid S, differences between (i) and (iii) will be described. The evaluation according to (i) has to be done in two steps: First, average all local matrices SU(t) along a typical trajectory. Obviously, each of the n 2 elements of SU(t) has to be averaged separately. Hence this averaging results in a final symmetric matrix Sfinal. In a second step, only for this final matrix Sfinal the eigenvalues and the eigenvectors have to be evaluated corresponding to global Lyapunov exponents and Lyapunov directions, respectively.
The recipe for performing (iii) is much simpler: The averages of the n diagonal elements (SU)ii(t) of the local matrix SU(t) along a typical trajectory are the global Lyapunov exponents. Again obviously, each of the n diagonal elements of (SU)ii(t) has to be averaged separately. Comparing both methods, the results of (iii) are already incorporated in the final matrix Sfinal as diagonal elements. Hence, to test the equivalence between (i) and (iii) it is sufficient to compare the values of the diagonal elements of Sfinal with its eigenvalues. This test has been performed numerically for the 'stable' Lorenz chaos. There are small but distinct deviations for the free running rotation of Oseledec. The test is successful for the constraint rotation of RK.
So far the literature has mainly focused on local and global Lyapunov exponents. Both quantities have been defined above within a new approach, and tests of this concept will be described later on in this paper. So what else? Previous concepts of Lyapunov instability are only appropriate for 'stable' chaos, where a trajectory does not change its character in time. They are not suitable to analyse transients, crises, or continuous changes of parameters in time in equations of motion. However, an adequate method is easy to find: Instead of only considering 'global quantities' defined for an infinite interval of time τ, try a series of successive finite time intervals τn, each starting at tn and ending at tn+1. It is essential to define the successive times tn such that the resulting data correspond to the character of the trajectory. To give an example, for a 'stable' Lorenz chaos time intervals τn between successive Poincaré points are convenient, which have been chosen such that each interval describes a loop on the left or on the right hand side of the strange attractor. Again both methods (i) and (iii) are compared for different frames of rotation. Only the RK frame produces rather small but distinct deviations of the values for (i) and (iii). These discrepancies are caused by residual acceleration parts parallel to the flow dx. It is easy to eliminate all these directions by forming the n-1 dimensional subspace S┴ (2…n) orthogonal to the flow of each instantaneous local SU(t) with U as the RK frame, a straightforward evaluation. For (i) the eigenvalues of the averaged subspace are then compared with (iii) as the averages of the exponents orthogonal to dx of RK in the same time interval. Indeed, both values are equal within numerical accuracy for all time intervals. The method (iii) based on a series of finite intervals is then successfully applied to a transient Lorenz chaos. It distinguishes well between an initial period of weak chaos, an intermediate period of strong chaos, and a final spiralling onto a stable fixed point. Also here the RK frame is superior to the Oseledec frame.
As a by-product, each instantaneous local subspace S┴ (2…n) -to be calculated at any point x of the phase space R{x} -has n-1 eigenvalues. The largest eigenvalue within this subspace accounts for the extreme rate of divergence between neighbouring trajectories, not obscured by partial acceleration. This yields a novel indicator for extreme local divergence, which can be used to find 'hot spots' of maximum local dynamical instability in the whole phase space.
These applications demonstrate that the new 'Jacobian deformation ellipsoid' is a powerful tool, which furthermore opens a pedestrian approach to defining both local and global Lyapunov exponents.

Comparison of four types of instantaneous local exponents
The novelty of this article is that it compares four types of local directions and exponents with each other. The first two are strictly local to x. They use only the knowledge of the local Jacobian J and the value of dx, there is no need to evaluate any trajectory by integration: 1. The instantaneous local extreme expansion exponent of a local sphere and its corresponding direction are found as the largest eigenvalue and eigenvector of S. Note that this yields the maximal possible value of the exponent for all four different methods, and the direction of this maximal deformation is rarely parallel to the direction obtained by the other three methods.
2. The extreme exponent for the instantaneous local divergence of neighbouring trajectories measured orthogonal to dx and its direction are found as the largest eigenvalue and eigenvector of the subset S┴ .Note that this exponent is the maximal possible value of the largest exponent for divergence in what we will call the W frame of RK as described in method 4 below.
The other two types are evaluated by method (iii). and need integrations both of a trajectory and the directions of the frames. 3. The standard method of Oseledec (O) with an integrated free running frame called V frame produces n instantaneous local exponents and directions. 4. The method of RK using an integrated W frame, where the first direction is always constrained to be parallel to dx thus assessing the instantaneous local acceleration parallel to dx (needs no integration). The remaining instantaneous local exponents are orthogonal to dx describing the divergence of neighbouring trajectories, with their corresponding directions found by integration of the W frame.
1.4 A test of how well the exponents of methods 3. and 4. are correlated with the exponents of methods 1. and 2.
A heuristic test uses Poincaré points and the corresponding local exponents. The distribution of these local exponents is tested as a function of the distance between the corresponding Poincaré points. The Lorenz model reveals that the O method has a much larger chaotic character than the RK method. This confirms again that the RK method is more efficient than the O method. Only the RK method is adequate to evaluate meaningful Lyapunov exponents, which furthermore indicate where the instabilities start locally.

Geometric interpretation of local expansion
Local expansion of a vector u to a neighbouring point of x(t) during a time interval Δt will be described geometrically. In first order, using the Jacobian matrix J the new vector u(t+Δt) will have rotated and changed its length from d to d+Δd. If only the change Δd of the length is considered, the rotation can be eliminated by projecting Δu onto u by using the scalar product Δd ≈ (u , Ju) Δt/ │u│ 2 valid for Δt→0. With the ratio r of stretching r = (d+Δd) / d = 1 + Δd/d the local exponent λu for the direction of u is found according to the form, see Greene and Kim [7], Using ln(1+α) ≈ α for │α│ << 1, the result is the scalar product if the vector u has been normalized to unity. Greene and Kim [7] published in their eqs. (26-28) the same result in form of the diagonal elements Kii for an orthonormalised reference frame V and the Jacobian matrix J. Note that in the scalar product (u,Ju) the values Jik and Jki appear in pair sums (Jik + Jki). This symmetry is enforced by the special structure of the scalar product (u,Ju): on both sides there occurs the same vector u. Therefore, the same result is obtained for a symmetrised Jacobian matrix S constructed by adding the transpose J T .
In more detail, the Jacobi matrix J can be decomposed into J = S + D with the symmetric matrix S defined by eq. (7) and the anti-symmetric matrix D with Dik = -Dki defined by the difference D = ½ (J -J T ). The scalar product of eq. (5) thus results in (u,Ju)=(u,Su) + (u,Du). The second term (u,Du) is zero, because the same vector on both sides of the scalar product enforces terms which add up to zero due to the anti-symmetry of D. Therefore, the reduction of (u,Ju) to (u,Su) as performed in eq. (8) is justified. Furthermore, note that (u,Du)=0 implies that the vector Du is orthogonal to u describing only a rotation of u. At first glance, for capturing the rotation of u in eq. (3) it is thus tempting to reduce the term Ju therein to Du. However, the vector Su can have any direction thus incorporating both elongation and rotation before it is projected onto u by the scalar product (u,Su) to find only elongation. Hence, for describing the full rotation of u it is essential to keep both terms (S+D)u in eq. (3).
The symmetrised Jacobian matrix S will pave the way to identify the Jacobian deformation ellipsoid as demonstrated in the next section.

The Jacobian deformation ellipsoid
A local n dimensional sphere around a point x of eq. (1) will be rotated and deformed after a time interval Δt into a complicated geometrical object. Only in linear approximation its form can be described by a symmetric n dimensional ellipsoid with orthogonal principal axes of the principal deformation exponents. Here, a straightforward simple way to find this ellipsoid will be described: The expansion exponent λu in any direction u is found by the scalar product (u ,Su), see eq. (8).
Defining an orthonormal set U and its transpose U T , the transformed matrix SU= U T SU can be evaluated. It is now important to note that the scalar product of eq. (8) implies that only the diagonal elements (SU)ii are equal to the local expansion exponents λi corresponding to the n orthogonal normalised directions ui in U. Note further that the off-diagonal elements (SU)ik with i≠k can have any non-zero value.
Let us first consider the special case that all offdiagonal elements (SU)ik with i≠k are zero. In this case, the diagonal elements (SU)ii= λi could be eigenvalues of S. Since S is by definition symmetric, its eigenvalues αi are real and the corresponding eigenvectors ai are both real and orthogonal. It is then obvious from eq. (8) that the eigenvalues αi can indeed be local expansion exponents. By expanding the arbitrary vector u in eq. (8) into the basis of eigenvectors of S, in complete analogy to the expectation value problem of a quantum mechanical operator [17], one concludes that the eigenvalues αi are the extreme local expansion exponents, denoted by λi (e) . This implies that the eigenvectors ai correspond to the extreme vectors ui (e) of U (e) , which are orthogonal by definition. They transform the matrix S into its diagonal eigenvalue form Sdiag=(U (e) ) T S U (e) . The vectors ui (e) thus define the principal axes of the ellipsoid into which the local sphere is deformed.
Note that in the general case of arbitrary directions of U the corresponding local expansions exponents λi=(SU)ii are not eigenvalues of S. In summary, according to eq. (8) the matrix S can be considered as the generator of deformations in all possible local directions including the principal axes which yield the extreme values [18]. It is therefore justified to call S the n-dimensional 'Jacobian deformation ellipsoid'. Note that S does not describe the deformed sphere. It accounts for the rate of expansion (positive values) or contraction (negative values) in any direction; hence, this 'deformation ellipsoid' can have values of both signs.
It will be convenient to order the exponents α according to their values αk > αk+1, thus the largest first, and the eigenvectors ak, written as columns in the matrix Ax accordingly, its components expressed in the frame of R{x}.
The Jacobian deformation ellipsoid can be written as a symmetric n dimensional tensor T, its explicit form depending on the frame of reference. The simplest way is in the local frame of principal axes, the matrix Tdiag with only the exponents αk in the diagonal. Its trace as the sum of the exponents is clearly the rate of change of the volume of the sphere around x. The more convenient form Tx would be expressed in the reference frame of R{x}, to be found by back transformation, resulting, obviously, in the form of S expressed usually in the reference frame of R{x}, with the trace unchanged as the rate of change of the volume, An alternative frame of reference will be described in the next section.

The orthogonal Jacobian reduced ellipsoid S┴ (2, …n) for divergence
Although the Jacobian deformation ellipsoid contains all local information about neighbouring points, it is essential to find the principal exponents and their directions orthogonal to the flow. Only these directions indicate true divergence, because they are not disturbed by partial acceleration.
First the reason will be explained why these directions are important. Then an arbitrary local frame of reference is introduced serving to find the local frame of extreme divergence. The ellipsoid is transformed into this frame. After its truncation to the n-1 dimensional subspace orthogonal to the flow dx, the new extreme exponents and principal directions are found by solving the n-1 dimensional eigenvalue problem. Although these procedures are straightforward, the matrix operations will be described in more detail in order to facilitate their implementation in programs for Lyapunov exponents.

The mixing of divergence and acceleration
The Jacobian deformation tensor will be described now by Sx(x) at the point x, the subscript denotes that the matrix is written in the coordinate system R{x}. Sx(x) describes the rate of stretching along all possible directions u to a neighbouring point of any of the points x of the phase space without the need to execute any integration of eq. (1). Whereas most of these points belong to a neighbouring trajectory, there is one point exactly on the trajectory through x. This point is found for Δt→0 along the vector X of eq. (1) which describes the flow through x. Defining a flow unit vector f║ =X /│X│, the corresponding expanding exponent φ║ does not describe any divergence of a neighbouring trajectory. Instead, it is the relative acceleration (dv/dt)/v related only to the change of velocity v along the trajectory through x [7].
Hence, any vector u can be decomposed into a vector sum u=c║f║+c┴g┴, where g┴ is the appropriate unit vector in the space orthogonal to f║. Therefore, u has its expansion exponent λu composed of the exponent φ║ related only to acceleration and λ┴g describing only the divergence of the neighbouring trajectory line in the g┴ direction, λu = c║ 2 (f║,Sf║) + c┴ 2 (g┴,Sg┴) = c║ 2 φ║ + c┴ 2 λ┴g (11)

Constructing the orthogonal Jacobian reduced ellipsoid for divergence
Hence, if the sole interest is to find locally the largest divergence of neighbouring trajectory lines, and not to be disturbed by interference with acceleration, the n-1 dimensional subspace orthogonal to the flow direction f║ has to be constructed. This will be performed using an arbitrary local orthonormal set of reference F with the first column as the unit vector f║. Then, construct the remaining n-1 vectors fk by permutation of the components of f║. Finally, use a Gram-Schmidt procedure to make the fk orthogonal as columns of the matrix Fx(x) expressed in the frame of reference R{x} (denoted by subscript x ) at the local point x.
The Jacobian deformation ellipsoid Sx is transformed into the F frame by The first row and column of SF describe the expansion along the flow. The remaining n-1 dimensional reduced square submatix S┴ (2,…,n) contains all expansions in the orthogonal subspace.

The principal local exponents for divergence
The reduced square submatix S┴ (2,…n) is symmetric and has n-1 eigenvalues β┴k and eigenvectors b┴k as the principal perpendicular local exponents and directions, respectively, again both to be ordered according to their values β┴k>β┴k+1, the largest first, and B is the n-1 dimensional matrix containing the ordered vectors b┴ as columns.
It is worthwhile to construct a local reference frame {h} as the matrix H with the first vector as flow direction f║ and the remaining directions as local perpendicular extreme expansion directions b┴k. In the F frame, the first row and column are zero except HF11=1. The remaining submatrix is filled with the matrix B. This frame HF can be transformed into the reference frame of R{x} by the arbitrary matrix Fx Hx = FxHF (13) Although the exponents β┴k were evaluated already by solving the eigenvalue problem of the reduced submatrix S┴ (2,…n) , a general relation as a numerical control could be written using the frame matrix Hx in a matrix product with the diagonal elements (…)ii giving φ║ as the first, and β┴k=1,…,n-1 as the remaining numbers

Exploring the whole phase space for extreme local divergence
Before a specific trajectory is evaluated, it seems worthwhile to explore the phase space by producing an n dimensional map of the principal exponent of local divergence β┴k=1 in order to find 'hot' regions of large divergence or 'cool' regions of missing divergence. However, applying this procedure to the Lorenz attractor the interpretation is not trivial; there are strong 'hot' and very 'cool' regions well outside the strange attractor. Examples will be given later.
Moreover, changing the parameters in eq. (1) could result in a very different behaviour, more, stronger, less or no 'hot spots' in the whole phase space. To test this would be very elaborate by evaluating various trajectories. The procedure described here is much faster.

Exponents following a trajectory
A specific trajectory x(t) is found after choosing a starting point x0 = x(t=0) by integrating eq. (1) starting at x0. In the spirit of Lyapunov, n local exponents could be evaluated if a specific orthogonal set of directions is defined for each x(t). The average of these local exponents will for t→∞ lead to the Lyapunov exponents Λk characterizing the type of trajectory.
Therefore, the problem of finding these local directions is essential. First, the well established method of Oseledec (O) without a constraint will be shortly described. Then the new method of Rateitschak and Klages (RK) with their constraint will be introduced. A comparison and a test will how later on that only the second method should be used, implying a fundamental change for the description of instability.

The Oseledec (O) method for the local frame V
According to Oseledec [2], any arbitrarily chosen orthogonal frame V0 at the origin x0, rotated according to eq. (1), then orthogonalised and normalised after each time step Δt, will become a unique local frame V[x(t)] for each point x(t) of the trajectory after an initial transient time τtransV. The rotation can be evaluated by applying eq. (1) to neighbouring points or by using eq. (3) as proposed by Greene and Kim [7]. The idea is that the first direction, which is never adjusted by the orthogonalising process, will turn to the direction of strongest divergence from the trajectory, independently from its starting direction in V0. The instantaneous local expansion exponents λVk can be found by eq. (5) or (8), and their averages will be the global Lyapunov exponents ΛVk for t→∞. The first exponent ΛVk=1 will be the largest.
At this point it is interesting to note that and Meier [19] found in their numerical analysis of the angles of the V frame that the direction corresponding to the smallest ('most negative') local exponent is always nearly orthogonal to the flow. Small deviations are probably due to finite step integration and numerical limitations.

The Rateitschak and Klages (RK) constraint frame W
Rateitschak and Klages (RK) [10] introduced a new concept for a local frame W. The rotation can be made by eq. (3), and the first vector wk=1 = f║ is always set parallel to X/│X│. The remaining directions are then orthogonalised always in the same order. Also here, after a transient time τtransW a unique frame W[x(t)] for each point of the trajectory x(t) will be defined. The local expansion exponents λWk can be found by eq. (5) or (8) by using W instead of V, and their time averages will be the global exponents ΛWk for t→∞. Now, the first ΛWk=1 =Λ║ will not be the largest. It has a different function: it measures the mean acceleration of the flow, which is zero for many chaotic models. The remaining ΛWk>1 = Λ┴Wk all describe only divergence of neighbouring trajectory lines, hence the second ΛWk=2 will be the largest.
A numerical test with the Lorenz model [14] confirmed that the third direction of the V frame is not only always nearly orthogonal to the flow [20], but also always nearly parallel to the third direction of the W frame, with deviations of the order of the deviations within the V frame. How fast does an arbitrarily chosen frame approach the locally unique frame directions, described by a transient time τtransient? This question is analysed in the Lorenz system by measuring the recovery time τrecover of a suitably rotated frame. The tests starts with a frame V at time t0. At a time t1 >> τtransient a new frame Vrot(t1) is constructed, which is rotated by 90 degrees with respect to V(t1). Both frames are integrated in time, and the difference of the angles between the maximal expansion direction of the flipped Vrot(t) and the unflipped V(t) are plotted in Fig. 1 as a function of time t. Clearly, the recovery time τrecover is much longer for O (top) than for RK (bottom).

Comparing recovery times of Oseledec with Rateitschak and Klages
For RK, the direction of the flow is obviously not rotated at t1 Therefore, only n-1 directions have to readjust. These recovery times τrecover are an indication for the transient times τtrans after t0.
In addition, it seems worthwhile to note that the exponent of the acceleration of the flow has zero recovery time, since its direction X/│X│ is the local value of eq. (1) for each point x of the phase space.

Oseledec vs. Rateitschak and Klages Lyapunov exponents for t→∞
The corresponding values of the Lyapunov exponents Λ for t→∞ are numerically tested for three and four dimensions using the Lorenz and the Rössler model, respectively. Fig. 2 (Lorenz), and   The first three exponents integrated vs. time following the free V frame according to O. [3]. Centre: the same integrated exponents vs. time following the constraint W frame according to RK [10].
The first two exponents are positive indicating hyper-chaotic behaviour; the third is zero for zero mean acceleration. The final magnitudes are nearly equal and within the fluctuations of both frames, see Fig.4.
The fluctuations are more pronounced for the W frame (below), since its recovery time is much shorter, therefore the strong peaks of x3 (see Fig. 3 bottom) are not so well integrated out as for the V frame (above). Note that the third exponent of the free V frame (above) is slightly below zero although it should be zero according to the theory. These interesting results confirm that the largest Oseledec exponent ΛVk=1 describes the divergence, although its local values are mixed with acceleration. For the Lorenz model the second ΛVk=2 and for the Rössler model the third exponent ΛVk=3 tend to a zero value, thus describing the mean acceleration ΛWk=1 of the RK method by averaging out local acceleration and divergence.
Hence, for infinite time and 'stable' chaos both methods are equivalent since the contribution of acceleration is cancelled by the integration. However, as shown later, for transients and chaotic regimes changing their character only the RK method is adequate.

The four sets of directions and exponents
The following local sets of directions have been studied, each with n corresponding local exponents: (1.1.) (A) The principal axes a in A of the 'full' Jacobian deformation ellipsoid S, with local exponents α.
(1.2.) (B) The 'constrained' flow direction f║ and the extreme divergence directions b┴ of the 'reduced' Jacobian submatrix S┴ (2,…n)   Obviously, the extreme local exponents of 'full' (α) and 'reduced' (β) frames are the limiting values of the corresponding exponents of 'free full' (λV) and 'constrained' (λW) frames, respectively. Deviations between these exponents will be large when the corresponding directions of the different frames are very different. Therefore, angles between the main directions in the different frames for the Lorenz model with {x}={x,y,z} are shown in Fig. 5 as functions of time t (left) and of the variable y (right), with the variable y(t) and z(y) of the trajectory at the bottom. On top, the 'full' case indicates that the angle between the principal axes (a1) of the 'full' frame A and the main axis (v1) of the 'free' frame V is a complicated function of the trajectory.

Comparison of 'full' exponents and directions with the 'constrained' case
At the centre of Fig. 5, however, the 'constrained' case is shown with the angle between the extreme orthogonal divergence direction and the main orthogonal axis of the frame W. This angle has only a small variation (centre left). It is also somewhat closer (centre right) to the shape of the attractor (bottom right).

Eliminating the ambiguity of the starting frame
In both the O and RK method the starting frames V0 and W0 are ambiguous. This ambiguity can be removed if the corresponding local frame defined by the Jacobian deformation ellipsoid is used. For the Oseledec frame the set Ax refers to the principal deformation directions. For RK the frame Hx has as a first vector the flow direction, the remaining ones are the local extreme directions of orthogonal divergence. For RK this setting shortens the transient time until the Wx frame is close to the unique local frame, see Fig. 5 centre, whereas Ax and Vx might be very different, see Fig. 5 top, resulting in a long transient.

Comparing the local exponents as functions of time t
The Lorenz model will be used to study local exponents as functions of time t. Fig. 6 (bottom) shows the x component as a function of time for a short time interval. Above, local exponents are displayed for the same time interval, the left side for the 'free' 'full' case, the right side for the 'constrained' 'reduced' case (see Table 1).
For the 'full' case, Fig. 6 (left, top) plots the first 'free' exponent λV1 together with the main principal exponent α1 (fine line) as the local extreme value. The first exponent λV1 is not at all times the largest and rarely has the extreme possible value of α1. The difference (λV1-α1), (left side, second row) has no relation to the x component (bottom). The third row, left, displays the second 'free' exponent λV2, again with α1 (fine line). Remember that for t→∞ the average of the first exponent is positive and the second average approaches zero. Both local exponents λV1 and λV2 are a mixing of divergence and acceleration.
For the 'constrained' case, Fig. 6 (right), which discriminates between divergence and acceleration, Fig  6 (right, top) shows the main principal exponent orthogonal to the flow β┴1, thus indicating the possible maximum of divergence. The first exponent for divergence λ┴W1 is very close to that maximum β┴k, with the small difference (λ┴W1-β┴) displayed below, related to the small deviation of the relative directions, as shown in Fig. 5 centre.
However, the exponent for the local acceleration λ║W= φ║ of Fig. 6 right, third row, can have values nearly twice as large as the exponent λ┴W1 for divergence, but clearly never exceeds the main principal 'full' exponent α1 (fine line).
At this point it is interesting to ask why the global exponents are equal for both frames as demonstrated in Fig. 2, since their local behaviour is quite different. Assuming validity of the conjecture of Moser and Meier [20] that the direction for the smallest exponent is always orthogonal to dx for the V frame and the numerical result that it is parallel to the corresponding direction of the W frame, both smallest exponents are equal locally and globally. Knowing further that one global exponent must be zero for the Lorenz model and that the sum of all three exponents is equal, also the largest global exponents should be equal for the V and W frames.

Why discussing local exponents as functions of phase space?
The temporal behaviour of the four local exponents has been extensively discussed in the previous sections. For a better understanding it seems worthwhile to explore their behaviour in the phase space R{x}. The directions of the Oseledec V-frame were already shown in 3-D plots and discussed by Wolf et al. [5] and Green and Kim [7] for a short part of the Lorenz trajectory. Here, longer parts of trajectories changing the loop several times will be displayed.. First the magnitude of the exponents will be shown as a function of the x component. Then the trajectories will be shown projected onto the yz plane if the values of the associated different exponents exceed a certain limit c. The reason for the spread of values of λV1 in this diverging range is its mixing of divergence with acceleration, which is strongly negative in this range. The local exponent λ║W= φ║ is shown on the second row, right, together with the second exponent λV2 left, again exhibiting a complex structure with respect to the x component. How is the relation between the first and the second exponent? In order to indicate this relation, the row at the bottom displays λV1 as a function of λV2 (left), and λ┴Wk as a function of λ║W= φ║ (right).

Lorenz local exponents as a function of the x component
In summary, in Fig. 7 all plots of the 'free' V frame (left) show a rather complicated structure compared with the plots of the 'constrained' W frame (right). The reason for this will be illustrated in the next figure. In contrast, the first vector w┴1 of the 'constrained' W frame for divergence is by definition always orthogonal to the flow direction. In Fig 8 (top right) the angle of this direction w┴1 within the orthogonal plane relative to the direction b┴1 of the extreme divergence exponent β┴1 is shown to be always small as a function of the x component with a range between -16 to 8 degrees. Therefore, the values of the divergence exponent λ┴W1 are never far from the value of β┴1. Fig 8 (bottom right) shows the values of λ┴W1 as a function of β┴1 with the straight diagonal indicating the limit λ┴W1 = β┴1. A similar diagonal λ┴V1 = β┴1 is plotted in Fig. 8 (bottom left), showing the first 'free' exponent λV1 also as a function of β┴1. All the excess on the left side as compared to the right side has to be averaged out during integration in time t to finally describe only divergence.

Lorenz local exponents of integrated frames vs. local maximum exponent α1
It seems worthwhile to compare in Fig. 9 the local exponents of the integrated frames V (left) and W (right) as functions of the strictly local maximum exponent α1. Obviously, α1 is the upper limit, which is attained at specific points of the trajectory for the first (top) and the second (bottom) exponents, indicating that the local axes of the frames coincide with the main principal axis of the Jacobi deformation ellipsoid S at these specific points of the trajectory. Again the left side (V) has a more complex structure than the right side (W). An interesting question is: Where in the phase space occur the points where the exponents are close to the limiting value? A yz plot reveals that about the same number of such points is widely distributed for the V frames but is concentrated in a few small regions for the W frame.  Fig. 11 Again the trajectory is plotted at equal time intervals in order to show the change of the velocity. The small circle denotes the maximal value. On both figures, the left side displays the 'free' 'full' case, the right side the 'constrained' 'reduced' case (see table 1). The top rows show the strictly local exponents α1 (general maximum) and β┴1 (maximum of divergence only). The centre display the largest exponents of the integrated frames V and W as λV1 (mixing divergence and acceleration) and λ┴W1 (divergence only), respectively. The bottom rows show the second exponents λV2 (mixing acceleration and divergence) and λ║W= φ║ (acceleration only). Comparing the 'free' exponents of the left side with each other, it is obvious that only α1 has clear cut borders of the limit c, whereas centre and bottom exhibit a peculiar pattern of the limit. Furthermore, centre and bottom have regions where both exponents are above the limit.  Finally, and most importantly, only the 'constrained' case right, top and centre, has large values where the real divergence of the two loops occurs. The discrepancy to the left side is easy to understand: Since acceleration is larger than divergence, and the first exponent λV1 is a strong mixing of divergence and acceleration, the pattern bottom right λ║W= φ║ is easy to see in the pattern centre left λV1, and the maximum (small circle) is nearly at the same place. Hence, only the 'constrained' case seems to analyse directly the local chaotic behaviour. The directions of both frames V and W are dependent on integration and not only on the location x. This brings about a dependence on the former sections of the trajectory. Although the exponents are unique for each location x, they are not simple functions of the phase space, since the former section of each location is different. The integration is a nonlinear procedure and, therefore, might have in itself a chaotic behaviour sensitive to small changes in the previous conditions. This implies a chaotic behaviour in addition to the chaotic behaviour of the analysed chaotic trajectory, resulting in the fact that another trajectory nearby might have a very different local exponent caused by a tiny difference at earlier points of that trajectory. However, this additional complexity is very different for the two frames, dependent on their 'transient times' as a measure to 'forget' earlier sections of the trajectory and on the range of changing angles of the frames. This complex behaviour is also present when angles of selected directions of these frames are plotted, as is performed in the next section. Figure 12 shows angles within certain bounds in the yz plane: left the angle of the first vector of the 'free' V frame relative to the flow direction; right the angle between the first orthogonal vector of the 'constrained' W frame relative to the direction of the extreme divergence direction orthogonal to the flow. The same angles were shown in Fig.8, top, but only as functions of the x component.

Local angles as functions of phase space R{x} in the yz plane
The bounds imposed on the left are smaller than 10, between 60 and 120, and larger than 120 degrees, from top to bottom, respectively.
The bounds on the right are much narrower, since this angle is between -16 and -4, between -2 and 2, and between 4 and 8 degrees, from top to bottom, respectively, with the total range being from -16 to 8 degrees, see Fig. 8, top right. Since this range is small and the angle fluctuates on both sides (defined via the zero degree) of the extreme direction, the resulting exponents fluctuate much less than the exponents of the V frame. Furthermore, although the V frame changes sometimes through 90 degrees, see Fig. 12 left centre and Fig. 8 left top, and thus through the plane orthogonal to the flow, its direction within this plane has not to coincide with the orthogonal direction of the W frame at this point. However, since the third directions of both frames have been found to be nearly parallel, this plane is not far from the second direction of the W frame.

Local divergence in the phase space R{x} outside trajectories
Since the 'constrained' exponent β┴1 of true divergence is only a function of the position x in the phase space and does, therefore, not depend on a trajectory, its value can be evaluated directly in the full phase space R{x}. These values are shown in Fig. 13 as a mesh plot for the xy plane at z values 60 (top left) and 20 (bottom left). The strange Lorenz attractor is about along the diagonal x ≈ y where the values β┴1 are low, but with a higher pass between the sections of the two loops.
A simple local indicator of the curvature of trajectories in the phase space can be found easily: For the evaluation of the local acceleration, the new vector J f║ after the time interval dt is projected by the scalar product ( f║, J f║) onto the flow direction f║. The absolute magnitude of the new vector minus the absolute magnitude of the projection is a measure for the deviation at t+dt from the direction f║ at t of the local trajectory and hence a measure of the curvature of the trajectory, This 'deviation' d is shown in Fig. 13 on the right side. In more detail, the Jacobian matrix J can be written as sum of the symmetric matrix S and the antisymmetric D: J =S+D, see sect. 2.1. The evaluated separate action of these matrices onto d reveals a very different behaviour: The symmetric S causes a strongly varying peculiar pattern in the region of the strange attractor, but its action is nearly negligible outside the attractor. In contrast, the anti-symmetric D creates a rather smooth structure, small at the attractor, but increasingly large with larger distance from the attractor. The combined action of these different sources are visible in Fig. 13 on the right side as varying peculiar pattern along the diagonal (S), and smooth increase outside (D).

The interest of integration during shorter time intervals than infinity
At the beginning of analysing chaos, the main interest was to find general numbers defined in the limit of t→∞, which presupposes a 'stable' chaos. Later on, also chaotic motions were analysed where the strength of the chaotic behaviour was changing. Therefore, shorter time intervals were used. The large variation of the local exponents in time made the resulting numbers strongly dependent on the limits where the time intervals start and end. In order to test different methods, a chaotic behaviour with strongly changing character would be welcome. Such a 'transient' chaos will be described in the next section.

'Stable' and 'transient' Lorenz chaos
Until now, the Lorenz chaos with parameters (σ,ρ,β) = (16,40,4) was 'stable' with repelling unstable fixed points at X=0, where trajectories starting nearby would spiral out, join the wellknown double loop strange attractor and stay there in theory to infinity, in practice until the build-up of computing errors will be too high. The top of Figure 14 displays a section of the 'stable' chaos where the trajectories are plotted on the yz plane at equal time intervals to show the variation of the speed as in earlier plots.
A change of the parameter ρ in the Lorenz equation (dy/dt = -x z + ρ x -y) from 40 to 28.165 replaces the unstable fixed points by attracting stable fixed points. Figure 14 bottom displays a peculiar 'transient' Lorenz trajectory starting with several loops on the right side with increasing radius, followed by a chaotic interval with loops on both sides, and a final decay spiralling on the left side to a stable fixed point.

Time intervals τp between Poincaré points for Lorenz trajectories
For the Lorenz chaos, integrating over one loop would be reasonable. Since there are no closed loops, well chosen Poincaré points will be used to define the start and end of a single loop. Poincaré points will be defined here when a trajectory is crossing a plane at y=cp 'from above', i.e. for decreasing values of y. The time interval τp is then defined between two consecutive Poincaré points. The plane parameter cp is chosen at a level where the trajectories are about normal to the Poincaré plane. Figure 14 top shows this plane (cp=40) as a The integration between Poincaré points is therefore performed over one loop. The results are shown in the next sections for the 'free' and 'constraint' case, first for the 'stable' Lorenz chaos. These results will be compared with the integration of all possible local directions, which are evaluated with the local Jacobian deformation ellipsoid.

Description of the idea of integrating all directions
It seems worthwhile to recall here the basic idea behind the definition of Lyapunov exponents, which consists of integrating along a trajectory the subsequent deformations of the surface of a sphere, surrounding a given initial condition, along all direction in phase space. At the final time, the principal axes and eigenvalues of this deformed ellipsoid are considered to be the relevant characteristics of the problem, thus only n exponents are sufficient to describe the type of dynamical instability. It was proved by Oseledec that only the knowledge of the deformations of a frame of n orthogonal directions is necessary to find these final axes. The evaluation of the deformation with the Jacobian ellipsoid S enables to test this established fact numerically.

The problem of subsequent rotation
It is straightforward to evaluate S[x(t)] along a trajectory x(t). However, for calculating Lyapunov exponents from it the rotations of the directions after each time step dt still have to be eliminated. This important fact is not too obvious, since the expression Ju of eq. (3) implies elongation and rotation. However, in order to find instantaneous local exponents describing only elongation, the rotation must be projected out by using the scalar product of eq. (5). Moreover, because this scalar product has the same vector on both sides, it enforces the symmetrisation of the Jacobian matrix J in form of the symmetric Jacobian deformation ellipsoid S. Therefore, a rotation of S has to be performed after each time step dt along the trajectory. In principle, each direction should be rotated according to eq. (1) separately. However, as a crude oversimplification all orientations might be rotated equally according to the rotation of an n dimensional orthogonal frame evaluated with the equation of motion eq. (1). Along these lines, transformations S into the frames V or W will be used. Again, the average of values computed at discrete time steps will approximate the continuous time integration. Let the interval τ be divided into m parts Δt. With xk = x(t0 + k Δt), and first for the frame V, the resulting arithmetic average S (τ) V can be written as Secondly, in complete analogy the average S (τ) W is defined for the RK frame W with the first vector constrained along the flow direction, Thirdly, another approach consists of considering only all the directions orthogonal to the flow as described by the reduced Jacobian ellipsoid S┴ (2,…n) , obviously to be transformed by a reduced matrix W┴ (2,…n) , resulting in the average S┴ (τ) , At the final time of the integration, for each of the above averages the eigenvalues are evaluated and compared with the integrals of the exponents of the corresponding frames. For convenience, the eigenvalues will be ordered with the largest first for the principal axes of the averaged sums of the ellipsoids.
In contrast, the averaged sums of the local exponents of the frames are ordered according to the order during the orthogonal normalization process, which might result in a reversed order, such that the value for the first exponent is smaller than for the second exponent. For the first case S (τ) V (top), large discrepancies are visible. For the second case S (τ) W (second row), there exist rather small deviations. Only for the third case S┴W (τ) (third row), the eigenvalues of the integral of the reduced deformation ellipsoid give the same values as the integral of the instantaneous local exponents of the W┴ frame. The discrepancies might be caused by the oversimplification of rotating all orientations equally. However, the most important result is that the RK method yields values that are equivalent to the method based on computing eigenvalues by integrating all directions in the n-1 space orthogonal to the flow, even when these orientations are rotated equally but within the orthogonal subspace. The second exponent is by definition equal for RK and S║W (τ) . The small but distinct discrepancies between the full S (τ) W (second row) and S┴W (τ) (third row) are easy to understand: Already at each point x the eigenvectors of the reduced matrix S┴ and the direction of the flow dx only accidentally coincide with the principal axes of the full local S

Integrations of all deformations for long times compared to global exponents
The same comparison is made after a long integration time τ when the sums are close to the asymptotic values. First, using eq. (16), the ellipsoid S (τ) V of the Oseledec type V frame of free rotations has the diagonal elements (1.59, 0.002, -22.59). According to the conjecture for Lyapunov exponents, the global ellipsoid should have principal axes parallel to the directions of the V frame. Equally important, the eigenvalues of S (τ) V should be equal to the global exponents evaluated by the integration of the instantaneous local exponents in the local directions of the frame V, resulting in the values (1.38, 0.002, -22.38). Indeed, there is a small but distinct difference for two values. Moreover, the directions of the principal axes of S (τ) V deviate by about 5 degrees from the directions of the V frame.
Secondly, using eq. (17) the ellipsoid S (τ) W has the eigenvalues (1.37, 0.0001, -22.37), exactly as the corresponding global exponents of the W frame (1.37, 0.0001, -22.37). Moreover, the corresponding directions differ by less than 1 degree. Note that for short intervals as described in the previous section, for S (τ) W the corresponding values are slightly different. Hence, for the global values these small discrepancies are levelled out for long integration time τ.
At this point it seems interesting to discuss the question: Are the global 'Lyapunov eigenvectors' fixed in time? The answer is yes and no, depending on the observer. For a static observer at the origin of the coordinate system of R{x}, the eigenvectors are fixed only if no rotating system U(t) is used. In reality, these eigenvectors are rotating as fast as U(t) itself. Only a rotating observer at the origin of U(t) sees the Lyapunov eigenvectors fixed in time. Moreover, although the final ellipsoid is fixed accordingly for this rotating observer, he would recognise that even after a very large interval τ the instantaneous local ellipsoid is still deforming as rapidly as after a short time.
In conclusion, the approach by Lyapunov as already been described in sect. 1.2 (iii) proposes that the global exponents evaluated only for n local orthogonal directions are the eigenvalues of the global deformation ellipsoid, which contains all directions. This idea has been tested numerically. However, it is only confirmed if the constrained frame W of RK is used. the reaction at the beginning and after the flip. Furthermore, note that during the strong chaotic behaviour the integral of the second exponent is sometimes larger than the first exponent.

'Transient' with a sudden change of the frame: Oseledec V frame
Clearly, the Oseledec V frame has problems with delayed reorientation, which could give wrong exponents for finite time intervals.  Fig. 16. Here the frame is also suddenly rotated by 90 degrees at t=7, as shown in the centre for the angle between the first orthogonal direction of the W frame and the direction of the local orthogonal direction of extreme divergence. This flip changes the exponents only for the integral over one loop, as is shown on top for the second exponent, since the recovery time of the angle is small, and one direction along the flow remains fixed according to the definition of the 'constrained' W frame. This implies that the direction for the first exponent is independent of the orientation of the other directions of the W frame, so the flip at t=7 does not affect the first exponent, see + in Fig. 17 top panel,   Fig. 17. Lorenz 'chaotic transient to a fixed point'; Bottom: x(t). Centre: Angle between the first orthogonal direction of the W frame and the direction of the local orthogonal direction of extreme separation. The W frame is artificially rotated by 90 degrees at t = 7. Top: Poincaré integrals of the first two exponents of the W frame denoted by the symbols o and +, respectively. The dots plotted above the o symbols are the integrals of the local extreme orthogonal exponents, which are independent of orientation and flipping of the W frame.

Integrating the local extreme divergence exponents
It is also interesting to compare the RK Poincaré integrals with the integrals of the local extreme orthogonal exponents evaluated as largest eigenvalues of the local reduced Jacobian ellipsoid, see the dots above the symbols for the second exponent in Fig. 17. Since these integrals are for local exponents, they are independent of the orientation of any frame. Their values are always larger than the RK values, but they provide also an excellent description of the trajectory. They are small during the last part of the spiralling dynamics to a fixed point but still slightly positive there, which indicates chaotic behaviour. However, the negative integral of the exponent along the acceleration simultaneously indicates non-chaotic behaviour.
Hence the combination of both the local extreme exponent with the exponent along the acceleration leads already to a good estimate of the type of regime, chaotic or non-chaotic. Furthermore, both exponents have zero transient time, they are correct already at the staring time t0, Since no frame V or W is needed, the ambiguity of a starting frame, their integration and their transient times are avoided.

Purpose of the method
Qualitatively, the previous figures seem to favour the RK method: The structure of the timedependent dynamics of the exponents is simpler than the one of the established Oseledec method. Hence, it seems worthwhile to quantitatively compare both methods. There are ways to test the complexity of their local exponents, see Chlouverakis et al. [20]. However, the local exponents themselves are merely tools to analyse the chaos of a given nonlinear dynamics. The question is therefore rather: How well are the given trajectories assessed by the local exponents? How well are the exponents correlated to the trajectories? How can the degree of the associated complexity be evaluated? A method similar to an established method to compare the volatility of financial values will be constructed, namely the evaluation of the standard deviation of the distribution of the values.

Description of the 'Poincaré point correlated deviation' (PPCD)
It seems sufficient to analyse only a typical set of points in a reduced n-1 dimensional subspace. For Lorenz chaos, an xy-plane will be defined by z=constant. Pairs of neighbours of the points in this plane are selected and their distance r is evaluated. Then, the absolute difference Δλ of the corresponding local exponents is evaluated as a function of the distance r of these neighbouring points i,k.
Δλik(r) = │ λi -λk │ with r 2 = (xi-xk) 2 + (yi-yk) 2 (19) Finally, plots of the strictly local exponents are compared to the exponents of the V and W frames which incorporate integrals over previous times. This comparison can be performed in a quantitative way and results in a measure of how well the V and W frames correspond to the trajectories to be analysed. The strength of chaotic behaviour is manifest in the 'pseudo'-irregularity of the differences Δλik(r), which is measurable by its standard deviation, after a possible underlying regular behaviour has been eliminated (we call it 'pseudo'irregular, since all data points are regularly determined by eq. (1), including the rotations of the frames).
The following quantitative procedure will be used to evaluate the distribution of the data depending on the function Δλik(r<rm) up to radius rm. Since there might be a systematic bias as a function of r, it seems appropriate to first get rid of this bias. Therefore, the values Δλik(r) are fitted by a linear approximation resulting in aik(r) for each point ik. Then the standard deviation Q of the resulting difference [Δλik(r) -aik(r)] is considered to be an adequate measure for the correlation of the exponent λ with the trajectories, called here 'Poincaré point correlated deviation' (PPCD), The quantity QPPCD(λlocal) of the strictly local exponents serves as a basis to assess the strength of the chaos of the trajectory. Then, QPPCD(λV-frame) and QPPCD(λW-frame) are compared to QPPCD(λlocal). The resulting fractions fV and fV analyse quantitatively the correlation of both frames with the trajectories, If fframe >>1, the respective frame adds chaotic behaviour to the complexity of the trajectories, which has to be eliminated by long-time integration. For fframe≈1, this frame could also be used to analyse short portions of different chaotic transients.

A peculiar structure of the points in the Poincaré plane
For Lorenz chaos, fixing z=40 defines an xy-plane, which serves for finding Poincaré points where trajectories cross the plane. Such Poincaré points p(t) have already been used for the previous 'Poincaré integrals', defined in terms of time limits between consecutive points for decreasing z-values. Here, both points p (+) (t) and p (-) (t) for increasing and decreasing z, respectively, are evaluated. However, in addition to the times ti also the positions xi,yi are considered. A plot of these Poincaré points in the Poincaré xy-plane reveals a particular structure, see Fig. 18 (top): There are two branches of p (-) points around the centre, one with y>x, denoted by c1, and the other with y<x, denoted by c2. The p (+) points have branches on the left (ℓ) torus and on the right (r) torus. Each branch consists of a narrow band of points with only a small scattering around a smooth curve. It is now straightforward to determine for each point p (-) (ti) if the corresponding trajectory comes from ℓ or r and goes to ℓ or r by analysing the time series p (+) (ti-1), p (-) (ti), p (+) (ti+1). The result is striking: the points c1 display only the sequences ℓ→ℓ or ℓ→r plotted in Fig. 18 (bottom) as circles ○ and x signs, respectively. Moreover, these two possibilities are separated in the Poincaré plane. Similarly, the c2 points consist only of an r→r or r→ℓ sequences, plotted as stars * and + signs, respectively. This symmetry allows restricting the analysis to the c1 branch.
At this point an interesting question can be asked: Where on the trajectory could additional noise most easily change the dynamical behaviour? As an example, perturb the trajectory at z=40 where the ℓ→ℓ sequence is close to the ℓ→r sequence, thus the circles ○ and x are very close, such that the trajectory jumps from one to the other sequence. At this position a noise-induced transition as described by Gassmann [21] would be most easily possible To provide an overview, the local exponents of the c1 points are plotted as a function of the x component in Fig. 19. The separation between the ℓ→ℓ and ℓ→r sequences is marked by a vertical line. The strictly local exponents are: extreme expansion α1 (top left), extreme orthogonal divergence β1 (top right). The second exponent λ║W of the W frame of local acceleration is also a strictly local exponent (bottom right). The central row displays the first exponents λV1 and λ┴W, the bottom row the second exponents λV2 and λ║W of the V frame (left) and of the W frame (right), respectively. It is interesting that all three plots on the right side (extreme orthogonal divergence β1 and W frame) share the same feature: The values of the exponents all differ for the two sequences, the barrier is marked by a horizontal line indicating a complete correlation to the behaviour of remaining on the same left loop, or changing the loop, which certainly implies a larger divergence. This feature is missing for both exponents of the V frame (left centre and left bottom) showing a large spreading.  For each point of the c1 branch all local exponents are stored. Then the distance to neighbouring points is evaluated and registered for the radius r<rmax=0.4, together with the absolute difference │Δλik(r)│ of the corresponding local exponents λ. These values are plotted as functions of r in Fig. 20, arranged similarly to Fig. 19. The pairs ik between ℓ→ℓ points are depicted as circles ○, pairs between ℓ→r as x symbol and the rare pairs between ℓ→ℓ and ℓ→r as stars *. The straight lines in Fig. 20 represent linear approximations to the data points. Only the differences of the acceleration │Δλ║W│(bottom right) exhibit a distinct bias as a function of the radius r, due to the strong dependence of these local exponents along the branch c1, as seen for λ║W in Fig. 19 (bottom right). Clearly, all plots of Fig. 20 show a strongly irregular pattern of the data. However, the magnitude of this spread is very different, see the corresponding scales, resulting in the following QPPCD(λ) values: left, from top to bottom: 0.013, 1.8, 1.8, and right, from top to bottom: 0.017, 0.015, 0.022. The average for the three strictly local exponents (0.013, 0.017, 0.022) is 0.017. The average of the two V exponents is 1.8, and the only W exponent │Δλ┴W │subject to integration (centre right) has a value of 0.015. Although the heuristic PPCD analysis provides a strong oversimplification by using only a small selected portion of the data, its outcome is surprisingly clear.

The PPCD analysis of the Lorenz chaos
In order to visualize this final result, Fig. 20 is potted again as Fig. 21, but now with the same scale for all plots. Obviously, only the Oseledec V frame (left centre and bottom) is adding a strong chaotic complexity to the complexity of the trajectories to be analyzed.

Conclusions
This work suggests that the evaluation of Lyapunov exponents should be revisited. The introduction of the symmetric deformation Jacobian ellipsoid and its submatrix orthogonal to the flow allows the direct determination of the principal exponents and of the extreme local exponent for diverging trajectories for every point in the phase space without the need to integrate along a specific trajectory to find the local 'Lyapunov directions' according to the procedure of Oseledec. Moreover, a fundamentally different approach evaluates the divergence of nearby trajectories by avoiding a mixing with local acceleration. This is performed by adding a simple constraint to the Oseledec frame as already done in Ref. [10]: The first vector remains fixed along the flow direction dx, and any further vectors are subsequently orthogonalised after each rotation, which implies that only the divergence and not any partial acceleration are assessed at every point on the trajectory.
This avoids the additional complexity introduced by the Oseledec method, and the local exponents are in accord with the local divergence of the trajectories. The main problem of dynamical instability: Remain neighbouring trajectories in a tube? is thus solved. The largest exponent indicates locally already where trajectories diverge and where the origins of instability occur, Only the evaluation of the Jacobian deformation ellipsoid is treated analytically, the remaining conclusions are based on studying numerical examples. Hence, a formal revision of the original idea of Lyapunov to evaluate the divergence of nearby trajectories instead of considering neighbouring points remains to be carried out.
We finally remark that there might exist interesting crosslinks between our work and the very recent active field of calculating and understanding so-called Lyapunov modes in interacting many-particle systems, see Ref. [22] for a short review and further references therein. Lyapunov modes refer to the eigenmodes associated with the spectrum of Lyapunov exponents which are closest to zero, projected onto the single particles from which they originate. They were found to form interesting spatio-temporal periodic patterns. The methods developed in our paper, based on solving the eigenvalue problem for the local Jacobian deformation ellipsoid in a suitable local coordinate system, could possibly serve for developing alternative techniques of computing such Lyapunov modes. Using our methods, it might also be interesting to check for spatio-temporal structures in the corresponding distribution of local Lyapunov exponents in such systems.