# 1 Introduction

## 1 Introduction

The Lyapunov spectrum is recognized as an important diagnostic of chaotic dynamical systems. As such, it has been studied intensely in the context of extended systems . It has been observed that in the thermodynamic limit the spectrum seems to approach a continuous density, and some theoretical studies have focused on this phenomenon . However, despite the large amount of available data, there is an unsatisfactory degree of understanding of the relation between the Lyapunov spectrum of extended systems, and their global or collective properties.

In connection with these problems, the recent study of  presents an interesting development (see also  for further results and references). In the context of molecular dynamics simulations, they find hydrodynamic, i.e., slow, long-wavelength behavior in the tangent space dynamics. Namely, they observe that the Lyapunov vectors associated with the Lyapunov exponents of small absolute value have ordered, wave-like structure, and that the exponents themselves follow an ordered pattern. Hydrodynamic behavior in phase space is of course present in every extended system with a continuous symmetry. In the models considered in  the symmetries in question are translation and Galilei invariance, precisely those which enable the hydrodynamic description of a fluid in terms of the Navier-Stokes equations . However, it is for the first time that a similar phenomenon is observed in tangent space.

In this paper, we study theoretically the slow Lyapunov modes (vectors and exponents) of extended systems with translation invariance. We focus attention on a simplified model which shares the essential features with the more elaborate model of . This simplified model is constructed only in tangent space without an accompanying real space dynamics, and is based on a random matrix approximation. As has been often found before, in systems with strong chaos, qualitative features of the Lyapunov spectrum are well reproduced by approximating the tangent matrices by independent random matrices with appropriately chosen distributions . We prove several statements on the slow Lyapunov modes of this model in the thermodynamic limit, which show that in this limit the Lyapunov vectors and exponents are indeed well described as being hydrodynamic.

The basic reason for the existence of these hydrodynamic modes is evidently the translation invariance. Its presence dictates that the dynamics are indifferent to a uniform shift of all the particles (or their momenta), so that the associated Lyapunov vectors are decoupled from the rest of the dynamics, and the associated Lyapunov exponents vanish. We show that slowly growing large wavelength disturbances are nearly decoupled for the same reason, and use this property to show how the clean wave structure is obtained as a result of the orthogonalization procedure which involves all the faster growing Lyapunov vectors. It should be emphasized that the wave-like structure characterizes the Lyapunov vector at any given instant and is not an average property. Our arguments depend essentially on the local and hyperbolic character of the interactions, in addition to translation invariance. The absence of translation invariance has been recognized to ruin the hydrodynamic modes . In translation invariant anharmonic chains, the absence of short time hyperbolicity seems to ruin the hydrodynamic modes . We present theoretical arguments for the existence of hydrodynamic modes in the simplified model, which are complemented by numerical verifications. The outline of the paper is as follows. In section 2 the hydrodynamic phenomenology is described in some more detail, the definition of the random matrix model is presented and motivated, and the results are stated. They are derived in sections 35. Section 6 is devoted to numerical studies of the hydrodynamic properties of Lyapunov vectors.

## 2 Hydrodynamic behavior in tangent space

The systems studied in  consist (among others) of a large set of disks moving in a two dimensional box $`\mathrm{\Omega }`$ with periodic boundary conditions (torus geometry), with elastic scattering. In this case, the phase space is $`4N`$-dimensional where $`N`$ is the number of disks. The Lyapunov vectors have $`4N`$ components which we label as

$$(\delta x_n,\delta y_n,\delta p_{x,n},\delta p_{y,n})1nN,$$

with evident notation. To give them a geometrical meaning the components of the Lyapunov vectors are drawn in  at the instantaneous position of the particles which carry a given specific index. That is, one constructs a vector field $`\stackrel{}{v}(t,\stackrel{}{x})`$ with values in $`𝐑^4`$, which are defined only at the instantaneous positions $`\stackrel{}{x}_n(t)`$ of the particles, for example

$$v_x(t,\stackrel{}{x}_n(t))=\delta x_n(t),$$

and similarly for the other components.

The vector fields of the slow Lyapunov vectors, defined as the Lyapunov vectors with small corresponding Lyapunov exponents, are very well approximated by the long wavelength eigenmodes of a ‘reverse wave equation’ in the domain $`\mathrm{\Omega }`$:

$$_t^2\stackrel{}{v}(t,\stackrel{}{x})=\frac{1}{N^2}^2\stackrel{}{v}(t,\stackrel{}{x})$$
(1)   
(note the unusual sign in front of $`^2`$). That is, the vectors look like long wavelength waves with, say, $`n`$ nodes in the $`x`$ direction and $`m`$ nodes in the $`y`$ direction, and the corresponding Lyapunov exponent is proportional to

$$\pm \frac{1}{N}\sqrt{\left(\frac{m}{L_x}\right)^2+\left(\frac{n}{L_y}\right)^2}.$$

Note that the translation modes—constant Lyapunov vectors with zero exponents, which are trivially present in any system with translation invariance, correspond to the special case $`m=n=0`$. This phenomenology was observed in simulations with widely varying parameters, such as aspect ratio, density, and the shape of the particles .

The tangent flow of the molecular dynamics system can be written as

$$_t\left(\begin{array}{c}\delta \stackrel{}{x}\\ \delta \stackrel{}{p}\end{array}\right)=G(t)\left(\begin{array}{c}\delta \stackrel{}{x}\\ \delta \stackrel{}{p}\end{array}\right),$$
(2)   
where the components $`\delta \stackrel{}{x},\delta \stackrel{}{p}`$ are column vectors with $`N`$ entries each of which is a vector in $`𝐑^2`$. The quantity $`G(t)`$ is the action on the tangent space induced by the flow $`\mathrm{\Phi }(t)`$ of the dynamical system: If $`\psi _0`$ is the instantaneous state of the system then $`G(t)`$ is given by $`G(t)f=D\mathrm{\Phi }(t)_{\psi _0}f`$, where $`\mathrm{\Phi }(t)(\psi _0+\epsilon f)=\mathrm{\Phi }(t)\psi _0+\epsilon D\mathrm{\Phi }(t)_{\psi _0}f+O(\epsilon ^2)`$. The evolution operator of eq. (2) may be written formally as

$$U=𝖳\mathrm{exp}G(t)𝑑t.$$
(3)   
We proceed to construct a simplified model of the tangent dynamics, by making a series of modifications and assumptions about the nature of $`G`$ and $`U`$.

We first replace the hard-core interaction with a short range ‘soft’ potential. In that case $`G`$ will have a block structure of the form

$$G(t)=\left(\begin{array}{cc}0& 1\\ A(t)& 0\end{array}\right),$$
(4)   
where the symmetric $`N\times N`$ matrix $`A`$ depends on the instantaneous positions of the particles and couples only nearest neighbors. Since the interactions are purely repulsive, the flow is hyperbolic, which implies that $`A`$ should be taken positive.

At this stage one may note that if the matrix $`A(t)`$ in eq. (4) were replaced by the negative of the discrete Laplacian, the Lyapunov spectrum of $`G`$ would be precisely that described in eq. (1). However, the matrices $`A(t)`$ are in fact generated by chaotic dynamics, and therefore fluctuate rapidly. Furthermore, the particles in the gas rearrange in time, so that the positions of the non-zero elements in the matrix also evolve. In our study we concentrate on the first feature. That is, we show how tangent dynamics of the type (4) result in slow hydrodynamic modes in spite of the fluctuations in $`A(t)`$; the effects of particle rearrangement may in principle dealt with similarly, but need to be studied further.

The above discussion allows us to conclude that the matrices $`A(t)`$ should have non-zero elements only at those positions which are nonzero in the discrete Laplacian. Furthermore, momentum conservation implies that the sum of elements in any row and column of $`A`$ must vanish. This specifies completely the matrix structure of $`A`$, and it remains to model the time dependence of the off-diagonal non-zero elements of $`A`$. For this we invoke the hypothesis of strong chaos : The elements of $`A`$ may be treated as independent random processes, with a correlation time $`\tau `$ which is short with respect to other time scales of the system. It is commonly found that this approximation yields results which are in good qualitative agreement with those of the actual tangent flow.

With this in mind we model the evolution operator $`U`$ by a product of independent random matrices $`S_n`$

$$U=\underset{n}{}S_n,$$
(5)   
where

$$S_n𝖳\mathrm{exp}_{(n1)\tau }^{n\tau }G(t)𝑑t.$$
(6)   
During the time interval of length $`\tau `$, $`A`$, and therefore $`G`$, may be considered constant, so the simplest model for $`S`$ would be $`S=\left(\begin{array}{cc}1& \tau \\ \tau A& 1\end{array}\right)`$. However it is more convenient to correct this form by a second order term in $`\tau `$ in order to preserve the symplectic property which holds for $`U`$. We thus arrive at our model :

$$S=\left(\begin{array}{cc}1& \tau \\ \tau A& 1+\tau ^2A\end{array}\right).$$
(7)   
The matrices $`A`$ are independent, and their off-diagonal elements are independent and identically distributed. The actual probability distribution of the off-diagonal elements can be chosen arbitrarily, subject to the constraint of uniform hyperbolicity, namely that the support of the distribution is strictly negative, and bounded away from zero.

The model as defined above makes sense in any space dimension, but for the sake of simplicity we study it in one dimension. There it bears similarity to the tangent dynamics of an anharmonic chain. However, in the latter case the matrices $`S`$ would be elliptic rather than hyperbolic. As we show below this is an essential ingredient in the mechanism for hydrodynamic modes, which are not present in the Lyapunov spectra of anharmonic chains . Unfortunately, we have not been able to find a model dynamics in real space whose tangent space dynamics would resemble that generated by the matrices of type (7). On the other hand our results do not use explicitly the dimensionality of the system and seem to be generalizable to higher dimensions.

Since the individual matrices are symplectic, the Lyapunov exponents of (5) come in pairs of equal absolute value and opposite signs. Translation and Galilei invariance imply the existence of two vanishing exponents. We concentrate our attention on the Lyapunov exponents $`\lambda _{N1}`$ and $`\lambda _{N2}`$ of smallest positive value, and the corresponding Lyapunov vectors $`v_{N1}`$ and $`v_{N2}`$.

Before we go on, it is necessary to make precise what we mean by Lyapunov vectors. As is well-known, essentially all numerical methods for calculating the tangent space dynamics of the the Lyapunov vectors $`v_n`$ are defined as follows: One starts with an orthogonal matrix $`Q`$, mutiplies it from the left with the tangent matrix (in the present case $`S`$) and decomposes the result as $`SQ=Q^{}R`$, where $`Q^{}`$ is orthogonal and $`R`$ is upper triangular. This procedur is iterated to yield a sequence of $`Q_t`$. The columns of the orthogonal matrices $`Q_t`$ are what we will call the Lyapunov vectors. The reader should note that these vectors are not the ones whose existence is proved in the multiplicative ergodic theorem.

We can now state the main result of this paper. Existence of hydrodynamic Lyapunov modes: As the size $`N`$ of the matrices tends to infinity the exponents $`\lambda _{N1}`$ and $`\lambda _{N2}`$ as well as the vectors $`v_{N1}`$ and $`v_{N2}`$ are asymptotic to the exponents and vectors that would be obtained if the matrices $`A`$ where replaced everywhere by the negative of the discrete Laplacian (properly rescaled). The statement holds for the Lyapunov vectors, which are random objects, in probability.

The statement is spelled out only for two Lyapunov modes, which have nearly equal exponents, and where the deviation from hydrodynamic behavior is the smallest. However, as will become apparent from the arguments below, the result can be extended to a number of Lyapunov modes near the middle of the spectrum which is proportional to $`\sqrt{N}`$. Our numerical studies also indicate that this is in fact true.

As already explained, the basic reason for the existence of the hydrodynamic modes is translation invariance. However, this general observation is not sufficient, and the actual proof is not trivial. It depends on the random nature of successive matrices, i.e., on strong chaos. Our strategy will be to show the existence of hydrodynamic modes first in the spectrum of a single matrix of type $`A`$, i.e., a negative random Laplacian (sec. 3); then this will be used to show that such modes exist in the Lyapunov spectrum of non-symplectic products of type $`(1+A_n)`$ in sec. 4, which in turn will be used to show the same property for symplectic products in sec. 5.

## 3 Spectral properties of a single matrix

The matrix $`A`$ defined in section 2 takes in one dimension the explicit form

$$A=\left(\begin{array}{ccccccc}a_1+a_2& a_2& 0& 0& \mathrm{}& 0& a_1\\ a_2& a_2+a_3& a_3& 0& \mathrm{}& 0& 0\\ \mathrm{}& & & & & & \mathrm{}\\ a_1& 0& 0& 0& \mathrm{}& a_N& a_N+a_1\end{array}\right),$$
(8)   
where the $`a_n`$ are positive identically distributed independent random numbers. The distribution of the $`a_n`$ is arbitrary, subject to the condition $`0<a_{\mathrm{min}}<a<a_{\mathrm{max}}`$ with $`a_{\mathrm{min}}<a_{\mathrm{max}}`$, and normalized such that $`a^1=1`$, for later convenience. The $``$ always denote expectation with respect to the probability distribution of the $`a`$. We are not going to assume that the width of the distribution is small.

The matrix $`A`$ may be written as a product

$$A=\underset{¯}{}𝒜\overline{}$$
(9)   
where $`\underset{¯}{}`$ and $`\overline{}`$ are the discrete derivatives whose action on a vector $`v𝐑_N`$ is

$$(\underset{¯}{}v)_n=v_nv_{n1},(\overline{}v)_n=v_{n+1}v_n,$$
(10)   
and $`𝒜`$ is a diagonal matrix with diagonal elements $`a_n`$. (The indices are extended periodically so that $`a_{N+1}a_1`$). If all the $`a_n`$ were equal to one, $`A`$ would reduce to the discrete Laplacian matrix $`^2\underset{¯}{}\overline{}`$.

We define the Fourier transform matrix $`F`$ with elements

$$F_{kn}=\frac{1}{\sqrt{N}}e^{i\frac{2\pi }{N}kn},$$
(11)   
which is a unitary transformation taking $`𝐑^N`$ to $`\stackrel{~}{𝐂}^N`$, the subset of $`𝐂_N`$ (with standard basis vectors $`e_n`$) consisting of vectors $`\stackrel{~}{v}`$ for which $`v_k=v_k^{}`$, which is an $`N`$ dimensional vector space over $`𝐑`$. The components of $`A`$ in the new basis are

$$\stackrel{~}{A}_{kl}(FAF^{})_{kl}=\mu _k^{}(a\delta _{k,l}+\stackrel{~}{a}_{kl})\mu _l,$$
(12)   
where $`\mu _k=1\mathrm{exp}\left(\frac{2\pi i}{N}k\right)`$, and $`\stackrel{~}{a}`$ is related to the Fourier transform of $`a`$ considered as a vector in $`𝐑^N`$ by

$$\stackrel{~}{a}=N^{1/2}F(\stackrel{}{a}\stackrel{}{a}).$$
(13)   
The random variables $`\stackrel{~}{a}_k`$ are centered, and as sums of independent random numbers their ‘single-point’ distribution is nearly Gaussian with variance

$$|\stackrel{~}{a}|^2=\frac{a^2a^2}{N},$$
(14)   
so that they are typically small, of order $`𝒪(N^{1/2})`$. The joint distribution is not Gaussian.

Note that $`\mu _0=0`$, so that row $`0`$ and column $`0`$ of $`\stackrel{~}{A}`$ are zero, with the translation vector $`e_0`$ being trivially a zero eigenvector. We define the slow subspace

$$V_\mathrm{s}=\mathrm{Span}(\{e_0,e_1,e_1\})\stackrel{~}{𝐂}^N,$$
(15)   
and its orthogonal complement, the fast subspace $`V_\mathrm{f}`$. We will consider often below the block decomposition of $`\stackrel{~}{A}`$ and other matrices into the fast and slow subspaces, e.g.,

$$\stackrel{~}{A}=\left(\begin{array}{cc}A_{\mathrm{ff}}& A_{\mathrm{fs}}\\ A_{\mathrm{sf}}& A_{\mathrm{ss}}\end{array}\right).$$
(16)   
Note that $`V_\mathrm{f}`$ contains slow as well as fast modes.

The block $`A_{\mathrm{ss}}`$ has small norm of order $`𝒪(N^2)`$, and the off-diagonal blocks have norm of order $`𝒪(N^1)`$. However, there are more specific properties of $`A`$ which are needed to establish the existence of hydrodynamic eigenmodes. Consider the eigenvalue problem $`Av=\lambda v`$. Letting $`v=\underset{¯}{}u`$, and using the representation (9) gives an equation for $`u`$

$$^2u=\lambda 𝒜^1u.$$
(17)   
It is convenient to proceed by writing eq. (17) in Fourier component form

$$(|\mu _k|^2\lambda )\stackrel{~}{u}_k=\lambda \underset{q}{}\stackrel{~}{b}_{kq}\stackrel{~}{u}_q,$$
(18)   
with the $`\stackrel{~}{b}_k`$ bearing a relation to $`a_n^1`$ analogous to that between $`\stackrel{~}{a}_k`$ and $`a_n`$, namely

$$\stackrel{~}{b}_k=N^{1/2}F_{kn}(\frac{1}{a_n}1).$$
(19)   
Since $`\frac{1}{a^2}<1/a_{\mathrm{min}}^2`$ is $`𝒪(1)`$ we find that the $`\stackrel{~}{b}_k`$ are $`𝒪(N^{1/2})`$ for the same reason that the $`\stackrel{~}{a}_k`$ are.

We claim that given a fixed $`m`$, and for $`N\mathrm{}`$ the system (18) has two linearly independent solutions $`u^{(\pm m)},\lambda _{\pm m}`$ such that

$$\frac{1}{|u_m^{(\pm m)}|}\underset{|k|m}{}|u_k^{(\pm m)}|=𝒪(N^{1/2}),$$
(20)   
and

$$\left|\frac{\lambda _{\pm m}}{|\mu _m|^2}1\right|=𝒪(N^{1/2}).$$
(21)   
We justify the claim by showing that eqs. (20) and (21) are consistent with the eigenvalue equation (18). For this we rewrite (18) as

$$\stackrel{~}{u}_k^{(m)}=\frac{\lambda _m}{|\mu _k|^2\lambda _m}\underset{q}{}\stackrel{~}{b}_{kq}\stackrel{~}{u}_q^{(m)}.$$
(22)   
We assume that (2021) hold; this implies that the sum over $`q`$ in (22) is dominated by the two terms with $`q=\pm m`$, that is,

$$\stackrel{~}{u}_k^{(m)}=\frac{|\mu _m|^2}{|\mu _k|^2|\mu _m|^2}(b_{km}\stackrel{~}{u}_m^{(m)}+b_{k+m}\stackrel{~}{u}_m^{(m)}),\text{for }|k|m\text{.}$$
(23)   
On substituting this expression in the left-hand-side of (20) the sum over $`k`$ is observed to be local, in the sense that it is dominated by terms with $`|k|m`$, where $`|\mu _k|^2\left(\frac{2\pi k}{N}\right)^2`$. Since $`\stackrel{~}{b}_k`$ is $`𝒪(N^{1/2})`$, assumption (20) is verified. On the other hand, using (20) in (22) for $`k=m`$ gives

$$\stackrel{~}{u}_m^{(m)}=\frac{\lambda _m}{|\mu _m|^2\lambda _m}(b_0\stackrel{~}{u}_m^{(m)}+b_{2m}\stackrel{~}{u}_m^{(m)}).$$
(24)   
Since $`b_0`$ and $`\stackrel{~}{b}_{2m}`$ are $`O(N^{1/2})`$ it follows that $`\lambda _m/(|\mu _m|^2\lambda _m)=O(N^{1/2})`$ verifying (21), which shows that (2021) are indeed consistent with (18).

In terms of the original variables $`v`$, eq. (23) reads

$$\stackrel{~}{v}_k^{(m)}=\frac{\mu _k^{}\mu _m}{|\mu _k|^2|\mu _m|^2}(b_{km}\stackrel{~}{v}_m^{(m)}+b_{k+m}\stackrel{~}{v}_m^{(m)}),\text{for }|k|m\text{,}$$
(25)   
so that the norm of $`v_{}^{(m)}`$, the component of $`v^{(m)}`$ orthogonal to $`e_{\pm m}`$ is small,

$$v_{}^{(m)}^2\frac{1}{N}\underset{|k|m}{}\frac{k^2m^2}{(k^2m^2)^2}=O\left(\frac{1}{N}\right).$$

In words, these eigenvectors are almost pure Fourier modes, i.e., eigenvectors of the discrete Laplacian.

For further developments we also need to show these modes are the only ones with eigenvalues of order $`𝒪(N^2)`$.

This is established easily by noting that the sharp cutoff on the probability distribution of the $`a`$ implies that every realization $`A`$ satisfies the bounds

$$a_{\mathrm{min}}^2<A<a_{\mathrm{max}}^2,$$
(26)   
and then by the minimax principle it follows that the $`p`$th eigenvalue of $`A`$ is larger than $`a_{\mathrm{min}}`$ times the $`p`$th eigenvalue of $`^2`$ (sorting the eigenvalues of both matrices in increasing order).

The results of this section can be summarized using the decomposition of $`\stackrel{~}{𝐂}^N`$ into slow and fast subspaces defined above. We have shown there exist small numbers $`\epsilon `$ and $`\lambda `$, and a number $`0<\alpha <1`$, such that the matrix $`\stackrel{~}{A}`$ can be block diagonalized,

$$\stackrel{~}{A}=RDR^T,RR^T=1,D=\left(\begin{array}{cc}D_\mathrm{f}& 0\\ 0& D_\mathrm{s}\end{array}\right),$$
(27)   
with the off-diagonal blocks bounded by $`R_{\mathrm{sf}},R_{\mathrm{fs}}<\epsilon `$, and the diagonal blocks obeying

$$D_\mathrm{f}>\lambda >\alpha \lambda >D_\mathrm{s}0,$$
(28)   
and furthermore

$$A_{\mathrm{ss}}<\alpha \lambda .$$
(29)   
The orders of magnitude for $`\lambda `$ and $`\epsilon `$ are $`\epsilon =𝒪(N^{1/2})`$ and $`\lambda =𝒪(N^2)`$, whereas $`\alpha 1/4`$. However, to keep the discussion reasonably general we are not going to use these specific values in our arguments below. Rather, we will make statements regarding arbitrary matrices which satisfy the conditions (2729).

Although this will not be used below, it is relevant to note that if we let $`m`$ vary, the small parameter in (20) and (21) becomes $`m/N^{1/2}`$. This means that we can expect a number of hydrodynamic eigenmodes which is proportional to $`\sqrt{N}`$. Another way to see this is related to the study of the vibrations of one-dimensional disordered lattices which are modeled precisely by the eigenmodes of matrices of type (8). There it is known that the localization length $`\xi `$ is proportional to $`\lambda _m^1`$. Since $`\lambda _m(2\pi /N)^2`$ the localization length will reach $`N`$ when $`m=𝒪(N^{1/2})`$. Thus, again, we only expect wave-like modes when $`m<𝒪(N^{1/2})`$.

## 4 Products of matrices of the form $`1+\tau A`$

In this section we use the properties derived in section 3 to derive the existence of hydrodynamic modes in the Lyapunov spectrum of the product $`_n(1+\tau A_n)`$, where the matrices $`A_n`$ are independent realizations of the random matrix defined in eq. (8). Beside providing a step towards proving the existence of hydrodynamic modes in symplectic products, such a product may be regarded as the discrete approximation to a continuous tangent flow given by

$$U=𝖳\mathrm{exp}A(t)𝑑t.$$
(30)   
\[Compare eqs. (3) and (4).\] Although this does not correspond to the tangent flow of a mechanical system, it is nonetheless the simplest example where hydrodynamic Lyapunov modes can be expected. For convenience of further analysis we absorb $`\tau `$ into the definition of $`A`$ and change to Fourier basis once and for all, so the problem becomes that of a product

$$\underset{n}{}(1+\stackrel{~}{A}_n).$$
(31)   
Since the Lyapunov exponents of the slow part are expected to be smaller than the rest we aim at showing that the first $`N3`$ Lyapunov vectors span a subspace $`L_\mathrm{f}`$ (of $`\stackrel{~}{𝐂}_N`$) which is almost orthogonal to $`V_\mathrm{s}`$ in the sense that for any two unit vectors $`u_\mathrm{f}L_\mathrm{f}`$ and $`v_\mathrm{s}V_\mathrm{s}`$ one has $`|u_\mathrm{f}v_\mathrm{s}|1`$. Although the subspace $`L_\mathrm{f}`$ changes after each step, we show below that the ‘almost orthogonality’ is propagated from step to step.

To show this, we propose the following scheme. Take an arbitrary vector $`uL_\mathrm{f}`$, whose components in $`V_\mathrm{f}`$ and $`V_\mathrm{s}`$ are $`u_\mathrm{f}`$ and $`u_\mathrm{s}`$ respectively, normalized so that $`u_\mathrm{f}=1`$, and assume that $`u_\mathrm{s}`$ is small. The action of $`1+A`$ generates a new normalized vector $`u^{}`$ by

$$u^{}=\frac{(1+A)u}{[(1+A)u]_\mathrm{f}},$$
(32)   
where $`[]_\mathrm{f}`$ is the projection onto the $`\mathrm{f}`$-component We would like to show that $`u_\mathrm{s}`$ remains small after repeated iteration of this process.

The block diagonalization (27) shows that the subspaces $`V_\mathrm{f}`$ and $`V_\mathrm{s}`$ are indeed almost invariant under the transformation $`\stackrel{~}{A}`$. However, in trying to apply this fact to the Lyapunov vectors of the product (31) we immediately encounter the danger that the small perturbations may accumulate. The basic problem is that a vector in $`V_\mathrm{s}`$ is contracted with respect to the ‘slowest’ direction in $`V_\mathrm{f}`$ by a factor of only $`1𝒪(\lambda )`$ (as can be seen from the bounds on $`D_\mathrm{f}`$ and $`D_\mathrm{s}`$), whereas the perturbations which tilt a vector in $`L_\mathrm{f}`$ with respect to $`V_\mathrm{f}`$ are of order $`\epsilon `$, which is the typical size of the off-diagonal blocks \[cf. eqs. (27)–(29)\] and since we are interested in the case $`\lambda \epsilon `$ this contraction is not strong enough to overcome the perturbation.

This order of magnitude argument can be made explicit by constructing a series of matrices with the properties given by eqs. (27)–(29), which take a vector in $`V_\mathrm{f}`$ and rotate it such that the outcome is a vector which has an angle with $`V_\mathrm{f}`$ of order 1. This counter-example is given in appendix A.

An essential ingredient in the construction of this counter-example is that the matrix $`\stackrel{~}{A}`$ has to be chosen specifically given $`u`$ which in turn depends on former realizations, in violation of the independence assumption. In other words, although such a ‘bad’ sequence is possible one naturally expects that this is an event with very low probability. Typically the perturbations to $`u_\mathrm{s}`$ generated by the off-diagonal part of the matrices $`R`$ do not have the same direction, and should serve to cancel one another. Therefore the statement one can hope to show is that in the sequence generated by iteration of eq. (32), the probability that $`u_\mathrm{s}>C\epsilon `$ for some fixed $`C`$ is very small, as was shown for a similar example in . Here we will only prove the weaker statement that the variance $`u_\mathrm{s}^2`$ is $`𝒪(\epsilon ^2)`$, and take that as an indication that the probabilistic statement is correct, since the behavior of higher moments can be treated in an analogous manner.

To prove this statement we look at the $`s`$-component of eq. (32),

$$u_\mathrm{s}^{}=\frac{A_{\mathrm{sf}}u_\mathrm{f}+(1+A_{\mathrm{ss}})u_\mathrm{s}}{(1+A_{\mathrm{ff}})u_\mathrm{f}+A_{\mathrm{fs}}u_\mathrm{s}}.$$
(33)   
The quantity $`u_\mathrm{s}^{}^2`$ is a sum of three terms $`E_1+E_2+E_3`$:

$$E_1=\frac{A_{\mathrm{sf}}u_\mathrm{f}^2}{\mathrm{}^2},E_2=2\frac{A_{\mathrm{sf}}u_\mathrm{f}(1+A_{\mathrm{ss}})u_\mathrm{s}}{\mathrm{}^2},E_3=\frac{(1+A_{\mathrm{ss}})u_\mathrm{s}^2}{\mathrm{}^2},$$
(34)   
where $`\mathrm{}=(1+A_{\mathrm{ff}})u_\mathrm{f}+A_{\mathrm{fs}}u_\mathrm{s}`$.

To bound these terms we first need a lower bound on the denominator $`\mathrm{}`$. Let $`v_\mathrm{f}=R_{\mathrm{ff}}^Tu_\mathrm{f}+R_{\mathrm{fs}}^Tu_\mathrm{s}`$, and define $`d`$ by

$$D_\mathrm{f}v_\mathrm{f}dv_\mathrm{f}.$$
(35)   
Note that $`d`$ can vary widely between $`𝒪(1)`$ values and $`𝒪(\lambda )`$. But, using the lower bound on $`D_\mathrm{f}`$ of (28), we see that

$$(1+D_\mathrm{f})v_\mathrm{f}^2=v_\mathrm{f}^2+2v_\mathrm{f}D_\mathrm{f}v_\mathrm{f}+D_\mathrm{f}v_\mathrm{f}^2(1+2\lambda +d^2)v_\mathrm{f}^2.$$
(36)   
Expanding $`\mathrm{}`$ as

$$\mathrm{}=R_{\mathrm{ff}}(1+D_\mathrm{f})v_\mathrm{f}+R_{\mathrm{fs}}D_\mathrm{s}v_\mathrm{s},$$
(37)   
and using the estimates $`R_{\mathrm{fs}}D_\mathrm{s}=𝒪(\epsilon \lambda )`$ and $`1R_{\mathrm{ff}}=𝒪(\epsilon ^2)`$ (cf. eq. (28)), we get from (36) the desired lower bound on $`\mathrm{}`$:

$$\mathrm{}^2>(1+2\lambda +d^2)v_\mathrm{f}^2(1𝒪(\epsilon )).$$
(38)   
We can now bound $`E_1`$, $`E_2`$ and $`E_3`$. First, we have

$$A_{\mathrm{sf}}u_\mathrm{f}=R_{\mathrm{sf}}D_\mathrm{f}v_\mathrm{f}+𝒪(\lambda \epsilon )<\epsilon dv_\mathrm{f}+𝒪(\lambda \epsilon ).$$
(39)   
Thus, neglecting higher order corrections in $`\epsilon `$, we get

$$E_1<\frac{\epsilon ^2d^2}{1+2\lambda +d^2}.$$
(40)   
The bound on the term $`E_2`$ makes essential use of the translation invariance. For this, we note that

$$\frac{A_{\mathrm{fs}}u_\mathrm{s}}{(1+A_{\mathrm{ff}})u_\mathrm{f}^2}$$
(41)   
transforms as a vector, that is, its $`k`$th component is multiplied by $`\mathrm{exp}(i\frac{2\pi }{N}kx)`$ under a relabeling of the coordinates $`nn+x`$. Therefore, because of translation invariance, the expectation value of (41) must remain invariant under such transformation, which means it must vanish. Since the denominator in $`E_2`$ is $`\mathrm{}^2`$ (which also depends on $`u_\mathrm{s}`$) and not $`(1+A_{\mathrm{ff}})u_\mathrm{f}^2`$, we need some gymnastics to exhibit the vanishing term. In order to see this we write

$$\begin{array}{c}E_2=2<\frac{A_{\mathrm{sf}}u_\mathrm{f}u_\mathrm{s}}{(1+A_{\mathrm{ff}})u_\mathrm{f}^2}+\frac{A_{\mathrm{sf}}u_\mathrm{f}A_{\mathrm{ss}}u_\mathrm{s}}{(1+A_{\mathrm{ff}})u_\mathrm{f}^2}\hfill \\ \frac{A_{\mathrm{sf}}u_\mathrm{f}(1+A_{\mathrm{ss}})u_\mathrm{s}\left[2(1+A_{\mathrm{ff}})u_\mathrm{f}A_{\mathrm{fs}}u_\mathrm{s}+A_{\mathrm{fs}}u_\mathrm{s}^2\right]}{(1+A_{\mathrm{ff}})u_\mathrm{f}^2\mathrm{}^2}>.\hfill \end{array}$$
(42)   
The first term in (42) vanishes because of translation invariance, as explained before. The second term is bounded by

$$2\alpha \lambda \epsilon u_\mathrm{s}$$
(43)   
and the dominant part of the third is

$$4\frac{\left(A_{\mathrm{sf}}u_\mathrm{f}u_\mathrm{s}\right)\left((1+A_{\mathrm{ff}})u_\mathrm{f}A_{\mathrm{fs}}u_\mathrm{s}\right)}{(1+A_{\mathrm{ff}})u_\mathrm{f}^2\mathrm{}^2}<\frac{4d\epsilon ^2u_\mathrm{s}^2}{1+2\lambda +d^2}.$$
(44)   
The last term is bounded by

$$E_3<\frac{1+2\alpha \lambda }{1+2\lambda +d^2}u_\mathrm{s}^2.$$

Collecting the bounds yields

$$u_\mathrm{s}^{}^2<\frac{d^2\epsilon ^2+(1+2\alpha \lambda +4d\epsilon ^2)u_\mathrm{s}^2}{1+2\lambda +d^2}+2\alpha \lambda \epsilon .$$
(45)   
It appears from (45) that although large perturbations are possible when $`d`$ is $`𝒪(1)`$, the contraction rate increases precisely enough to compensate this contribution. Thus if $`u_\mathrm{s}^2`$ is $`𝒪(\epsilon ^2)`$ to start with, it will stay so indefinitely.

In summary, assuming that the variance is indeed a measure of typical fluctuations, we have shown that, for $`N1`$, the subspace $`L_\mathrm{f}`$ spanned by the first $`N3`$ Lyapunov vectors of the product (31) is, with very high probability, almost orthogonal to $`V_\mathrm{s}`$. This implies that the last three Lyapunov vectors (including the translation) remain approximately in $`V_\mathrm{s}`$. This means by definition that they are hydrodynamic, in the sense that they are well approximated by eigenvectors of the discrete Laplacian. Since the action of $`\stackrel{~}{A}`$ on $`V_\mathrm{s}`$ has two eigenvalues close to $`(2\pi /N)^2`$ as shown in the previous section, it follows as a corollary that the two smallest non-trivial Lyapunov exponents have approximately this value, so that they are also hydrodynamic. This completes the demonstration.

## 5 Products of symplectic matrices

We now turn to products of matrices of the form

$$S=\left(\begin{array}{cc}1& \tau \\ \tau A& 1+\tau ^2A\end{array}\right).$$
(46)   
We disregard the two translation modes in $`S`$ for convenience and view the matrices $`S`$ as $`(2N2)\times (2N2)`$ matrices. Let us recall that since the matrices $`S`$ are symplectic and hyperbolic, the Lyapunov exponents are non-zero and come in pairs of opposite signs. We concentrate on modes number $`N2`$ and $`N1`$ which are the smallest positive ones.

We reduce the problem to an equivalent one to which the results of section 4 can be applied directly. We denote by $`L_+`$ the subspace spanned by the first $`N1`$ Lyapunov vectors. It is spanned by a set of $`N1`$ independent vectors, which we display in the form of a $`(N1)\times (2N2)`$ matrix $`𝒱`$. The $`N1`$ vectors can always be chosen in such a way that $`𝒱`$ is of the normal form

$$𝒱=\left(\begin{array}{c}\mathrm{𝟏}\\ V\end{array}\right),$$
(47)   
where both blocks are $`(N1)\times (N1)`$. Acting on $`𝒱`$ with $`S`$ gives a spanning set of the image subspace $`L_+^{}`$,

$$S𝒱=\left(\begin{array}{c}1+\tau V\\ \tau A+(1+\tau ^2\stackrel{~}{A})V\end{array}\right),$$
(48)   
and changing basis to normal form gives $`𝒱^{}=\left(\begin{array}{c}\mathrm{𝟏}\\ V^{}\end{array}\right)`$ where

$$V^{}=\tau A+\frac{V}{1+\tau V}.$$
(49)   
A convenient property of this matrix dynamical system is that if $`V`$ is symmetric to begin with, it stays so as a consequence of the symplectic property of $`S`$ .

By definition, any vector $`vL_+`$ has a block representation $`v=\left(\begin{array}{c}u,\\ Vu\end{array}\right)`$. From eq. (48) it follows that its image is $`\left(\begin{array}{c}u^{},\\ V^{}u^{}\end{array}\right)`$, where

$$u^{}=(1+\tau V)u.$$
(50)   
Hence, the first $`N1`$ Lyapunov modes of the products of the $`S`$ are the same as those of the product $`_n(1+\tau V_n)`$ where the matrices $`V_n`$ are evolving according to eq. (49): $`V_{n+1}=\tau A_n+V_n/(1+\tau V_n)`$.

In view of this equivalence, it suffices to show that the matrix $`V`$ has the properties formulated in eqs. (2729) and to apply the results of section 4. First note if $`A_n=^2`$ (minus the discrete Laplacian) for all $`n`$, then the $`V_n`$ converge to $`f(^2)`$, where

$$f(x)=\frac{\tau x}{2}+\sqrt{x+\left(\frac{\tau x}{2}\right)^2}$$

is the larger root of the quadratic equation $`f(x)=\tau x+\frac{f(x)}{1+\tau f(x)}`$. For small $`x>0`$ this is close to $`x^{1/2}`$, and therefore we assume that $`V`$ has a representation of the type given by eqs. (2729), with $`\epsilon =𝒪(N^{1/2})`$ as before and $`\lambda `$ is now $`f(4\pi /N^2)=𝒪(N^1)`$. The aim is to show that this property is carried on to $`V^{}`$.

In order to avoid the necessity of presenting even more technical details, we present the argument for the case where the slow subspace contains a single mode, rather than a pair of nearly degenerate modes. Since $`V^{}`$ is symmetric its smallest eigenvalue is given by

$$\lambda _V^{}=\underset{u=1}{\mathrm{min}}u(\tau A+\frac{V}{1+\tau V})u.$$
(51)   
It follows from our assumptions that there exist (normalized) eigenvectors of $`V`$ and $`A`$:

$$Ae_A=\lambda _1(1+c_A\epsilon )e_A,Ve_V=f(\lambda _1)(1+c_V\epsilon )e_V,$$
(52)   
where $`\lambda _1`$ is the smallest positive eigenvalue of $`^2`$ and $`e_A`$ and $`e_V`$ are close to $`e_1`$, the corresponding eigenvector. The variational principle then gives immediately a lower bound on $`\lambda _V^{}`$,

$$\lambda _V^{}>f(\lambda _1)(1+\epsilon c_V^{}),$$
(53)   
for some $`c_V^{}`$ between $`c_A`$ and $`c_V`$. To get an upper bound on $`\lambda _V^{}`$ recall that it was shown in sec. 3 eq. (25) that the $`k`$ component of $`e_A`$ is of order $`\epsilon k^1`$, and note that the bound (26) implies

$$f(a_{\mathrm{min}}^2)<V<f(a_{\mathrm{max}}^2).$$
(54)   
We now use $`u=e_A`$ in eq. (51) and get

$$\lambda _V^{}<\lambda _1(1+c_A\epsilon )+e_Af(a_{\mathrm{max}}^2)e_A<f(\lambda _1)(1+\overline{c}\epsilon ),$$
(55)   
for some constant $`\overline{c}`$. Eqs. (53) and (55) establish the desired property of the eigenvalues.

The corresponding eigenvector $`e_V^{}`$ is the one which minimizes eq. (51). Because of the minimax principle applied to $`A`$ and $`V`$, letting $`u=e_1+w`$ with $`we_1=0`$ and $`w`$ small, the quadratic form $`uAu`$ may be approximated by

$$uAu(ww_A)A(ww_A)+\lambda _A,$$

and similarly

$$uVu(ww_V)V(ww_V)+\lambda _V.$$

Therefore, in order to find $`w_V^{}`$ we need to minimize the quadratic form

$$\tau (ww_A)A(ww_A)+(ww_V)\frac{V}{1+\tau V}(ww_V).$$

The minimum occurs at

$$w_V^{}=(1+B)^1(Bw_A+w_V),$$
(56)   
where

$$B=\tau \frac{1+\tau V}{V}A.$$

We can use again the bounds (26) and (54) to show that

$$b_{\mathrm{min}}g(^2)<B<b_{\mathrm{max}}g(^2),$$

for some positive numbers $`b_{\mathrm{min}},b_{\mathrm{max}}`$ and a positive function $`g`$, and thus bound the components of $`w_V^{}`$,

$$|(w_V^{})_k|<\frac{(w_V)_k+b_{\mathrm{max}}g(|\mu _k|^2)(w_A)_k}{1+b_{\mathrm{min}}g(|\mu _k|^2)},$$

where $`|\mu _k|^2`$ is the $`k`$th eigenvalue of $`^2`$ (see sec 3). This shows that $`(w_V^{})_k\epsilon k^1`$. This completes the demonstration of the desired properties of $`V^{}`$, and, on applying the results of section 4, the existence of hydrodynamic Lyapunov modes in the symplectic case.

## 6 Numerical tests

The purpose of this section is to verify numerically some of the statements given above, and on to further study numerically the dependence of hydrodynamic behavior of several Lyapunov modes as a function of noise level as well as system size.

The simplest system we discuss is a product $`_nA_n`$ of independent matrices of the form (8). Since the relative gap between the first two non-zero eigenvalues is $`𝒪(1)`$ \[see eqs. (2729)\], the contraction in this case is strong, and the potential problems of the accumulation of errors as discussed in sec. 4 and appendix A are absent. Nevertheless, even in this case, there are some qualitative differences between the behavior of the Lyapunov modes, and the corresponding eigenmodes of a single matrix.

We quantify the degree of hydrodynamic behavior in the Lyapunov modes as follows. For the Lyapunov vectors $`v_i`$ we computed the residuals $`r_i`$, that is, the norm of the orthogonal complement

$$r_i=v_i(v_ie_k)e_k(v_ie_k)e_k$$

where $`k=k(i)`$ is the wave vector associated with vector $`i`$. (For example, $`k=1`$ for the vectors $`v_{N1}`$ and $`v_{N2}`$ discussed above.) In fact, to get more precise results we subtracted from $`v_i`$ all the components with lower-lying $`k`$:

$$r_i=v_i\underset{kk(i)}{}\left((v_ie_k)e_k(v_ie_k)e_k\right)$$

(The results are not very different for the two definitions of $`r_i`$.)

Fig. 1 presents these residuals for systems with different sizes of the matrices and different values of the noise variance $`\sigma `$. The vertical axis measures $`(r_k/\sigma )^2`$, and the horizontal axis gives $`k/N`$. The approximate collapse of the graphs for small $`k`$ implies that the dependence of the residuals on system size $`N`$ and noise strength (variance) $`\sigma `$) is given by the scaling form

$$r_k=\sigma f_1\left(\frac{k}{N}\right).$$

The behavior of $`f_1`$ for small $`x`$ is approximately $`f_1(x)=𝒪(\sqrt{x})`$, which implies that for a fixed $`k`$

$$r_k\frac{\sigma }{\sqrt{N}},$$

the same dependence as in the residuals of the eigenvectors of a single matrix. However, the dependence of the residuals as a function of $`k`$ is $`r_k\sqrt{k}`$, slower the linear dependence on $`k`$ in the case of a single matrix.

For the Lyapunov exponents $`\lambda _i`$ we measure the relative deviation $`\delta _k`$ from the respective eigenvalues $`|\mu _k|^2`$ of the discrete Laplacians (see sec. 3),

$$\delta _k=\frac{\mathrm{exp}(\lambda _k)}{|\mu _k|^2}1.$$

The results for the deviations $`\delta _i`$ of the Lyapunov exponents are displayed in fig. 2, where $`\delta _k/\sigma ^2`$ is plotted against $`k/N`$. The data collapse implies that

$$\delta _k=\sigma ^2f_2\left(\frac{k}{N}\right).$$

The function $`f_2(x)`$ is approximately linear for small $`x`$ which implies for $`kN`$:

$$\delta _kr_k^2.$$

This is not unreasonable, since the Lyapunov exponents, unlike the eigenvalues of a single matrix, are given as a result of an averaging process.

We next present a similar analysis for the product $`_n(1+A_n)`$ which was considered in section 4. The results for the residuals of the Lyapunov vectors and the deviations of the exponents are presented in figures 3 and 4 respectively. The scaling form for the residuals is in this case

$$r_k=\sigma f_3\left(\frac{k}{N}\right),$$

where $`f_3(x)x`$ for small $`x`$. Thus, in this case the residual for a fixed $`k`$ decreases as

$$r_k\frac{\sigma }{N},$$

that is, faster than the $`N^{1/2}`$ decrease in the residuals of the eigenvectors of a single matrix. The analysis presented in section 4 is too general to capture this behavior.

The relative deviations of the exponents scale in this case as

$$\delta _k=\sigma ^2f_4\left(\frac{k}{N}\right),$$

and $`f_4(x)x^2`$ for small $`x`$, so that as in the product of random Laplacians, the relative size of $`r_k`$ and $`\delta _k`$ is $`\delta _kr_k^2.`$

## Appendix A Counter-example

We want to show here that an ‘unfortunate’ choice of rotations can move the system out of the region where the Lyapunov vectors remain essentially aligned with the eigendirections of the Laplacian. The issue here is that, on one hand, the cones in which these vectors lie are slightly contracted and on the other hand slightly turned. The ‘counter-example’ shows that the turning wins over the contraction.

Let $`V_\mathrm{f}=\mathrm{Span}\{e_1,e_2\}`$ and $`V_\mathrm{s}=\mathrm{Span}\{e_3\}`$. Suppose that $`L_\mathrm{f}`$ contains a vector with block representation $`u=(u_\mathrm{f},u_\mathrm{s})`$, with $`u_\mathrm{s}>0`$, normalized such that $`u_\mathrm{f}=1`$, and let $`v_\mathrm{f}`$ span the orthogonal complement to $`u_\mathrm{f}`$ in $`V_\mathrm{f}`$. We construct the matrix $`\stackrel{~}{A}`$ by giving the components in the representation (27),

$$R\left(\begin{array}{cc}\mathrm{𝟏}& \epsilon v_\mathrm{f}\\ \epsilon v_\mathrm{f}^T& 1\end{array}\right),$$
(57)   
$`D_\mathrm{s}=\alpha \lambda `$, and $`D_\mathrm{f}`$ is such that

$$\begin{array}{c}D_\mathrm{f}u_\mathrm{f}=(\lambda +\epsilon ^2)u_\mathrm{f}+\epsilon v_\mathrm{f}\\ D_\mathrm{f}v_\mathrm{f}=\epsilon u_\mathrm{f}+v_\mathrm{f}\end{array}.$$
(58)   
The image of $`u`$ is

$$(1+\stackrel{~}{A})\left(\begin{array}{c}u_\mathrm{f}\\ u_\mathrm{s}\end{array}\right)=\left(\begin{array}{c}(1+\lambda +\epsilon ^2(1u_\mathrm{s}))u_\mathrm{f}+\epsilon v_\mathrm{f}\\ (1+\alpha \lambda \epsilon ^2)u_\mathrm{s}+\epsilon ^2\end{array}\right).$$
(59)   
After normalizing the $`f`$ component to 1, the $`s`$ component becomes

$$u_\mathrm{s}^{}\left[1+(\alpha 1)\lambda +\epsilon ^2(3/2u_\mathrm{s})\right]u_\mathrm{s}+\epsilon ^2.$$
(60)   
Evidently, even if $`u_\mathrm{s}=0`$ initially, by choosing $`\stackrel{~}{A}`$ as above, $`u_\mathrm{s}`$ can be increased to an $`𝒪(1)`$ value (as $`\epsilon ,\lambda 0`$) if $`\lambda =𝒪(\epsilon ^2)`$.

Acknowledgments: We have profited from very useful discussions with S. Lepri, C. Liverani, Z. Olami, A. Politi, and L.-S. Young. But we are particularly indebted to H. Posch for having provided and discussed with us his beautiful numerical experiments. This work was partially supported by the Fonds National Suisse, and part of it was done in the pleasant atmosphere of the ESI in Vienna and at the Istituto Nazionale di Ottica in Florence.