Lyapunov modes in soft-disk fluids

Lyapunov modes are periodic spatial perturbations of phase-space states of many-particle systems, which are associated with the small positive or negative Lyapunov exponents. Although familiar for hard-particle systems in one, two and three dimensions, they have been difficult to find for soft particles. We present computer simulations for soft-disk systems in two dimensions and demonstrate the existence of the modes, where also Fourier-transformation methods are employed. We discuss some of their properties in comparison with equivalent hard-disk results. The whole range of densities corresponding to fluids is considered. We show that it is not possible to represent the modes by a two-dimensional vector field of the position perturbations alone (as is the case for hard disks), but the momentum perturbations are simultaneously required for their characterization.


Introduction
For the last 50 years, molecular dynamics simulations have decisively contributed to our understanding of the structure and dynamics of simple fluids and solids [1]. More recently, the concepts of dynamical systems theory have also been applied to study the tangent-space dynamics of such systems [2]. Of particular interest is the extreme sensitivity of the phase-space evolution to small perturbations. On average, such perturbations grow, or shrink, exponentially with time, which may be characterized by a set of rate constants, the Lyapunov exponents. The whole set of exponents is referred to as the Lyapunov spectrum. This instability is at the heart of the ergodic and mixing properties of a fluid and offers a new tool for the study of the microscopic dynamics. In particular, it was recognized very early that there is a close connection with the classical transport properties of systems in nonequilibrium stationary states [3]. For fluids in thermodynamic equilibrium, an analysis of the Lyapunov instability is expected to provide an unbiased expansion of the dynamics into events, which, in favourable cases, may be associated with qualitatively different degrees of freedom, such as the translation and rotation of linear molecules [4], or with the intramolecular rotation around specific chemical bonds [5].
Hard spheres are considered the simplest model for fluids. With respect to the structure, they even serve as a reference system for highly successful perturbation theories of liquids [6,7]. Recently, we studied the Lyapunov instability of such a model in two and three dimensions [8]- [12] and found the following: (i) The slowly-growing and decaying perturbations associated with the non-vanishing Lyapunov exponents closest to zero are related to, and in some cases may even be represented as, periodic vector fields coherently spread out over the physical space and with well-defined wave vectors k. Because of their similarity to the classical modes of fluctuating hydrodynamics, we refer to them as Lyapunov modes. Depending on the boundary conditions, the respective exponents are degenerate, and the spectrum has a step-like appearance in that regime. (ii) The fast-growing or decaying perturbations are localized in space, and the ratio of particles actively contributing at any instant of time converges to zero in the thermodynamic limit.
Experimentally, the Lyapunov modes were found for hard-particle systems in one, two and three dimensions and for various boundary conditions, provided that the system size L in at least one direction is large enough for the discrete particles to generate a recognizable wavelike pattern with a wave number k. Their theoretical understanding is based on the spontaneous breaking of the translational symmetry of the zero modes (k = 0), which are associated with the vanishing Lyapunov exponents and are a consequence of the conservation laws in the system [12]- [17].
Based on this experimental and theoretical evidence for hard-particle systems, it is natural to expect that the Lyapunov modes are a general feature of many-body systems with short-range interactions and that the details of the pair potential should not matter. However, already our very first simulations with soft disks in two dimensions [9,10,18] revealed a much more complicated scenario. Whereas the spatial localization of the fast-growing or decaying perturbations could be easily verified, the mode structure for the slowly growing or decaying perturbations was elusive and could, at first, not be unambiguously detected. Here, we demonstrate that the modes for soft sphere fluids do indeed exist. However, sophisticated Fourier-transformation techniques are needed to prove their existence. Interestingly, the degeneracy of the Lyapunov exponents and, hence, the step structure of the spectrum so familiar from the hard-disk case is recovered, but only for low densities. This suggests that kinetic theory is a proper theoretical framework in that case. For intermediate and large-density soft-disk fluids, the degeneracy disappears. Most recently, Lyapunov modes were demonstrated by Radons and Yang [19,20] for one-dimensional Lennard-Jones fluids at low temperatures and densities.
In this work, we analyse the Lyapunov instability of repulsive Weeks-Chandler-Anderson (WCA) disks in two dimensions and compare it to analogous hard-disk results. The main difficulty is the large box size required for the modes to develop and, hence, the large number of particles, N. Parallel programming is reqired for the computation of the dynamics both in phase and tangent space and for the renormalization of the perturbation vectors according to the classical algorithms of Benettin et al [21] and Shimada et al [22]. In section 2, we characterize the systems and summarize the numerical methods used. In section 3, we discuss the surprising differences found between the Lyapunov spectra for the soft-and hard-particle fluids, particularly at intermediate and large densities. Section 4 is devoted to an analysis of the zero modes associated with the vanishing Lyapunov exponents. The Lyapunov modes are analysed in sections 5 and 6. We close with some concluding remarks in section 7.

Characterization of the system and methodology
We consider N disks with equal mass m in a two-dimensional rectangular box with extensions L x and L y and aspect ratio A ≡ L y /L x . Periodic boundary conditions are used throughout. The particles interact with a smooth repulsive potential, for which we consider two cases: (i) a Weeks-Chandler-Anderson potential, with a force cutoff at r = 2 1/6 σ. As usual, σ and are interaction-range and energy parameters, respectively, which are ultimately set to unity. Actually, such a potential is not the best choice, since the forces and their derivatives (required for the dynamics in tangent-space) are not continuous at the cutoff. This introduces additional noise into the simulation and spoils the conservation laws. It also affects the computation of the Lyapunov spectrum. To avoid this problem, we sometimes also use: (ii) a power-law potential, which looks very similar but does not suffer from this deficiency. We mention, however, that the results discussed below turn out to be insensitive to this weakness of the WCA potential. For reasons of comparison, we also consider: (iii) a hard-disk potential, The phase space X has 4N dimensions, and a phase point ∈ X is given by the 4Ndimensional vector (t) = {q i , p i ; i = 1, . . . , N}, where q i and p i denote the respective positions and momenta of the disks. It evolves according to the time-reversible equations of motioṅ which are conveniently written as a system of first-order differential equations. Here, F( ) = { p i m , − ∂ ∂q i ; i = 1, . . . , N} follows from Hamilton's equations, where ≡ i j>i φ(|q i − q j |) is the total potential energy.
Any infinitesimal perturbation δ l = {δq (l) i , δp (l) i ; i = 1, . . . , N} lies in the tangent space TX, tangent to the manifold X at the phase point (t). Here, δq (l) i and δp (l) i are two-dimensional vectors and denote the position and momentum perturbations contributed by particle i. δ l evolves according to the linearized equations of motioṅ Oseledec's theorem [23] assures us that there exist 4N orthonormal initial tangent vectors δ l (0), l = 1, . . . , 4N, whose norm grows or shrinks exponentially with time such that the Lyapunov exponents exist. However, in the course of time, these vectors would all get exponentially close to the most unstable direction and diverge, since the unconstrained flow in tangent space does not preserve the orthogonality of these vectors. In the classical algorithms of Benettin et al [21] and Shimada et al [22], which are also used in this work, this difficulty is circumvented by replacing these vectors by a set of modified tangent vectors, δ l , l = 1, . . . , 4N, which are periodically re-orthonormalized with a Gram-Schmidt procedure. The Lyapunov exponents are determined from the timeaveraged renormalization factors. The vectors δ l represent the perturbations associated with λ l and are the objects analysed below. Henceforth, we use the notation δq (l) i , and δp (l) i also for the individual particle contributions to the orthonormalized vectors, δ l = {δq (l) i , δp (l) i ; i = 1, . . . , N}. For later use, we introduce also the squared norm in the single-particle µ space, This quantity is bounded, 0 γ (l) i 1, and obeys the sum rule N i=1 γ (l) i = 1 for any l. Since the equations (5) are linear in δ l , the quantities γ (l) i indicate how the activity for perturbation growth, as measured by λ l , is distributed over the particles at any instant of time.
For the hard-disk systems, the equations of motion (4) and (5) need to be generalized to include also the collision maps due to the instantaneous particle collisions. For algorithmic details, we refer to our previous work [8,24,25].
The exponents are taken to be ordered by size, λ l λ l+1 , where the index l numbers the exponents.According to the conjugate-pairing rule for symplectic systems [26,27], the exponents appear in pairs such that the respective pair sums vanish, λ l + λ 4N+1−l = 0. Only the positive-half of the spectrum needs to be calculated. Furthermore, six of the exponents, {λ 2N−2 l 2N+3 }, vanish as a consequence of the constraints imposed by the conservation of energy, momentum, centre of mass (for the tangent-space dynamics only) and of the non-exponential time evolution of the perturbation vector in the direction of the flow.
In the soft-disk case, we use reduced units for which the disk diameter σ, the particles' mass m, the potential parameter and Boltzmann's constant k B are all set to unity. As usual, the temperature is defined according to K = p 2 2m = (N − 1)k B T , where K is the total kinetic energy and . . . is a time average. For hard disks, K is a constant of the motion and K/N is taken as the unit of energy. Since Lyapunov exponents basically scale with the particle velocitiesstrictly for hard disks and approximately for soft disks, all comparisons between soft and hard particles must be taken for equal temperatures. We use T = 1 throughout.
One might think that, instead of σ = 1, an effective (temperature-dependent) hard-disk diameter σ eff (T ) should be used for a comparison of Lyapunov exponents of WCA and harddisk particles. With the method of Barker and Henderson [28], one finds σ eff (1) = 1.016, which increases the hard-disk exponents by the same insignificant factor. We have therefore used σ = 1 throughout.
The computation of full Lyapunov spectra for soft particles is a numerical challenge, since the number of simultaneously integrated differential equations increases with the square of the particle number. Therefore, parallel processing with up to six nodes is used, both for the integration and for the Gram-Schmidt re-orthonormalization procedure [29]. A fourth-order Runge-Kutta algorithm is used for the integration of the equations of motion.

Phenomenology
We start our comparison of the spectra for soft (WCA) and hard-disks with low-density gases, ρ ≡ N/L x L y = 0.2, at a common temperature T = 1. The system consists of N = 80 particles in a rather elongated rectangular simulation box, L x = 100, L y = 4; A = L y /L x = 0.04, with periodic boundary conditions. An inspection of figure 1 (left panel) shows, as expected, that the two systems have similar spectra and, thus, similar chaotic properties: the maximum exponents λ 1 and the shapes of the spectra are reasonably close. For even lower densities (not shown) the agreement is even better. Most important, however, is the observation that the step structure due to the degeneracy of the small Lyapunov exponents is observed both for the hard disks, for which it is well known to exist [12,24,30,31], and for the soft-disk gases at low densities. We conclude from this comparison in figure 1 that the Lyapunov modes exist for low-density soft disks. In view of the smallness of L y and the quasi-one-dimensional nature of the system, these modes may have only wave vectors parallel to the x-axis of the box. This general picture also persists for intermediate gas densities, ρ = 0.4, as shown in the right panel of figure 1. The box size is the same as before, but the number of particles (and, as a consequence, also the number of exponents) is doubled. The spectral shapes for the WCA and hard-disk fluids, and the maximum exponents λ 1 in particular, become significantly different. Chaos, measured by λ 1 , is obviously weakened for the WCA particles with a finite collision time as compared to the instantaneously colliding hard disks. Interestingly, the Kolmogorov-Sinai (or dynamical) entropy h KS , which is equal to the sum of all positive exponents [32,33], differs less due to the compensating effect of the intermediate exponents in the range 0.2 < l/2N < 0.8. This may be verified in figure 4 below.
The systems of figure 1 are quasi-one-dimensional and do not represent bulk fluids. We show in figure 2 analogous spectra for bulk systems in a rectangular simulation box with aspect ratio A = 0.6 and containing N = 375 particles. The density varies between 0.1 and 0.4 and is specified by the labels. For all spectra, the lowest step is clearly discernible, but   is progressively less pronounced if the density is increased. For ρ = 0.1, the four smallest exponents (744 l 747) are identified as belonging to longitudinal modes (see below) and the two next-larger exponents (742 l 743) to transverse modes. In figure 3, we compare a WCA system to a hard-disk system in a square simulation box at a density ρ = 0.4 and N = 400 particles. The length of the simulation box in every direction is L x = L y = 31.62. The step structure that is so prominent in the hard-disk spectrum vanishes altogether for the soft potential case. If the density of such a square system with N = 400 particles is decreased by increasing the size of the simulation box, the step structure with the same degeneracy reappears.

Measures for chaos
The maximum Lyapunov exponent λ 1 and the Kolmogorov-Sinai entropy h KS are generally accepted as measures for chaos. The latter corresponds to the rate of information on the initial conditions generated by the dynamics and, according to Pesin's theorem [32,33], is equal to the sum of the positive Lyapunov exponents, h KS = λ l >0 λ l . In figure 4, we compare the isothermal density dependences of these quantities for 375-disk WCA systems with corresponding hard-disk results. Also included in this set of figures is a comparison for the smallest positive exponent, λ 2N−3 , which is directly affected by the mechanism generating the Lyapunov mode with the longest possible wave length (if it exists). For very dilute hard-disk gases, λ 1 and h KS /N have been estimated from kinetic theory by Dorfman et al [34], who have provided expansions to leading orders in ρ, λ 1 = Aρ(− ln ρ + B) + · · · and h KS /N =Āρ(− ln ρ +B) + · · ·, with explicit expressions for the constants A, B and A,B. These predictions were very successfully confirmed by computer simulations [35,36]. Computer simulations by Dellago et al also provided hard-disk results over the whole range of densities [24]. λ 1 and h KS /N were found to increase monotonically with ρ with the exception of a loop in the density range between the freezing point of the fluid (ρ f = 0.882 [37,38]) and the melting point of the solid (ρ s = 0.912 [37,38]). This loop disappears if λ 1 and h KS /N are plotted as a function of the collision frequency ν instead of the density [24]. If ρ approaches the close-packed density ρ 0 = 1.1547/σ 2 , both λ 1 and h KS /N diverge as a consequence of the divergence of ν.
The results for the WCA disks reported here are for 0.1 ρ 1.1, which covers the whole fluid range and even extends to the solid. The densities for fluid-solid coexistence are similar to those of hard disks but have not been determined with precision. Further simulations in that regime are currently under way [39]. Not unexpectedly, we infer from figure 4 (top panel) that λ 1 -WCA approaches λ 1 -HD for small densities ρ < 0.1. However, the density dependence for larger ρ is qualitatively different. λ 1 approaches a maximum near the density corresponding to the fluid-to-solid phase transition, which confirms our previous results for Lennard-Jones fluids and solids [40]. The cooperative dynamics at such a transition causes maximum chaos in phase space. In the solid regime, λ 1 decreases with increase in density, which is easily understandable in view of the fact that harmonic solids are not chaotic at all.
The smallest positive exponent λ 2N−3 for soft-disk gases follows the hard-disk density behaviour more closely than λ 1 for lower densities, as inferred from figure 4 (middle panel). This means that the associated collective dynamics is affected less by the details of the pair potential. Significant deviations start to occur for ρ > 0.3. For h KS /N, the agreement between soft and hard disks extends even to a larger range of densities, as may be seen in figure 4 (bottom panel). However, this observation is deceptive and does not necessarily signal a similar dynamics. It is a consequence of significant cancellation due to the exponents intermediate between the maximum and the small exponents of the spectra as in figure 1 (right panel).

Localization of tangent-space perturbations in physical space
One may interpret the maximum (minimum) Lyapunov exponent as the rate constant for the fastest growth (decay) of a phase-space perturbation. Thus, it is dominated by the fastest dynamical events, binary collisions in the case of particle fluids. Therefore, it does not come as a surprise that the associated tangent vector has components which are strongly localized in physical space [4]. Similar observations for other spatially extended systems have been made previously by various authors [41]- [44]. With the help of a particular measure for the localization [11], we could show that, for both hard-and soft-disk systems, the localization persists even in the thermodynamic limit, such that the fraction of tangent-vector components contributing to the generation of λ 1 at any instant of time converges to zero with N → ∞ [9,31]. The localization becomes gradually worse for larger Lyapunov indices l > 1, until it ceases to exist and (almost) all particles collectively contribute to the perturbations associated with the smallest Lyapunov exponents, for which coherent modes are known (or believed) to exist [9,31].
Here, we adopt another entropy-based localization measure due to Taniguchi and Morriss [45], which was successfully applied to quasi-one-dimensional hard-disk gases. Recalling the definition of the individual particle contributions γ (l) i to δ 2 l (≡ 1) in equation (7), one may introduce an entropy-like quantity where . . . denotes a time average. The number may be taken as a measure for the spatial localization of the perturbation δ l associated with λ l , 1 l 2N. It is bounded according to 1 W (l) N. The lower bound 1 indicates the most localized state with only one particle contributing, and the upper bound N signals uniform contributions by all the particles {γ (l) i = 1/N, i = 1, . . . , N}. We shall refer to the set of normalized localization parameters {W (l) /N, l = 1, . . . , 2N} as the localization spectrum. It should be noted that the Euclidean norm is used in the definition of γ (l) i and that W (l) /N depends on this choice. Qualitatively, it suffices to demonstrate localization.
In figure 5, we compare the localization spectra for the same hard-and soft-disk systems for which the Lyapunov spectra are given in figure 1 (L x = 100, L y = 4, A = 0.04). For the lower density gases in the left panel (ρ = 0.2, N = 80), the localization spectra for soft and hard disks are almost indistinguishable. They indicate strong localization for the perturbations belonging to the large exponents (small l), and collective behaviour for large l. Of particular interest is the comb-like structure also magnified in the inset, which is a consequence of the Lyapunov modes (which will be discussed in more detail below). Stationary transverse modes lead to large values of W (l) /N, indicating strong collectivity. The propagating LP modes [12], however, are characterized by a significantly reduced value of W (l) /N. For intermediate gas densities in the right panel (ρ = 0.4, N = 160), differences between the hard-and soft-disk results become apparent. Most prominently, the comb structure in the localization spectrum of the WCA system has almost disappeared, although the Lyapunov spectrum clearly displays steps in the small-exponent regime, which is a clear indication of the existence of modes.

Zero modes
The dynamics in phase space and tangent space is strongly affected by the inherent symmetries of a system, i.e. infinitesimal transformations leaving the equations of motion invariant. They are intimately connected with the conservation laws obeyed by the dynamics. If these transformations act as infinitesimal perturbations of the initial conditions, the latter do not grow/shrink exponentially (but at the most linearly) with time and give rise to as many vanishing Lyapunov exponents as there are symmetry-related constraints. The generators of these transformations are taken as unit vectors in tangent space, which point to the direction of the respective perturbation and are called zero modes. They span an invariant subspace N ( ) of the tangent space at any phase point , which is referred to as the zero subspace. They are essential for an understanding of the Lyapunov modes dealt with in the following. For a planar system of hard or soft disks with periodic boundaries, there are six vanishing Lyapunov exponents. There must be six symmetry-related initial perturbations leading to nonexponential growth [12,16]: (1) Translation invariance in x and y directions. This gives rise to the generators e 1 and e 2 explicitly given below. It is a consequence of momentum conservation. (2) Invariance with respect to Galilei transformations. This is consistent with generators e 3 and e 4 below, and may be viewed as a consequence of the conservation of the initial phase. (3) Energy conservation. This has been associated with a generator e 5 orthogonal to the energy hyper-surface. This turns out to be correct for noninteracting particles such as streaming hard disks (with vanishing forces F x,i = F y,i = F z,i = 0 below). However, for soft disks, there is a difficulty, as will be shown below.
(4) Homogeneous time shift with generator e 6 , corresponding to no perturbation growth in the direction of the phase flow.
If we use the explicit notation for an arbitrary tangent vector, one obtains where α is a normalizing factor. Each vector has 4N dimensions. To check whether the orthogonal vectors e 1 to e 6 form a basis for the invariant zero space, we consider the expansion of the six orthonormal tangent vectors δ 2N−2 , . . . , δ 2N+3 , responsible for the six vanishing exponents in the simulation, in that basis, where α l,j ≡ (δ l · e j ). Since all vectors involved are of unit length, α l,j may either be interpreted as the projection of δ 2N−3+l onto e j or vice versa. For hard disk fluids, one can easily show [12] that the tangent vectors δ 2N−2 , . . . , δ 2N+3 are fully contained in N ( ). The projection cosines strictly obey the sum rules  (e 1 , . . . , e 6 ). Note that there are no forces acting on the particles in this case, and the particles are moving on straight lines between successive instantaneous collisions.
For soft-disk systems, however, the situation is more complicated. In table 1, we list instantaneous matrix elements α l,j for a simple example, N = 4 particles interacting with the particularly smooth, repulsive potential of equation (2). The density, ρ = 0.1, is low enough for the isolated binary collisions to be easily identified. Two instances are considered: one just Table 1. Projection cosines, α l,j , according to equation (16) for N = 4 soft particles in a periodic box and interacting with the power-law potential (2 at the beginning of a collision of particles 1 and 2 (table 1(A)), and a time half way through this collision (table 1(B)). Considering the columns first, the squared sums 6 l=1 (α l,j ) 2 always add to unity for j = 1, . . . , 4, which means that e 1 to e 4 are fully contained in the subspace Span(δ 2N−2 , . . . , δ 2N+3 ) of the tangent space. The same is not true, however, for e 5 and e 6 , which contain in their definition the instantaneous forces acting on the colliding particles. If no collisions take place anywhere in the system (table 1(A)), the subspaces Span(δ 2N−2 , . . . , δ 2N+3 ) and Span(e 1 , . . . , e 6 ) are nearly the same, but not quite, as the squared projection-cosine sums indicate. The moment a collision occurs, the sums 6 l=1 (α l,j ) 2 for j = 5 and 6 may almost vanish, as it happens in table 1(B), and the vectors e 5 and e 6 become nearly orthogonal to Span(δ 2N−2 , . . . , δ 2N+3 ). In figure 6, an attempt is made to illustrate this relationship for such a high-dimensional tangent space. For systems containing many particles and/or for larger densities, there will always be at least one collision in progress, and the horizontal sums, 6 j=1 (α l,j ) 2 for l = 1, . . . , 6, and the vertical sums, 6 l=1 (α l,j ) 2 for j = 5 and 6, fluctuate and assume values significantly less than unity. The subspaces Span(δ 2N−2 , . . . , δ 2N+3 ) = D + ⊕D − (see the definition in figure 6) and N = Span(e 1 , . . . , e 6 ) do not agree in general. The culprit is e 5 , which is not a correct zero vector for soft-disk systems. We know from the symplectic nature of the equations of motion that a zero mode replacing e 5 must exist. So far, we have not been able to construct it.

Lyapunov modes
Lyapunov modes are periodic spatial patterns observed for the perturbations associated with the small positive and negative Lyapunov exponents. They are characterized by wave vectors k where a rectangular box with periodic boundaries is assumed, and n x and n y denote the number of nodes parallel to x and y, respectively. For modifications due to other boundary conditions, we refer to [12,46]. Experimentally, Lyapunov modes were observed for hard particle systems in one, two and three dimensions [8,10,12,17,45], for hard planar dumbbells [4,11,47] and, most recently, for one-dimensional soft particles [19,20]. Theoretically, they are interpreted as periodic modulations (k = 0) of the zero modes, to which they converge for k → 0. Since a modulation, for k > 0, involves the breaking of a continuous symmetry (translational symmetry of the zero modes), they have been identified as Goldstone modes [16], analogous to the familiar hydrodynamic modes and phonons. For the computation of the associated wavevector-dependent Lyapunov exponents governing the exponential time evolution of the Lyapunov modes, random matrices [13,30], periodic-orbit expansion [48] and, most successfully, kinetic theory have been employed [14]- [16]. For hard-disk systems, the representation of the modes is simplified by the fact that, for large N, the position perturbations, {δq (l) i , i = 1, . . . , N}, and momentum perturbations, . . , N}, may be viewed as two-dimensional vector fields ϕ (l) q (x, y) and ϕ (l) p (x, y), respectively, which turn out to be nearly parallel (for λ > 0) or anti-parallel (for λ < 0) [12]. To illustrate this point, we plot in figure 7 cos where l is the angle between the 2N-dimensional vectors comprising the position and momentum perturbations of δ l , when they are viewed in the same 2N-dimensional space.
where C (l) i is a positive number, which is almost the same for all particles i. Once the δq (l) i are known, the δp (l) i may be obtained from this relation. For a characterization of the modes, it suffices to consider only the position perturbations, which are interpreted as a vector field, ϕ (l) q (x, y), over the simulation cell [12]. To this level of approximation, the momentum perturbations are omitted.
For soft-spheres, the situation is more complicated. For low-density systems, cos l behaves as for hard disks, which is expected in view of the similarity of the Lyapunov spectra. For dense fluids, however, the behaviour of this quantity is qualitatively different: it does not converge to unity but to zero in the interesting mode regime. This is demonstrated in figure 7 by the lower red curve, which is for N = 375 soft WCA disks at the same density ρ = 0.4. The larger the collision frequency, the larger will be the deviations of cos l from 1 (for λ l > 0) or −1 (for λ l < 0). This result indicates that an approximate representation of a Lyapunov mode in terms of a single (periodic) vector field of the position (or momentum) perturbations of the particles is possible only for low-density systems. For larger densities and/or larger systems with many collisional events taking place at the same time, one needs to simultaneously consider all the perturbation components δx i , δy i , δp x,i , δp y,i of a particle. To understand this unexpected behaviour in more detail, we display in figure 8 (left panel) the time dependences of cos l belonging to the maximum (l = 1) and to the smallest positive (l = 2N − 3 = 5) exponents of a four-particle system in a square periodic box with a density ρ = 0.1. During a streaming phase without collisions, the 2N-dimensional vectors δq (l) and δp (l) tend to become parallel as required by the linearized free-flight equations in tangent space. However, this process is disrupted by a collision, which reduces (and in some cases even reverses the sign of) cos l . In the following streaming phase, i.e. forward in time, l (t) relaxes towards zero. This suggests that the numerical time evolution of cos l and, hence, of δ l is not invariant with respect to time reversal. This is indeed the case, as demonstrated figure 8 (right panel). To construct this figure, the phase space trajectory was stored for another 10 000 time units and, after a time-reversal transformation {q i → q i , p i → −p i ; i = 1, . . . , N}, it was consecutively used (in reversed order) as the reference trajectory for the reversed tangent-space dynamics. In figure 8 (right panel), we show cos 16 and cos 12 , which belong to the minimum and largestnegative exponents (of the time-reversed dynamics), respectively, and which should be compared to cos 1 and cos 5 for the forward evolution. Although the same collisions are involved, the curves look totally different. This confirms previous results concerning the lack of symmetry for the forward and backward time evolution of the Gram-Schmidt orthonormalized tangent vectors δ l for Lyapunov-unstable systems [49,50]. Furthermore, we have numerically verified that cos( l (t)) = cos( 4N−l+1 (t)) is always obeyed, both forward and backward in time.

Fourier-transform analysis
A simple inspection of the vector field ϕ (l) q (x, y) very often is not sufficient to unambiguously determine the wave vector of a mode, particularly for smaller wavelengths than for maximum possible wavelengths. Therefore, Fourier transformation methods have been used [51], where we regard ϕ a (x, y) as a spatial distribution The a i are identified with the particle perturbations δx (l) i , δy (l) i , δp (l) x,i , δp (l) y,i . In view of the periodicity of the box, the Fourier coefficients have wavenumbers given by equation (17). They are computed from The power spectrum is defined by We have also applied the algorithm for unequally spaced points by Lomb [52,53], suitably generalized to two-dimensional transforms. As a first example, the smooth lines in figure 9 show the power spectra P(k) for the transverse modes indexed by l = 317, 311, 305, 299, 293 and 287 of the WCA system described in figure 1 (right panel). They correspond to the successive wavenumbers k 1,0 to k 6,0 as predicted by equation (17). The dashed lines depict the power spectra for the longitudinal modes l = 315, 309 with wavenumbers k 1,0 and k 2,0 . In view of the discussion in the following section, however, the characterization of the modes as transversal or longitudinal is difficult and not unambiguous. In a second example, we show in figure 10 the Lyapunov spectra (left panel) of WCA particles in a fixed elongated box, L x = 100, L y = 4, for various densities varying from ρ = 0.1 to ρ = 0.7 as indicated by the labels. In figure 10 (right panel) the power spectra P(k) for the modes corresponding to the smallest positive exponent, λ l=2N−3 , of each spectrum are shown. Since L x is the same in all cases, all power spectra have a peak at the allowed wave number k 1,0 = 0.063. This peak is also well resolved for large densities, for which the step structure in the spectrum is blurred.

Concluding remarks
In the foregoing sections, we have demonstrated the existence of Lyapunov modes in soft-disk fluids with the help of various indicators such as the localization measure in section 3.3, and the Fourier analysis in the previous section. However, the classification and characterization of the modes is more complicated than for hard disks. This is demonstrated in figure 11, where the position perturbations δx 746 i (bottom) and δy 746 i (top) are plotted at the particle positions in the simulation cell for the 375-particle system with a density ρ = 0.4 familiar from figure 2. Näively, the upper surface would have to be classified as a transverse mode with the perturbation component δy perpendicular to a wave vector k parallel to x and the lower surface as a longitudinal mode with δx parallel to k. Since the momentum perturbations are not nearly as parallel to the position perturbations (see figure 7) as is (approximately) the case for hard disks, it is not sufficient to represent a mode by a single two-dimensional vector field such as ϕ (l) q (x, y) or ϕ (l) p (x, y). Position and momentum perturbations need to be considered simultaneously for a characterization of the modes. Furthermore, transversal and longitudinal modes do not remain uncoupled. The disappearance of the step structure in the Lyapunov spectrum and, hence, of the degeneracy for larger densities points to a complicated tangent-space dynamics we have not yet been able to unravel. We do not know of any theory which accounts for all the numerical results presented in this paper.