Covariant Lyapunov vectors for rigid disk systems.

We carry out extensive computer simulations to study the Lyapunov instability of a two-dimensional hard-disk system in a rectangular box with periodic boundary conditions. The system is large enough to allow the formation of Lyapunov modes parallel to the x-axis of the box. The Oseledec splitting into covariant subspaces of the tangent space is considered by computing the full set of covariant perturbation vectors co-moving with the flow in tangent space. These vectors are shown to be transversal, but generally not orthogonal to each other. Only the angle between covariant vectors associated with immediate adjacent Lyapunov exponents in the Lyapunov spectrum may become small, but the probability of this angle to vanish approaches zero. The stable and unstable manifolds are transverse to each other and the system is hyperbolic.

Introduction

Lyapunov exponents measure the exponential growth, or decay, of infinitesimal phase-space perturbations of a chaotic dynamical system. For a D-dimensional phase space, there are D exponents, which, if ordered according to size, λ ⩾ λ, are referred to as the Lyapunov spectrum. The classical algorithm for the computation is based on the fact that almost all volume elements of dimension d ⩽ D in tangent space (with the exception of elements of measure zero) asymptotically evolve with an exponential rate which is equal to the sum of the first d Lyapunov exponents. Such a d-dimensional subspace may be spanned by d orthonormal vectors, which may be constructed by the Gram–Schmidt procedure and, therefore, are referred to as Gram–Schmidt (GS) vectors. The GS vectors are not covariant, which means that at any point in phase space they are not mapped by the linearized dynamics into the GS vectors at the forward images of that point [1]. As a consequence, they are not invariant with respect to the time-reversed dynamics. Due to the periodic re-orthonormalization of the GS vectors only the radial dynamics is exploited for the computation of the exponents, whereas the angular information is discarded.

Although the angular dynamics is not a universal property and may depend, for example, on the choice of the coordinate system [2], it would be advantageous for many applications, to span the subspaces mentioned above by covariant vectors and to study also the angular dynamics of and between these vectors. It has the additional advantage to preserve the time-reversal symmetry for these tangent vectors, a property not displayed by the GS vectors. Recently, an efficient numerical procedure was developed by Ginelli et al. [1] for the computation of covariant Lyapunov vectors. Here, we apply their algorithm to a two-dimensional system of rigid disks.

The choice of hard elastic particles is motivated by the fact that their dynamics is comparatively simple, and their ergodic, structural and dynamical properties are well known and are thought to be typical of more realistic physical systems [3]. Secondly, hard-particle systems in two and three dimensions serve as reference systems for the most successful perturbation theories of dense gases and liquids [4,5]. Finally, the combination of a Lyapunov analysis with novel statistical methods for rare events [6] seems particularly promising for the study of such rare transformations in systems, for which hard core interactions are at the root.

The paper is organized as follows. After an introduction of the basic concepts for the dynamics of phase-space perturbations in Section 2, we summarize in Section 3 the features and our numerical implementation of the algorithm of Ginelli et al. [1] for the computation of covariant vectors and covariant subspaces. In Section 4, the Hénon map serves as a simple two-dimensional illustration. The hard-disk model is introduced in Section 5. In this work we restrict ourselves to 198 disks, a number which is dictated by computational economy, but still large enough to allow the study of Lyapunov modes. In Section 5.1 we study the relative orientations of Gram–Schmidt and covariant vectors, which give rise to the same Lyapunov exponents. Next, in Section 5.2, we compare the localization properties in physical space for these two sets of perturbation vectors. The configuration and momentum space projections of the perturbation vectors – Gram–Schmidt or covariant – are the topic of Section 5.3. The central manifold (or null subspace) and its dependence on the intrinsic continuous symmetries – translation invariance with respect to time and space – is discussed in Section 5.4. Although the null subspace is completely orthogonal to the unstable and stable subspaces, it is essential for a proper understanding of the Lyapunov modes [7,8]. Section 5.5 is devoted to a discussion of these modes and how they are represented by the covariant vectors. In Section 5.6 we compute the angles between the covariant modes and test for tangency between covariant Oseledec subspaces. In Section 6 we conclude with a summary.

Phase space and tangent-space dynamics

The dynamics of a system of hard disks is that of free flight, interrupted by elastic binary collisions. If Γ0 denotes the state of the system at time 0, the state at time t is given by Γ = ϕ(Γ0), where ϕ: X  → X defines the flow in the phase space X. Similarly, if δΓ0 is a vector in tangent space TX at Γ0, at time t it becomes , where Dϕ defines the tangent flow. It is represented by a D × D matrix, where D is the dimension of phase space. A subspace ( of the phase space is said to be covariant ifThis definition also applies to covariant vectors, if ( is one-dimensional. Loosely speaking, covariant subspaces (vectors) are co-moving (co-rotating in particular) with the tangent flow. An analogous relation holds for the time-reversed flow.

Next we consider the decomposition of the tangent space into subspaces according to the multiplicative ergodic theorem of Oseledec [9-11]. Here, we closely follow Ref. [7].

The first part of the multiplicative ergodic theorem asserts that the real and symmetric matricesexist for (almost all) phase points Γ. Here, T denotes transposition. The eigenvalues of Λ+ are ordered according to exp(λ(1))>⋯>exp(λ(ℓ)), where the λ( are the Lyapunov exponents, which appear with multiplicity m(. For symplectic systems as in our case, λ( = −λ(ℓ + 1−, which is referred to as conjugate pairing. Similarly, the eigenvalues of Λ− are exp(−λ(ℓ))>⋯>exp(−λ(1)). The eigenspaces of Λ± associated with exp(± λ() are denoted by . They are pairwise orthogonal but not covariant. If the λ( are degenerate with multiplicity , all multiplicities sum to D, the dimension of the phase space. Since the matrices Λ± are symmetrical, each of the two sets of eigenspaces, , completely span the tangent space,The eigenspaces are not covariant, but the subspacesare. They are, respectively, the most stable subspace of dimension of Λ+, and the most-unstable subspace of dimension of Λ− (corresponding to the most stable subspace of that dimension in the past).

The second part of Oseledec’ theorem asserts that for (almost) every phase-space point Γ there exists another decomposition of the tangent space into covariant subspaces ( (Γ) referred to as Oseledec splitting,For δΓ ∈ ((Γ) the respective Lyapunov exponent follows fromThe subspaces ( are covariant (see Eq. (1)) but, in general, not orthogonal. According to Ruelle [10], they are related to the eigenspaces of Λ±:This equation is at the heart of the construction of covariant vectors according to Ginelli et al. as described in the next section. Furthermore, one can show thatare covariant subspaces.

Numerical considerations

Numerical methods probe the tangent space by a set of D tangent vectors, such that the Lyapunov exponents are repeated with multiplicities, λ1 ⩾ ⋯ ⩾ λ. Here, the lower index is referred to as the Lyapunov index. The relation between the λ( and λ is given bywhere f( = m(1) + ⋯ + m( is the sum of all subspace dimensions up to j.

For notational convenience in the following, the vectors , spanning the tangent space at time t, are arranged as column vectors of a D × D matrix . The same convention is used below for other spanning vector sets such as and .

In the classical algorithm of Benettin et al. [12] and Shimada and Nagashima [13] for the computation of Lyapunov exponents, an orthonormal set of tangent vectors at time t is evolved to a time t ≡ t + τ, (τ > 0),where is the Jacobian of the evolution map taking the phase space point Γ at time t to Γ at time t. The column vectors of at time t generally are not orthonormal any more and need to be re-orthonormalized with a Gram–Schmidt procedure. This gives the matrix with column vectors {}, which form the next orthonormal Gram–Schmidt (GS) basis at time t. These vectors are pairwise orthogonal but not covariant. Each GS renormalization step is equivalent to a so-called QR decomposition of the matrix , where the matrix  is upper triangular [14]. The diagonal elements of  are required for the accumulative computation of the Lyapunov exponents. This procedure is iterated until convergence for the Lyapunov exponents is obtained.

For the computation of a covariant set of vectors {}0 spanning the tangent space for the phase point Γ0 ≡ Γ(0) at, say, time t0, Ginelli et al. [1] start with a well-relaxed set of GS vectors at t0 and follow the dynamics forward for a sufficiently long time up to t = t0 + ωτ, storing and (or, equivalently,  ) for t = t0 + nτ, n = 0, …, ω along the way. At t a set of unit tangent vectors {} is constructed according towhich serve as starting vectors for a backward iteration from t to time t0. The vector will stay in at any intermediate time t, because is the most stable subspace of dimension j for the time-reversed iteration. Arranging these vectors again as column vectors of a matrix and expressing them in the GS basis at time t, one has  =  , where the matrix is again upper triangular with elements . If, at any step n, is constructed from according to  = []−1, Ginelli et al. have shown that  = J and, hence, the respective column vectors of this matrix follow the natural tangent-space dynamics without re-orthogonalization. They are covariant but not orthogonal in general. At this stage of the algorithm, renormalization of is still required to escape the exponential divergence of the vector norms without affecting their orientation. After reaching t0 at the end of the iteration, the vectors point into their proper orientations in tangent space such that, according to Eq. (7), is the most-unstable subspace of dimension f( ≡ m(1) + ⋯ + m( of the tangent space at the space point Γ0, going forward in time. If there are degeneracies (as in the presence of Lyapunov modes to be discussed below), the Oseledec subspace ( is spanned according towhere, as in the following, we omit the arguments for the phase-space point. If there are no degeneracies, . Similarly, the Gram–Schmidt vectors may be expressed in terms of the eigenspaces of Λ−,For non-degenerate subspaces one finds [7,15,16].

The drawback of this algorithm for many-particle systems is the large storage requirement for the matrices and (or, equivalently,  ) for the intermediate times t = t0 + nτ, n = 0, …, ω, because τ must not be chosen too large (containing not more than, say, 20 particle collisions). At the expense of computer time, this can be bypassed by storing the matrices only for times separated by, say, 100τ intervals and recomputing the forward dynamics in between when required during the time-reversed iteration. In this case, also the phase-space trajectory needs to be stored.

A simple example: the Hénon map

To illustrate the foregoing algorithm, we apply it to a simple two-dimensional example, the Hénon map [17],with a = 1.4 and b = 0.3. In Fig. 1 the Hénon attractor is shown (black line), which is known to coincide with its unstable manifold. An approximation of the stable manifold is shown by the dotted lines. At the point 0 the initial GS basis is indicated by the two orthogonal vectors in blue, where one, as required, points into the direction of the unstable manifold. If these vectors are evolved forward in time with the GS method for a few hundred steps, the two orthogonal GS vectors at the point ω are obtained. Taking these vectors as the initial vectors and , the consecutive backward iteration yields the covariant vectors at point 0 indicated in red. As expected, one is parallel to the unstable manifold, the other parallel to the stable manifold at that point.

Fig. 1

The Hénon attractor (black line) and a finite-length approximation of its stable manifold (dotted line) are shown. The red vectors are the covariant vectors at the phase point 0 as explained in the main text. The blue vectors are Gram–Schmidt vectors. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Systems of hard disks

Now we turn to the study of a two-dimensional system of hard disks in a box with periodic boundaries, where the particles suffer elastic hard collisions (without roughness), and move along straight lines in between collisions. The case of rough hard disks is the topic of a forthcoming publication [18].

The Lyapunov instability of hard-disk systems has been studied in detail in the past [19-23]. Here, we are mainly concerned with the differences encountered with the GS and covariant vectors, which, as we have seen, give rise to identical Lyapunov spectra. To facilitate comparison with our previous work, we consider reduced units for which the particle diameter σ, the particle mass m and the kinetic energy per particle, K/N, are unity. Here, K is the total energy, which is purely kinetic, and N denotes the number of particles. Lyapunov exponents are given in units of . If not otherwise stated, our standard system consists of N = 198 particles at a density ρ ≡ N/(LL) = 0.7 and a simulation box with an aspect ratio L/L = 2/11, which is periodic in x and y. The choice of such a small aspect ratio facilitates the observation of the Lyapunov modes to be discussed later. As usual, the total momentum is set to zero.

The state of the system is given by the coordinates and momenta of all the particles,Similarly, an arbitrary tangent vector δΓ – either a Gram–Schmidt vector or a covariant vector – consists of the respective coordinate and momentum perturbations,The time evolution of these vectors and the construction of the map from one Gram–Schmidt step to the next has been discussed before [19,24].

Fig. 2 shows the Lyapunov spectrum for this system computed both in forward direction with the GS vectors (blue line) and in backward direction with the covariant vectors (red line). The typical time of the simulation in the forward direction is for t = 2.5 × 105τ, where τ = 0.6 is the largest interval between two successive Gram–Schmidt re-orthonormalizations, which does not affect the spectrum. The backward simulation is for a time t  − t0 = 2.5 × 104τ. The time t0 (at least of the order of 1 × 104τ) is required for the preparation of the relaxed initial state at t0. It can be observed in the figure that the unstable directions in the future correspond well to the stable directions in the past and vice versa. Of course, if the sequence of covariant vectors is followed in the forward direction of time, the spectrum is identical to the classical GS results (blue line in Fig. 2).

Fig. 2

Lyapunov spectrum for the 198 disk system described in the main text. The spectrum calculated in the forward direction with the GS method is shown by the blue line, the one calculated in the backward direction with the covariant vectors by the red line. Reduced indices i/4N are used on the abscissa. Although the spectrum is defined only for integer i, solid lines are drawn for clarity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Covariant versus Gram–Schmidt vectors

Whereas the time evolution of the GS vectors is determined by the exponential growth of infinitesimal volume elements belonging to subspaces 1 ⊕ ⋯ ⊕  for i ∈ {1, …, D} according to , the growth of an infinitesimal perturbation representing a covariant vector is directly proportional to exp(tλ), for all i. Thus, it is interesting to compare the relative orientations of respective vectors giving rise to the same exponent. In the left panel of Fig. 3 the difference in orientation of the two types of vectors is demonstrated by a plot of ∣·∣ as a function of i. The black line is an average over 100 frames separated by time intervals of 250τ. Since for tangent vectors only their direction and not the sense of direction is important, an absolute value is taken (here and for analogous cases below), otherwise the scalar product might average to zero over long times, with equal numbers of vectors pointing into opposite directions. For the unstable directions in the left half of the left panel, one observes a rapid decrease of the scalar product with i and, hence a rapid increase of the angle between respective covariant and GS vectors. This decrease is repeated for the stable directions in the right half of the figure. These two parts are separated by the mode region, an enlargement of which is shown in the right panel of Fig. 3 and which will be dealt with in more detail below.

Fig. 3

Plot of the scalar product norm, 〈∣·∣〉, for GS and covariant vectors giving rise to the same Lyapunov exponent λ, as a function of i. The line is a time average as described in the main text. Left panel: full range of Lyapunov exponents; right panel: enlargement of the central part.

In Fig. 4 we show similar projections (time averages of absolute values of scalar products as before) for selected covariant vectors with the whole Gram–Schmidt vector set. One observes that the covariant vectors belong to the GS subspace 1 ⊕ ⋯ ⊕ , for all j ∈ {1, … , 4N} and, thus, give rise to the upper-triangular property of the matrix R in the QR-decomposition mentioned above. The curves in the figure strongly depend on the choice of j:

Fig. 4

Time averaged absolute value of the scalar product of selected covariant vectors (as indicated by the labels) with the whole set of Gram–Schmidt vectors as a function of i.

If it belongs to the unstable subspace but does not represent a Lyapunov mode (top-left panel for j = 370), there is no obvious orientational correlation with any of the GS vectors with index i < j. For i = j = 1 corresponding to the maximum exponent, the covariant and GS vectors are identical. If, however, the covariant vector represents a Lyapunov mode as in the bottom-left panel for j = 393, then its angle with the respective GS vector may become smaller, giving rise to a scalar product closer to unity.

If the covariant vector belongs to the stable subspace but does not represent a mode as for j = 700 in the top-right panel of Fig. 4, it has non-vanishing components in the GS basis for all i ⩽ j with the exception of the zero subspace 2N − 2 ⩽ i ⩽ 2N + 3, which is strictly orthogonal. With the exception of the step at the conjugate index i = 4N + 1 − j = 93, the origin of which is not fully understood, there is no indication of orientational correlations between the covariant vector with any of the GS vectors for i ⩽ j. If, however, the covariant vector represents a mode as for j = 400 in the lower-right panel of the figure, there is strong orientational correlation not only with the respective GS vector with i = 400, but also with its conjugate pair at 4N + 1 − i (=393 in our example).

It is interesting to note that the leading GS and covariant vectors in the null subspace are always identical (up to an irrelevant sign): 2 = 2.

Localization

The maximum (minimum) Lyapunov exponent is the rate constant for the fastest growth (decay) of a phase-space perturbation and is dominated by the fastest dynamical events, a locally-enhanced collision frequency. It is not too surprising that the associated tangent-vector components are significantly different from zero for only a few strongly-interacting particles at any instant of time. Thus, the respective perturbations are strongly localized in physical space. This property persists in the thermodynamic limit such that the fraction of tangent-vector components contributing to the generation of λ1 follows a power law ∝ N−, η > 0, and converges to zero for N → ∞ [21,25-27]. The localization becomes gradually worse for larger indices i > 1, until it ceases to exist and (almost) all particles collectively contribute to the coherent Lyapunov modes to be discussed below. Similar observations for spatially extended systems have been made by various authors [22,23,28-30], which were consequently explained in terms of simple models [31,32]. We also mention Ref. [33], where the tangent-space dynamics of the first Lyapunov vector 1 for various one-dimensional Hamiltonian lattices is compared to that for the Kardar–Parisi–Zhang model of spatio-temporal chaos. The unexpected differences found for the scaling properties are traced back to the existence of long-range correlations, both in space and time, in the Hamiltonian chains, the origin of which, however, could not be fully disclosed. The same correlations are conjectured to be responsible for a slow convergence of λ1 towards its thermodynamic limit [33], which is also observed for hard-disk systems [19].

Up to now, all considerations concerning localization were based on the Gram–Schmidt vectors. Here, we demonstrate the same property for the covariant vectors. According to Eq. (11) we define the contribution of an individual disk n to a particular perturbation vector as the square of the projection of δΓ onto the subspace pertaining to this disk,Since δΓ is either a GS vector or a covariant vector both of which are normalized, one has , and μ may be interpreted as a kind of action probability of particle n contributing to the perturbation in question. It should be noted that for the definition of μ the Euclidean norm is used and that all localization measures depend on this choice. Qualitatively, this is still sufficient to demonstrate localization. From all the localization measures introduced [25,30], the most common is due to Taniguchi and Morriss [22,23],Here, S is the Shannon entropy for the “probability” distribution μ, and 〈⋯〉 denotes a time average. W is bounded according to 1/N ⩽ W ⩽ 1, where the lower and upper bounds apply to complete localization and delocalization, respectively. In Fig. 5, we compare W obtained for the full set of Gram–Schmidt vectors (blue curve) to that of all the covariant vectors (red curve). The spectra are obtained by identifying δΓ with all vectors of the respective sets, i = 1,…,4N. Not too surprisingly, the localization is stronger for the covariant vectors, whose direction in tangent space is solely determined by the tangent flow and is not affected by renormalization constraints. Another interesting feature is the symmetry W = W4, which is a direct consequence of the symplectic nature of the flow [34].

Fig. 5

Localization spectra W for the complete set of Gram–Schmidt vectors (blue) and covariant vectors (red). The details of the hard-disk system are given in Section 5. Reduced indices i/4N are used on the abscissa. In the inset a magnification of the central mode region is shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Tangent space projections

It is interesting to see how much the coordinate and momentum subspaces contribute to a particular tangent vector δΓ (see Eq. (11)), which may be a Gram–Schmidt vector or a covariant vector , both associated with the same Lyapunov exponent λ. The time-averaged squared projections of δΓ onto the coordinate and momentum subspaces and , respectively, are given byand are plotted in Fig. 6 for the whole set of Gram–Schmidt vectors, and in Fig. 7 for the whole set of covariant vectors, i = 1, …, 4N. One notes that for the Gram–Schmidt case the contributions of η and η to a vector and its conjugate 4 are interchanged, whereas for the covariant vectors and 4 they are the same. This is particularly noticeable for the expanded central regions in the respective right panels of Figs. 6 and 7.

Fig. 6

Mean squared projections for the full Gram–Schmidt vector set, i = 1, … , 4N, onto the coordinate subspace , (green line), and the momentum subspace , (black line), for the 198-particle system defined above. Left panel: full spectrum; right panel: magnification of the central mode-carrying region. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7

Mean squared projections for the full covariant vector set, i = 1, … , 4N, onto the coordinate subspace , (green line), and the momentum subspace , (black line), for the 198-disk system defined above. Left panel: full spectrum; right panel: magnification of the central mode-carrying region. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Central manifold and vanishing exponents

The dynamics of a closed particle system such as ours is strongly affected by the inherent continuous symmetries, which leave the Lagrangian and, hence, the equations of motion invariant. The symmetries relevant for our two-dimensional system with periodic boundaries are the homogeneity of time (or invariance with respect to time translation), and the homogeneity of space (or invariance with respect to space translations in two independent directions). Each of these symmetries is associated with two vector fields with sub-exponential growth (or decay) and, therefore, gives rise to two vanishing Lyapunov exponents [35]. At any phase-space point Γ, the six vectors span a six-dimensional subspace of the tangent space TX(Γ), which is referred to as null space or central manifold. This subspace is covariant. If the 4N components of the state vector are arranged asthe six orthogonal spanning vectors, which are the generators of the elementary symmetry transformations, are given by [7,21]1 corresponds to a change of the time origin, 4 to a change of energy, 2 and 3 to an (infinitesimal) uniform translation of the origin in the x and y-directions, respectively, and 5 and 6 to a perturbation of the total momentum in the x and y-directions, respectively. The six vanishing Lyapunov exponents are located in the center of the Lyapunov spectrum with indices 2N − 2 ⩽ i ⩽ 2N + 3. The first three of these vectors have non-vanishing components only for the position perturbations in the 2N-dimensional configuration subspace , the remaining only for the momentum perturbations in the 2N-dimensional momentum subspace . They are related by  = J for k ∈ {1, 2, 3}, where J is the symplectic (skew-symmetric) matrix.

Let us consider the projection matrices α and β of the GS and covariant vectors, respectively, onto the natural basis,For i ∉ {2N − 2, …, 2N + 3} these components vanish. Without loss of generality, we consider in the following example a system with only N = 4 particles in a periodic box, which is relaxed for t0 = 106 time units, followed by a forward and backward iteration lasting for t − t0 = 105 time units. Very special initial conditions for the backward iteration, for i = 1, …, 4N(=16), are used. The projections at the time t0 are given in Table 1 for the GS vectors, in Table 2 for the covariant vectors.

Table 1

Instantaneous projection matrix α of Gram–Schmidt vectors (for i ∈ {2N − 2,… ,2N + 3}) onto the natural basis {, 1 ⩽ k ⩽ 6} of the central manifold. The system contains N = 4 disks. The powers of 10 are given in square brackets.

i	αi,1	αi,2	αi,3	αi,4	αi,5	αi,6
2N − 2	−0.766	0.582	0.273	−0.766 [−6]	0.582 [−6]	0.273 [−6]
2N − 1	0.256	−0.114	0.960	0.256 [−6]	−0.114 [−6]	0.960 [−6]
2N	0.590	0.805	−0.062	0.590 [−6]	0.805 [−6]	−0.062 [−6]
2N + 1	−0.611 [−6]	0.782 [−6]	−0.121 [−6]	0.611	−0.782	0.121
2N + 2	0.575 [−6]	0.544 [−6]	0.611 [−6]	−0.575	−0.544	−0.611
2N + 3	−0.543 [−6]	−0.304 [−6]	0.783 [−6]	0.543	0.304	−0.783

Table 2

Instantaneous projection matrix β for the six central covariant vectors onto the natural basis {, 1 ⩽ k ⩽ 6} of the central manifold. The system contains N = 4 particles. The powers of 10 are given in square brackets.

i	βi,1	βi,2	βi,3	βi,4	βi,5	βi,6
2N − 2	−0.766	0.582	0.273	−0.766 [−6]	0.582 [−6]	0.273 [−6]
2N − 1	0.256	−0.114	0.960	0.256 [−6]	−0.114 [−6]	0.960 [−6]
2N	0.590	0.805	−0.062	0.590 [−6]	0.805 [−6]	−0.062 [−6]
2N + 1	−0.611	0.782	−0.121	0.611 [−5]	−0.782 [−5]	0.121 [−5]
2N + 2	0.575	0.544	0.611	−0.575 [−5]	−0.544[−5]	−0.611 [−5]
2N + 3	−0.543	−0.304	0.783	0.543 [−5]	0.304 [−5]	−0.783 [−5]

A comparison of the two tables reveals the following:

The six orthogonal GS vectors ; i = 2N − 2, … , 2N + 3, completely span the null subspace (the squared elements for each rows add up to unity in Table 1). The same is true for the six non-orthogonal covariant vectors ; i = 2N − 2, … , 2N + 3, in Table 2.

The first three covariant and Gram–Schmidt vectors completely agree. This is a consequence of the special initial conditions for the former at the time t as mentioned above. During the backward iteration the three covariant vectors stay in their respective subspaces and remain parallel to the GS vectors (which were stored during the forward phase of the algorithm). At t0 they are still identical to their GS counterparts. The first vectors always agree, , even if less special initial conditions conforming to Eq. (9) are used.

Equivalent components have the same mantissa but may differ by a factors of 105 or 106, which are related to the duration of the relaxation phase t0 and of the forward–backward iteration time t − t0.

The explanation for this behavior [34] is obtained by a repeated explicit application of the linearized maps for the free streaming and consecutive collision of particles [19,24] to the six basis vectors . One finds thatfor k ∈ {1,2,3}. Eq. (21) implies that any perturbation vector with non-vanishing components parallel to 4, 5, or 6 will rotate towards 1, 2, and 3, respectively. It follows (i) that the null subspace is covariant; (ii) that the subspaces ,  and are separately covariant (from Eq. (20)); that, as was already noted in Ref. [7], can be further decomposed into the three two-dimensional covariant subspaces , and .

Lyapunov modes

We have seen in Section 5.2 that the perturbation vectors are less and less localized, the smaller the Lyapunov exponents become, until they are coherently spread out over the physical space and form periodic spatial patterns with a well-defined wave vector . This collective patterns are referred to as Lyapunov modes. The modes were observed for hard-particle systems in one, two and three dimensions [7,20,22,23,36], for hard planar dumbbells [25,37,38] and for one and two-dimensional soft particles [27,39,40]. A formal classification of the modes is given in Ref. [7]. Physically, they are interpreted as periodic modulations with wave number (k ≠ 0) of the null modes associated with the elementary continuous symmetries and conservation laws. Since this modulation involves the breaking of such symmetries, the modes have been interpreted as Goldstone modes [8]. Theoretical approaches are based on random matrix theory [41,42], periodic orbit expansion [43], and kinetic theory [8,44,45].

So far the numerical work on Lyapunov modes has been exclusively concerned with the orthonormal Gram–Schmidt vectors {}, i = 1, … , 4N. The purpose of this section is to point out some differences one encounters, if the modes for the Gram–Schmidt and covariant vectors are compared.

Fig. 8 shows an enlargement of the mode-carrying region for the Lyapunov spectrum of Fig. 2. In order to emphasize the conjugate pairing symmetry λ = −λ4 for symplectic systems, conjugate exponent pairs are plotted with the same index i on the abscissa, where now i ∈ {1, … , 2N}. The open circles are computed from the Gram–Schmidt vectors in the forward direction of time, the dots from the covariant vectors during the time-reversed iteration. Considering the size of the system (N = 198), the agreement is excellent.

Fig. 8

Enlargement of the mode regime for the Lyapunov spectrum depicted in Fig. 2. The open symbols indicate exponents computed from the Gram–Schmidt vectors, the full dots are for exponents obtained from the covariant vectors.

The steps in the spectrum due to degenerate exponents is a clear indication for the presence of Lyapunov modes. According to the classification in our previous work [7], the steps with a twofold degeneracy are transverse (T) modes – T(1, 0), T(2, 0) and T(3, 0) from right to left in Fig. 2. Similarly, the steps with a fourfold degeneracy of the exponents are longitudinal-momentum (LP) modes – LP(1, 0), LP(2, 0) and LP(3, 0) again from the right. The arguments (n, n) account for the number of periods of the sinusoidal perturbations in the x and y-directions. Since our simulation cell is rather narrow, only wave vectors parallel to the x-axis of the (periodic) cell appear, leaving 0 for the second argument [7]. As usual, “transverse” and “longitudinal” refer to the spatial polarization with respect to of the wave-like pattern.

One of our early observations, which greatly facilitates the classification of the modes for the Gram–Schmid vectors [7], is that in the limit N → ∞ the cosine of the angle Θ between the 2N-dimensional vectors of the position perturbations and momentum perturbations converges to +1 for the smallest positive, and to −1 for the smallest negative exponents. See the blue line in Fig. 9. Furthermore, the relationholds with known constants C±. This means that these vectors are nearly parallel or anti-parallel for large N and that the mode classification may be based solely on δ. Somewhat surprisingly, this property does not strictly hold anymore for the covariant vectors. This is shown by the red line in Fig. 9, where cos(Θ) is seen to differ significantly from ±1 for all i outside of the null subspace (for which 394 ⩽ i ⩽ 399). Unfortunately, this has dire consequences for the representation of the covariant vector modes, since they cannot be purely understood as a vector field of the position perturbations only as in the GS case. For the purpose of this paper, however, we restrict to the GS-based classification of Ref. [7].

Fig. 9

Time averaged value of cos (Θ) =  (δ·δ)/(∣δ∥δ∣) as a function of the Lyapunov index i for a system with N = 198 hard disks. Here, δ ∈  and δ ∈  are the 2N-dimensional vectors of all position perturbations respective all momentum perturbation for the Gram–Schmidt vectors (blue line) and the covariant vectors (red line). The insets are magnifications of the mode-carrying region. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Transverse modes are two-dimensional subspaces (for the periodic boundary conditions and a rectangular box), for which two orthogonal basis vectors are given in Table 3. As an example, we show in the panels on the left-hand side of Fig. 10 snapshots of the mode T(1,0) for the index i = 393, namely plots of δq as a function of q (top-left), and of δp as a function of q (bottom-left). The respective plots for the x components fluctuate around zero, as expected, and are not shown. Analogous plots for the mode T(2,0) with i = 387 are shown in the panels on the right-hand side. The blue points are for GS vectors, the red squares for the respective covariant vectors. It is interesting to note that the scatter of the points for the position perturbations is smaller for the covariant modes (red squares) than for the GS modes (blue dots). A fit shows that the residuals for the covariant modes are smaller by about a factor of two in comparison to Gram–Schmidt. Quite the opposite is true for the momentum perturbations in the bottom row of panels. Although the proportionality of Eq. (22) still holds, the scatter of the red squares for the covariant vectors is larger than that of the blue dots for the GS vectors. Such a behavior is always observed and is not simple numerical noise. The reason for this behavior is related to the previous discussion in connection with Fig. 9 but still needs further clarification.

Table 3

Basis vectors of (n,0) modes for a hard-disk system in a rectangular box with periodic boundaries. We use the notation c = cos(kx), and s = sin (kx), where the wave vector is given by  = (k,k) = (2πn/L,0). Here n ∈ {1, 2, 3}.

n	Basis of T(n)	Basis of L(n)	Basis of P(n)
n_x0	0c_x,0s_x	c_x0,s_x0	p_xp_ys_x,p_xp_yc_x

Fig. 10

Instantaneous transverse Lyapunov modes T(1, 0) for index i = 393 (left panels), and T(2, 0) for index i = 387 (right panels) for the 198-disk system. In the panels at the top (bottom) the y-coordinate perturbations δq (y-momentum perturbations δp) of all particles are plotted as a function of their x coordinate, q, in the simulation cell. The wave vector is parallel to the x-axis. The blue dots are for Gram–Schmidt vectors, the red squares for covariant vectors. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Longitudinal (L) and associated momentum (P) modes share the same degenerate Lyapunov exponent λ(, and generally appear superimposed in experimental vectors. With a rectangular box and periodic boundaries, they form four-dimensional LP perturbations. The superposition varies periodically with time. This “dynamics” has been identified as a rotation of the pure L and P vectors in the standard frame. For details we refer to previous work in Ref. [7]. The patterns for the pure L mode are easily recognizable as sine and cosine functions, but those for the P modes are not. As is evident from the spanning vectors for L(1,0) and P(1,0) also listed in Table 3, the P modes are proportional to the instantaneous velocities of all particles, which does not at all constitute a smooth vector field. For a pattern to be recognizable, these velocities need to be “divided out”. A full mode reconstruction is required as is described for the case of Gram–Schmidt vectors in Ref [7]. Here we carry out an analogous reconstruction in terms of the covariant vectors and compare them to the GS modes. In Fig. 11 two of the reconstructed patterns for L and P modes belonging to the four-dimensional LP(1,0) subspace with indices i ∈ {388, 389, 390, 391} are shown. The blue dots are for GS modes, the red squares for covariant modes. To judge the quality of the reconstruction, we have included in the top-right panel also the δq versus q curve, which vanishes nicely as required.

Fig. 11

Reconstructed position perturbations of the pure L(1,0) mode (top panels) and P(1,0) mode (bottom panels). Only the patterns proportional to sin(2πq/L) are shown. Blue dots: Gram–Schmidt vectors; Red squares: covariant vectors. For details we refer to the main text. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

For comparison, Fig. 12 gives results for a completely analogous reconstruction, where instead of the position perturbations as in Fig. 11, the corresponding momentum perturbations are used. For this example cosine patterns were selected, whereas in Fig. 11 sine patterns were used. As before, blue dots refer to GS vectors, red squares to covariant vectors.

Fig. 12

Reconstructed momentum perturbations for the pure L(1,0) mode (top panels) and P(1,0) mode (bottom panels) depicted already in Fig. 11. Only the patterns proportional to cos(2πq/ L) are shown. Blue dots: Gram–Schmidt vectors; red squares: covariant vectors. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Transversality

From the Lyapunov spectrum of Fig. 2 and the magnification of its central part in Fig. 8, the following inequalities are read off,where the equal sign applies for the degenerate exponents belonging to modes. [λ0] = 0 is sixfold degenerate in our case. Conjugate pairing assures that λ = −λ4. The Oseledec splitting provides us with the following structure of the tangent space,where  = 1 ⊕ ⋯ ⊕ 2 and  = 2 ⊕ ⋯ ⊕ 4 are the covariant unstable and stable subspaces, respectively, and is the null subspace or central manifold. The question arises whether the system is hyperbolic, which implies that the angles between the stable manifold and the unstable manifold are bounded away from zero for all phase points (due to the existence of a central manifold this is referred to as partial hyperbolicity in the mathematical literature [46]). Even more, we may ask whether the angles between all Oseledec subspaces and, hence between all covariant vectors, are bounded away from zero for all phase-space points. To find an answer to that question, we compute in the following the scalar products for all covariant vector pairs and present representative results. (This procedure reminds us of the so-called coherence angles introduced by d’Alessandro and Tenenbaum [47,48], measuring the angular distance between a physically interesting direction and the direction of maximum perturbation expansion).

The lines in Fig. 13 depict the product norms 〈∣·∣〉 for selected covariant vectors with all other covariant vectors ,i ≠ j. As before, a time average is performed. The panels on the left-hand side provide three examples for from the unstable manifold outside of the mode regime (j = 1, 200, and 370 from top-left to bottom-left, respectively), and similarly on the right-hand side from the stable manifold outside of the mode regime (j = 420, 600, and 792 from bottom-right to top-right, respectively).

Fig. 13

The lines are time averaged (100 frames separated by 150 time units) absolute values of the scalar product of a specified covariant vector with all the remaining covariant vectors as a function of i. Left panels from top to bottom: j = 1, 200, and 370 belonging to the unstable manifold; right panels from bottom to top: j = 420, 600, and 792 from the stable manifold.

One immediately observes that the stable and unstable subspaces are not orthogonal. As has been mentioned in Section 5.4 and is also convincingly demonstrated in the following Fig. 14, the null subspace is orthogonal to both and . For two covariant vectors from the same subspace, or , however, the scalar product does not vanish indicating considerable non-orthogonality. But at the same time it is also well bounded away from unity, which means that the two vectors do not become parallel either. However, one possible exception may be the covariant vector pairs (, ) for adjacent Lyapunov exponents in the spectrum. In these cases, the scalar product reaches a pronounced maximum in all of the six panels of Fig. 13 which may still allow these vectors to become parallel occasionally. This will be discussed further below.

Fig. 14

The lines are time averaged (100 events separated by 150 time units) absolute values for the scalar product of covariant vectors (as specified by the label j) with all the remaining covariant vectors as a function of i. The abscissa is restricted to the mode regime. j is for modes from the unstable subspace with positive exponents only. For the respective conjugate modes from the stable subspace, the curves in all panels are just the mirror images around the center as in Fig. 13.

So far we have considered only vectors outside of the mode regime. The case of covariant vectors representing modes is treated separately in Fig. 14, where, as before, time-averaged scalar product norms 〈∣·∣〉 for i ≠ j are plotted as a function of i. The panels on the left-hand side are for j belonging to unstable transversal modes, the panels on the right-hand side for j belonging to unstable LP pairs. The curves for the conjugate stable modes just look like the mirror images around the central index. Each vector representing a T or LP-mode has significant contributions to the scalar product only for covariant vectors belonging to the same degenerate exponent and – to a lesser extent – the corresponding conjugate (negative) exponents (where the latter is not true anymore for the LP(3,0) modes in the bottom-right panel of Fig. 14, where no peak around i = 414 is discernible).

The covariant vectors belonging to transverse (or to LP) modes span covariant Oseledec subspaces ( with a dimension m( equal to 2 (respective 4). To ease the notation, we refer to them as ( in the following, where X is either T(n,0) or LP(n,0) with n ∈ {1, 2, 3}. The conjugate Oseledec subspaces, (, have the same dimension and are spanned by the respective conjugate covariant vectors. Fig. 14 shows that the covariant vectors spanning any of the subspaces ( or ( have a rather small but finite angular distance and, thus, are transversal. The Oseledec subspaces representing modes are themselves transversal to all other subspaces of the Oseledec splitting, but to a varying degree. The angular distances in tangent space are generally large except between conjugate subspaces ( and (, for which the scalar products of their spanning vectors may become surprisingly large.

To check more carefully for transversality even in this case, we show in Fig. 15 the probability distribution for the minimum angle between the conjugate subspaces ( and (, where X stands for the T and LP modes as indicated by the labels. This angle Φ is computed from the smallest principal angle between the two subspaces [49,50]. If the covariant vectors belonging to ( and ( are arranged as the column vectors of matrices and ∗, respectively, the QR decompositions  = QR and ∗ = Q∗R∗ of the latter provide matrices Q and Q∗, with which the matrix M = QQ∗ is constructed. The singular values of M are equal to the cosines of the principal angles, of which Φ is the minimum angle. Since Φ is never very small, this method works well and does not need more complicated refinements [49-51]. It is seen that all distributions are well bounded away from zero indicating transversality for the respective subspaces.

Fig. 15

Probability distributions for the minimum angle between the Oseledec subspace ( and its conjugate subspace (. Here, X ∈ {T(n,0), LP(n,0)} with n = 1, 2, 3 specifies the modes as indicated by the labels. The probability distribution Ψ for the minimum angle between the stable and unstable manifolds and is shown by the red line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Finally, we concentrate on the minimum angle between the full unstable subspace  = 1 ⊕ ⋯ ⊕ 393 and its conjugate stable counterpart  = 400 ⊕ ⋯ ⊕ 792, using the same method as before. These subspaces include the mode-carrying vectors studied before. The probability distribution for the minimum angle is denoted by Ψ and is also shown in Fig. 15 (red line). Also this distribution is well bounded away from zero and indicates transversality between and . We conclude that for finite N the hard-disk systems are (partially) hyperbolic in phase space.

In Fig. 13 it was observed that the scalar products between covariant vectors with adjacent indices are rather large and possibly may allow tangencies. To study this point more carefully, we follow a suggestion of G. Morriss and consider the angle Θ = cos−1∣·∣ between the vectors  ∈  ( and  ∈ (, for which i − j is a specified positive integer. The probability distributions for angles with i − j = 1, 2, 4, 20, 50 are shown in Fig. 16, and for i − j = 100, 200, 300 in the inset of the same figure. Whereas the probabilities for i − j > 1 are bounded away from zero, the distribution for i − j = 1 seems to converge to zero for Θ → 0.

Fig. 16

Probability distributions for the angles Θ between all covariant vectors and from (370) = 1 ⊕ ⋯ ⊕ 370 with prescribed separation i − j of their indices as indicated by the labels.

An even more demanding test is given in Fig. 17, where the minimum of Θ for given i − j > 0 is plotted as a function of i − j. The inset provides a magnification of the most interesting region. One observes that the minimum of the angle Θ between covariant vectors specifying Oseledec subspaces with i − j = 1 may indeed become very small, but this happens with extremely small probability. Our numerical evidence is consistent with the assumption that the angle becomes zero with vanishing probability.

Fig. 17

Plot of the minimum angle between different covariant vectors and from the unstable subspace without the mode-affected directions, (370) = 1 ⊕ ⋯ ⊕ 370, as a function of their index difference i − j.

Conclusion

A comparison of the covariant vectors with corresponding orthonormal Gram–Schmidt vectors reveal similarities, but also significant differences. The vectors associated with the maximum Lyapunov exponent are identical, 1 = 1, and also the leading vectors in the central manifold agree, 2 = 2. All the other corresponding vectors generally point into different tangent-space directions. Whereas the GS vectors are pairwise orthogonal by construction, the covariant vectors are not. Most notably, the perturbation contributions from the particles’ positions and momenta are significantly different and even exhibit a different symmetry between vectors from the stable and unstable manifold as in Fig. 6. For the covariant vectors these contributions agree in accordance with the time-reversal symmetry required for them, whereas for the Gram–Schmidt vectors these contributions are interchanged. Another significant difference is the degree of localization in physical space for the non-degenerate perturbations. As Fig. 5 shows, the covariant vectors are much more localized than the GS vectors in accordance with the fact that they are not dynamically constrained by re-orthogonalization.

From a theoretical point of view, an interesting result is that no tangencies occur between the respective unstable and stable manifolds and . In Fig. 15 the probability distribution Ψ of the minimum angle between stable and unstable subspaces (including the Lyapunov modes) is well bounded away from zero, and even more so for the vectors belonging to unstable respective stable modes. Thus, a hard-disk system with N = 198 particles as in our case is (partial) hyperbolic for all points in phase space. We even find that all Oseledec subspaces are pairwise transversal with non-vanishing angles between them.

We speculate that for N → ∞ the distribution for the minimum angle between and may possibly reach the origin in Fig. 15. To clarify this point further studies are required [34].

The concept of hyperbolicity is closely linked with the notion of dominated Oseledec splitting for all phase-space points [46]. We may rewrite Eq. (6) for the Lyapunov exponents, expressed in terms of the covariant vectors, according towhere t ≡ nτ, and τ is the short time interval between consecutive re-normalizations of the covariant vectors. Here, λ is expressed as a time average of a quantitywhich is referred to as local (or time-dependent) Lyapunov exponent, and is a function of the instantaneous phase point Γ(t). The Oseledec splitting is said to be dominated, if the local Lyapunov exponents, when averaged over a finite time Δ, do not change their order in the spectrum for any Δ larger than some finite Δ0 > 0 [34]. This is a very strong condition on the fluctuations of the local exponents [52,53]. For symplectic systems it is known that the domination of the splitting implies that the system is (partially) hyperbolic [46]. But it is not clear whether the converse is true in our case. The discussion of this point is deferred to a forthcoming publication [34].

The number and the dimension of the Oseledec subspaces are constant in phase space. There is no entanglement of subspaces, which has been identified as one of the main reasons for the occurrence of well-established Lyapunov modes [52]. We refer to Ref. [53] for a discussion of a simple but physically-relevant model, for which the dimensions of the stable and unstable manifolds frequently change along the trajectory.

An interesting extension of this work is the study of rough hard particles allowing for energy exchange between translational and rotational degrees of freedom [34]. Arguably, this is the simplest model of a molecular fluid. No Lyapunov modes are found in this case [54]. An analysis in terms of covariant vectors is presently under way and will be published separately.