What does dynamical systems theory teach us about fluids?

1.1 MB·English

by  Hadrien Bosetti

What does dynamical systems theory teach us about fluids?
Hadrien Bosetti1,2 and Harald A. Posch1
1 Computational Physics Group, Faculty of Physics,
University of Vienna, Boltzmanngasse 5, 1090 Wien, Austria and
2 Complex Energy Systems Group, Austrian Institute of Technology, Giefinggasse 6, 1210 Wien, Austria
(Dated: January 19, 2015)
We use molecular dynamics simulations to compute the Lyapunov spectra of many-particle sys-
temsresemblingsimplefluidsinthermalequilibriumandinnon-equilibriumstationarystates. Here
we review some of the most interesting results and point to open questions.
5
1 PACSnumbers: 05.10.-a,05.40.-a,05.45.-a,51.10.+y
0 Keywords: dynamicalsystemtheory;Lyapunovinstability;hydrodynamicmodes;phasetransitions;station-
2 arynonequilibriumstates

n
a I. INTRODUCTION Thestructureofasimplefluidisessentiallydetermined
J
by the steep (e.g. ∝ r−12) repulsive part of the pair
6
Whatdoesdynamicalsystemstheorycontributetoour potential, which can be well approximated by a discon-
1
understanding of matter? To a large extent, the royal tinuous hard core. In perturbation theories, this hard
] road to gain an understanding of fluids or solids has potential may be taken as a reference potential with the
h been statistical mechanics. Based on interaction poten- long-rangeattractivepotential(∝−r−6)actingasaper-
c tials obtained from experiments and quantum mechani- turbation [4]. Hard disk fluids in two dimensions, and
e
m cal simulations, sophisticated perturbation theories are hard sphere systems in three, are easy to simulate and
capable of providing a quantitative description of the are paradigms for simple fluids. The first molecular dy-
-
t static structural properties of such fluids on the atom- namics simulations for hard sphere systems were carried
a isticscale. Computersimulationshaveprovidedguidance outbyAlderandWainwright[5]in1957,thefirstcompu-
t
s andinvaluableinsighttounraveltheintricatelocalstruc- tation of Lyapunov spectra for such systems by Dellago,
.
t ture and even the non-equilibrium dynamics of “simple” Posch and Hoover [6] in 1996.
a
m liquids including hydrodynamic flows and shock waves. There exist numerous extensions of the hard-sphere
At present, this program is being extended to ever more model to include rotational degrees of freedom and var-
d- complex molecular fluids and/or to systems confined to ious molecular geometries [7]. Arguably the simplest,
n particular geometries such as interfaces and pores, and which still preserves spherical symmetry, are spheres
o to fluids out of thermal equilibrium. with a rough surface, so-called ”rough hard spheres”
c We have raised a related question, “What is liquid?” [8]. In another approach, fused dumbbell diatomics are
[ in1998inasimilarcontext[1]. Afluiddiffersfromasolid used to model linear molecules [9, 10]. Both models
1 by the mobility of its particles, and this ability to flow are used to study the energy exchange between trans-
v is a collective phenomenon. The flow spreads rapidly in lational and rotational degrees of freedom and, hence,
9 phase space, which constitutes the fundamental instabil- rotation-translation coupling for molecules with moder-
0 ity characteristic of a gas or liquid. About thirty years ateanisotropy. Otherapproachesmoresuitableforlarger
9
ago, dynamical systems theory provided new tools for molecular anisotropies involve spherocylinders [11, 12],
3
0 the characterization of the chaotic evolution in a multi- ellipsoids [13, 14] and even inflexible hard lines [15]. To
. dimensional phase space, which were readily applied to ourknowledge,onlyforthefirsttwoschemes,namelyfor
1 liquids shortly after [2, 3]. The main idea is to follow rough hard disks [16, 17] and for planar hard dumbbells
0
5 not only the evolution of a state in phase space but, si- [1,18,19],extensivestudiesoftheLyapunovspectrahave
1 multaneously,theevolutionofvarioustinyperturbations beencarriedout. Weshallcomebacktothisworkbelow.
: applied to that state at an initial time and to measure Therearenumerouspapersforrepulsivesoft-potential
v
the growth or decay of such perturbations. The study systems and for Lennard-Jones systems from various au-
i
X of the Lyapunov stability, or instability, of a trajectory thors, in which an analysis of the Lyapunov instability
r withrespecttosmallperturbationsalwayspresentinreal has been carried out. We shall make reference to some
a physicalsystemsishopedtoprovidenewandalternative of this work in the following sections.
insight into the theoretical foundation and limitation of The paper is organized as follows. In the next section
dynamical models for dense gases and liquids, of phase weintroduceglobalandlocal(timedependent)Lyapunov
transitions involving rare events, and of the working of exponents and review some of their properties. Of par-
the Second Law of thermodynamics for stationary non- ticular interest are the symmetry properties of the local
equilibrium systems. It is the purpose of this review to exponentsandofthecorrespondingperturbationvectors
assess how far this hope has become true and how much in tangent space. Both the familiar orthonormal Gram-
our understanding of fluids has gained from the study of Schmidt vectors and the covariant Lyapunov vectors are
the Lyapunov instability in phase space. considered. In Sec. III we study planar hard disk sys-
2
tems over a wide range of densities, and compare them where Dφt defines the flow in tangent space and is rep-
to analogous fluids interacting with a smooth repulsive resented by a real but generally non-symmetric D ×D
potential. There we also demonstrate that the largest matrix . The dynamical (or Jacobian) matrix,
(absolute) Lyapunov exponents are generated by pertur-
bations, which are strongly localized in physical space: ∂F
J(Γ)≡ , (4)
only a small cluster of particles contributes at any in- ∂Γ
stant of time. This localization persists in the thermo-
determines, whether the perturbation vector δΓ(t) has
dynamic limit. Sec. IV is devoted to another property of
thetendencytogroworshrinkataparticularpointΓ(t)
the perturbations associated with the small exponents,
in phase space the system happens to be at time t. Ac-
the so-called Lyapunov modes. In a certain sense, the
cordingly, the matrix element
Lyapunov modes are analogous to the Goldstone modes
offluctuatinghydrodynamics(suchasthefamiliarsound
δΓ†(Γ)J(Γ)δΓ(Γ) (5)
and heat modes). However, it is surprisingly difficult to
establish a connection between these two different view- turnsouttobepositiveornegative,respectively. Here,†
points. InSec.Vparticlesystemswithtranslationaland
means transposition. If, in addition, the perturbation is
rotational dynamics are considered. Two simple planar
normalized, ||δΓ|| = 1, and points into particular direc-
model fluids are compared, namely gases of rough hard
tions in tangent space to be specified below, this matrix
disksandofhard-dumbbellmolecules. Closesimilarities,
elementturnsouttobealocalrateforthegrowthorde-
but also surprising differences, are found. Most surpris-
cayof||δΓ(t)||andwillbereferredtoasalocalLyapunov
ingisthefastdisappearanceoftheLyapunovmodesand
exponent Λ(Γ) at the phase point Γ.
the breakdown of the symplectic symmetry of the local
Gram-Schmidt exponents for the rough hard disk sys-
tems, if the moment of inertia is increased from zero. In
A. Covariant Lyapunov vectors
Sec. VI we summarize some of the results for station-
ary systems far from thermal equilibrium, for which the
In 1968 Oseledec published his celebrated multiplica-
Lyapunov spectra have been valuable guides for our un-
tive ergodic theorem [20–23], in which he proved that
derstanding. For dynamically thermostatted systems in
under very general assumptions (which apply to molecu-
stationary non-equilibrium states, they provide a direct
lar fluids) the tangent space decomposes into subspaces
link with the Second Law of thermodynamics due to the
E(j) with dimension m ,
presence of a multifractal phase-space distribution. We j
conclude with a few short remarks in Sec. VII.
TM(Γ)=E(1)(Γ)⊕E(2)(Γ)⊕···⊕E(L)(Γ), (6)
foralmostallΓ∈M (withrespecttotheLesbeguemea-
II. LYAPUNOV EXPONENTS AND sure), such that (cid:80)L m = D. These subspaces evolve
PERTURBATION VECTORS j=1 j
according to
Let Γ(t) = {p,q} denote the instantaneous state of Dφt| E(j)(Γ(0))=E(j)(Γ(t)), (7)
Γ(0)
a dynamical particle system with a phase space M of
dimension D. Here, p and q stand for the momenta and are said to be covariant, which means that they co-
and positions of all the particles. The motion equations move – and in particular co-rotate – with the flow in
are usually written as a system of first-order differential tangent space. In general they are not orthogonal to
equations, each other. If v(Γ(0)) ∈ E(j)(Γ(0)) is a vector in the
Γ˙ =F(Γ), (1) subspace E(j)(Γ(0)), it evolves according to
whereFisa(generallynonlinear)vector-valuedfunction Dφt|Γ(0) v(Γ(0))=v(Γ(t))∈E(j)(Γ(t)). (8)
of dimension D. The formal solution of this equation is
written as Γ(t) = φt(Γ(0)), where the map φt : Γ → Γ The numbers
defines the flow in M, which maps a phase point Γ(0) at (cid:13) (cid:13)
1 (cid:13)v(Γ(t))(cid:13)
time zero to a point Γ(t) at time t. (±)λ(j) = lim ln(cid:13) (cid:13) (9)
Next, we consider an arbitrary infinitesimal perturba- t→±∞|t| (cid:13)v(Γ(0))(cid:13)
tion vector δΓ(t), which depends on the position Γ(t) in
for j ∈ {1,2,...,L} exist and are referred to as the
phase space and, hence, implicitly on time. It evolves
(global)Lyapunovexponents. Here,theupperindex(+)
according to the linearized equations of motion,
or(−)indicates,whetherthetrajectoryisbeingfollowed
δ˙Γ=J(Γ)δΓ. (2) forward or backward in time, respectively. If a subspace
dimension m is larger than one, then the respective ex-
j
This equation may be formally solved according to ponent is called degenerate with multiplicity m . If all
j
exponents are repeated according to their multiplicity,
δΓ(t)=Dφt| δΓ(0), (3)
Γ(0)
3
there are D exponents altogether, which are commonly B. Orthogonal Lyapunov vectors
ordered according to size,
AnotherdefinitionoftheLyapunovexponents,alsopi-
(+)λ ≥···≥(+) λ , (10)
1 D oneered by Oseledec [20–23], is via the real and symmet-
(−)λ1 ≥···≥(−) λD, (11) rical matrices
(cid:104) (cid:105) 1
lim Dφt|† Dφt| 2|t| , (16)
where the subscripts are referred to as Lyapunov index. t→±∞ Γ(0) Γ(0)
The vetors v(cid:96) generating λ according to
(cid:96)
whichexistwithprobabilityone(bothforwardandback-
1 (cid:13)(cid:13)v(cid:96)(Γ(t))(cid:13)(cid:13) wardintime). Their(real)eigenvaluesinvolvetheglobal
(±)λ(cid:96) =t→li±m∞|t| ln(cid:13)(cid:13)v(cid:96)(Γ(0))(cid:13)(cid:13); (cid:96)=1,2,...,D, Lyapunov exponents,
(12) exp((±)λ )≥···≥exp((±)λ ), (17)
1 D
are called covariant Lyapunov vectors. This notation,
whichtreatsthetangentspacedynamicsintermsofaset where,asbefore,degenerateeigenvaluesarerepeatedac-
of vectors, is more convenient for algorithmic purposes, cording to their multiplicities. For non-degenerate and
although the basic theoretical objects are the covariant degenerate eigenvalues the > and = signs apply, re-
subspaces E(j); j =1,2,...,L. spectively. The corresponding eigenspaces, U(j); j =
±
Degenerate Lyapuov exponents appear, if there exist 1,...,L, are pairwise orthogonal and provide two addi-
intrinsic continuous symmetries (such as invariance of tional decompositions of the tangent space at almost ev-
the Lagrangian with respect to time and/or space trans- ery point in phase space, one forward (+) and one back-
lation,givingrisetoenergyand/ormomentumconserva- ward (−) in time:
tion,respectively). Forparticlesystemssuchsymmetries
almostalwaysexist. Someconsequenceswillbediscussed TM(Γ)=U(1)(Γ)⊕U(2)(Γ)⊕···⊕U(L)(Γ). (18)
± ± ±
below.
The global exponents for systems evolving forward or In each case, the dimensions m of the Oseledec sub-
j
backward in time are related according to spacesU(j); j =1,...,L,haveasumequaltothephase-
±
(+)λ =−(−)λ ; (cid:96)∈{1,2,...,D}. (13) space dimension D: (cid:80)Lj=1mj =D. These subspaces are
(cid:96) D+1−(cid:96) not covariant.
The subspaces E(j) and, hence, the covariant Lyapunov Theclassicalalgorithmforthecomputationof(global)
vectors v(cid:96) generally are not pairwise orthogonal. Lyapunov exponents [24–26] carefully keeps track of all
If in Eq. (9) the time evolution is not followed from d-dimensionalinfinitesimalvolumeelementsδV(d) which
(almost always) evolve according to
zero to infinity but only over a finite time interval τ >0,
so-calledcovariantfinite-time dependentLyapunovexpo-
(cid:32) d (cid:33)
nents are obtained, (cid:88)
δV(d)(t)≈δV(d)(0)exp λ t . (19)
(cid:96)
(±)Λ¯τ(cid:96),cov = τ1 ln (cid:13)(cid:13)Dφ±τ|Γ(0) v(cid:96)(Γ(0)) (cid:13)(cid:13). (14) This is algorithmically achieved w(cid:96)=it1h the help of a
set of perturbation vectors (±)g (Γ); (cid:96) = 1,...,D,
In the limit τ → 0 so-called local Lyapunov exponents (cid:96)
which are either periodically re-orthonormalized with
are generated,
a Gram-Schmidt (GS) procedure (equivalent to a QR-
(±)Λc(cid:96)ov(Γ)=t→lim±0|1t| ln (cid:13)(cid:13) Dφt|Γ v(cid:96)(Γ) (cid:13)(cid:13) dsteacyomorptohsointioornm)a[l27w],itohraarteriaconngtuilnauromusalytricxonosftrLaaignreadngtoe
(cid:104) (cid:105)† multipliers [2, 28–31]. For historical reasons we refer to
= (±)v(cid:96)(Γ) J(Γ) (±)v(cid:96)(Γ) (15) them as GS-vectors. They have been shown to be the
spanning vectors for the Oseledec subspaces U(j) [32].
where ||v(cid:96)(Γ)||=1 is required. The second equality has Accordingly, they are orthonormal but not covar±iant.
a structure as in Eq. (5), applied to the covariant vec- AnalogoustoEqs.(14)and(15),onemaydefinefinite-
tors, and indicates that the local exponents are point time-dependent GS Lyapunov exponents (±)Λ¯τ,GS and
functions, which only depend on Γ. The finite-time- (cid:96)
dependent exponents may be viewed as time averages local GS exponents (±)ΛG(cid:96) S(Γ), where the latter are
of local exponents over a stretch of trajectory during the again point functions,
finite time τ, the global exponents as time averages over
aninfinitely-longtrajectory. Forthelattertobecomedy- (±)ΛG(cid:96) S(Γ)=[(±)g(cid:96)(Γ)]† J(Γ) (±)g(cid:96)(Γ). (20)
namicalinvariants,onerequiresergodicity,whichwewill
As before, the finite-time and global GS exponents are
alwaysassumeinthefollowingforlackofotherevidence.
time averages of the respective local exponents over a
finite respective infinite trajectory [31, 32].
4
C. Symmetry properties of global and local
exponents
For ergodic systems the global exponents are averages
of the local exponents over the natural measure in phase
space and, according to the multiplicative ergodic theo-
remforchaoticsystems, theydonotdependonthemet-
ricandthenormoneappliestothetangentvectors. Also
the choice of the coordinate system (Cartesian or polar,
for example) does not matter. For practical reasons the
Euclidian (or L2) norm will be used throughout in the
following. It also does not matter, whether covariant
or Gram-Schmidt vectors are used for the computation.
The global exponents are truly dynamical invariants.
This is not the case for the local exponents. They de-
pend on the norm and on the coordinate system. And
they generally depend on the set of basis vectors, covari-
ant or GS, as was mentioned before. Furthermore, some
FIG. 1. (Color online) Covariant and GS vectors for a point symmetry properties of the equations of motion are re-
O on the H´enon attractor (light-blue). The black line is a flected quite differently by the two local representations.
finite-time approximation of the (global) stable manifold.
• Local Gram-Schmidt exponents: During the con-
struction of the GS vectors, the changes of the dif-
From all GS vectors only (±)g (Γ), which is associ- ferential volume elements δV(d) following Eq. (19)
1
are used to compute the local exponents. If the
ated with the maximum global exponent, evolves freely
total phase volume is conserved as is the case for
without constraints which might affect its orientation in
tangent space. Therefore it agrees with (±)v (Γ) and is time-independentHamiltoniansystems,thefollow-
1
ing sum rule holds for almost all Γ:
also covariant. Also the corresponding local exponents
agree for (cid:96) = 1, but generally not for other (cid:96). However,
D
thelocalcovariantexponentsmaybecomputedfromthe (cid:88)ΛGS(Γ)=0. (21)
(cid:96)
localGSexponents,andviceversa,iftheanglesbetween
(cid:96)=1
them are known [31].
In this symplectic case we can even say more. For
As an illustration of the relation between covariant
each positive local GS exponent there exists a neg-
Lyapunov vectors and orthonormal GS vectors, we con-
ative local GS exponent such that their pair sum
sider a simple two-dimensional example, the H´enon map
vanishes [34]:
[33],
x =a−x2 +by , (±)ΛGS(Γ)=−(±)ΛGS (Γ), (cid:96)=1,...,D. (22)
n+1 n n (cid:96) D+1−(cid:96)
y =x ,
n+1 n As indicated, such a symplectic local pairing sym-
metry is found both forward and backward in
with a = 1.4 and b = 0.3. The light-blue line in Fig. 1
time. Non-symplectic systems do not display
represents the H´enon attractor, which is known to coin-
that symmetry. On the other hand, the re-
cidewithitsunstablemanifold. Theblacklineisafinite-
orthonormalization process tampers with the ori-
time approximation of its stable manifold. Let O denote
entation and rotation of the GS vectors and de-
apointontheattractor. TheGSvectorsatthispointare
stroys all consequences of the time-reversal invari-
indicated by g (O) and g (O), where the former is also
1 2 ance property of the original motion equations.
covariant and identical to v (O) (parallel to the unsta-
1
ble manifold). As required, the covariant vector v2(O) • Local covariant exponents: During their construc-
is tangent to the stable manifold (which was determined tion [35], only the norm of the covariant perturba-
by a different method). For the computation of the co- tion vectors needs to be periodically adjusted for
variant vectors at O, it is necessary to follow and store practical reasons, but the angles between them re-
the reference trajectory and the GS vectors sufficiently main unchanged. This process effectively destroys
far into the future up to a point F in our example. In all information concerning the d-dimensional vol-
Fig. 1 the GS vectors at this point are denoted by g1(F) ume elements. Thus, no symmetries analogous to
and g2(F). Applying an algorithm by Ginelli et al. [35], Eq.(22)exist. Instead,there-normalizedcovariant
to which we shall return below, a backward iteration to vectorsfaithfullypreservethetime-reversalsymme-
the original phase point O yields the covariant vectors, tryoftheequationsofmotion,whichisreflectedby
v (O) in particular.
2 (∓)Λcov(Γ)=−(±)Λcov (Γ), (cid:96)=1,...,D, (23)
(cid:96) D+1−(cid:96)
5
regardless,whetherthesystemissymplecticornot. systems.. In all applications below appropriate reduced
Thismeansthatanexpandingco-movingdirection units will be used.
is converted into a collapsing co-moving direction To ease the notation, we shall in the following omit to
by an application of the time-reversal operation. indicate the forward-direction in time by (+), if there is
no ambiguity.
The set of all global Lyapunov exponents, ordered ac-
cording to size, is referred to as the Lyapunov spectrum.
For stationary Hamiltonian systems in thermal equilib-
III. TWO-DIMENSIONAL HARD-DISK AND
rium the (global) conjugate pairing rule holds, WCA FLUIDS IN EQUILIBRIUM
(±)λ =−(±)λ ,
(cid:96) D+1−(cid:96) We are now in the position to apply this formalism to
various models for simple fluids [41]. First we consider
whichisaconsequenceofEq.(22)andofthefactthatthe
two moderately dense planar gases, namely a system of
global exponents are the phase-space averages of these
elastic smooth hard disks with a pair potential
quantities. In such a case only the first half of the spec-
trumcontainingthepositiveexponentsneedstobecom- (cid:26)
∞, r ≤σ,
puted. In the following we shall refer to this half as the φHD = 0, r >σ,
positive branch of the spectrum.
The subspaces spanned by the covariant (or GS) vec- and a two-dimensional Weeks-Chandler-Anderson
tors associated with the strictly positive/negative global (WCA) gas interacting with a repulsive soft potential
exponents are known as the unstable/stable manifolds. [42, 43]
Both are covariant and are linked by the time-reversal
transformation, which converts one into the other. (cid:40)4(cid:15)(cid:104)(cid:0)σ(cid:1)12−(cid:0)σ(cid:1)6(cid:105)+(cid:15), r ≤21/6σ,
φ = r r
WCA
0, r >21/6σ.
D. Numerical considerations
InFig.2thepositivebranchesoftheir(global)Lyapunov
spectra are compared. Both gases consist of N = 400
The computation of the Gram-Schmidt exponents is
particles in a square box with side length L and periodic
commonly carried out with the algorithm of Benettin
boundaries, and have the same density, ρ=N/L2 =0.4,
et al. [24, 26] and Shimada et al. [25]. Based on
and temperature, T =1. Although, as expected, for low
this classical approach, a reasonably efficient algorithm
densities the Lyapunov spectra are rather insensitive to
for the computation of covariant exponents has been re-
the interaction potential, the comparison of Fig. 2 for
cently developed by Ginelli et al. [35, 36]. Some com-
putational details for this method may also be found in
Refs.[16,31,37]. AnalternativeapproachisduetoWolfe
and Samelson [38], which was subsequently applied to
Hamiltonian systems with many degrees of freedom [39]. 3.5
Theconsiderationsoftheprevioussectionarefortime-
00..44
3
continuous systems based on the differential equations
(1). They may be readily extended to maps such as 2.5 λλ
systems of hard spheres, for which a pre-collision state WWCCAA
HHDD
of two colliding particles is instantaneously mapped to 2 00
an after-collision state. Between collisions the particles λ 778800 779900 880000
1.5 ll
move smoothly. With the linearized collision map the
time evolution of the perturbation vectors in tangent 1
space may be constructed [6].
0.5
In numerical experiments of stationary systems, the WCA
HD
initialorientationoftheperturbationvectorsisarbitrary. 0
There exists a transient time, during which the pertur- 0 0.2 0.4 0.6 0.8 1
bation vectors converge to their proper orientation on l/2N
the attractor. All symmetry relations mentioned above
FIG. 2. (Global) Lyapunov spectra of a planar hard-disk gas
refer to this well-converged state and exclude transient
and of a planar WCA fluid with the same number of parti-
conditions [40]. cles,densityandtemperature. Fordetailswerefertothemain
For the computation of the full set of exponents, the text. Onlythepositivebranchesofthespectraareshown. Al-
reference trajectory and D perturbation vectors (each of thoughthespectraconsistofdiscretepointsforintegervalues
dimensionD)havetobefollowedintime,whichrequires oftheindexl,smoothlinesandreducedindexesl/2N areused
D(D+1) equations to be integrated. Present computer for clarity in the main panel. In the inset a magnified view
technology limits the number of particles to about 104 of the small-exponents regime is shown, for which l is not
for time-continuous systems, and to 105 for hard-body normalized. The figure is taken from Ref. [41].
6
1.0
0.8
0.6
W
1.0
0.4
W
0.2
0.5
373 396 420
ℓ
0.0
0 0.2 0.4 0.6 0.8 1
ℓ / 4N
FIG.4. LocalizationspectraW forthecompletesetofGram-
Schmidt vectors (blue) and covariant vectors (red) for 198
FIG.3. Localizationoftheperturbationvectorg1 inphysical hard disks in a periodic box with an aspect ratio 2/11. The
space(theredsquare)foraplanargasof102,400softrepulsive density ρ=0.7 and the temperature T =1. Reduced indices
disks [44]. The quantity µ1, measuring the contribution of (cid:96)/4N are used on the abscissa of the main panel. The in-
individualparticlestotheperturbationvectorassociatedwith set shows a magnification of the central part of the spectra
themaximumGSexponent,isplottedintheverticaldirection dominated by Lyapunov modes. From Ref. [37].
at the position of the particles.
of v1 varies as a negative power of N [18, 41, 45]. For
moderately dense gases already reveals a rather strong example, Fig. 3 is obtained for a system of 102,400 (!)
sensitivity. In particular, the maximum exponent λ is smooth disks interacting with a pair potential similar to
1
much lower for the WCA fluid, which means that de- the WCA potential [44].
terministic chaosis significantly reducedin systems with Thisdynamicallocalizationoftheperturbationvectors
smooth interaction potentials. This difference becomes associated with the large Lyapunov exponents is a very
even more pronounced for larger densities, as will be generalphenomenonandhasbeenalsoobservedforone-
shown below. dimensional models of space-time chaos [46] and Hamil-
The maximum (minimum) Lyapunov exponent de- tonian lattices [47].
notestherateoffastestperturbationgrowth(decay)and For (cid:96) > 1, the localization for v(cid:96) differs from that of
is expected to be dominated by the fastest dynamical g(cid:96). Thismaybeseenbyusingalocalizationmeasuredue
events such as a locally enhanced collision frequency. To to Taniguchi and Morriss [48, 49],
provethisexpectationonemayaskhowindividualparti-
N
cles contribute to the growth of the perturbations deter- W = 1 exp(cid:104)S(cid:105); S(Γ(t))=−(cid:88)µ((cid:96))lnµ((cid:96)),
mined by v(cid:96) or g(cid:96). Writing the normalized perturbation N i i
i=1
vectors in terms of the position and momentum pertur-
bations of all particles, v(cid:96) = {δq((cid:96)),δp((cid:96));i = 1,...,N}, which is bounded according to 1/N ≤ W ≤ 1. The
i i
the quantity lower and upper bounds correspond to complete local-
ization and delocalization, respectively. S is the Shan-
µ((cid:96)) ≡(cid:16)δq((cid:96))(cid:17)2+(cid:16)δp((cid:96))(cid:17)2 (24) non entropy for the ‘measure’ defined in Eq. (24), and
i i i (cid:104)...(cid:105) denotes a time average. The localization spectra
for the two sets of perturbation vectors are shown in
ispositive,bounded,andobeysthesumrule(cid:80)N µ((cid:96)) = Fig. 4 for a rather dense system of hard disks in two
i=1 i
1 for any (cid:96). It may be interpreted as a measure for the dimensions. One may infer from the figure that localiza-
activity of particle i contributing to the perturbation in tion is much stronger for the covariant vectors (red line)
question. Equivalent relations apply for g(cid:96). If µ(1) is whose orientations in tangent space are only determined
i
plotted - vertical to the simulation plane - at the posi- by the tangent flow and are not constrained by periodic
tion of particle i, surfaces such as in Fig. 3 are obtained. re-orthogonalization steps as is the case for the Gram-
Theyarestronglypeakedinasmalldomainofthephysi- Schmidt vectors (blue line). One further observes that
cal space indicating strong localization of v1 ≡ g1. It thelocalizationofthecovariantvectorsbecomeslessand
means that only a small fraction of all particles con- less pronounced the more the regime of coherent Lya-
tributes to the perturbation growth at any instant of punov modes, located in the center of the spectrum, is
time. This property is very robust and persists in the approached.
thermodynamic limit in the sense that the fraction of NextweturnourattentiontothemaximumLyapunov
particles, which contribute significantly to the formation exponent λ and to the Kolmogorov-Sinai entropy. Both
1
7
quantities are indicators of dynamical chaos. The KS-
entropy is the rate with which information about an ini-
tial state is generated, if a finite-precision measurement
HHDD
at some later time is retraced backward along the sta- 55..00
ble manifold. According to Pesin’s theorem [50] it is
equal to the sum of the positive Lyapunov exponents, 44..00
hKS = (cid:80)λ(cid:96)>0λ(cid:96). It also determines the character- λλ11 33..00 WWCCAA
istic time for phase space mixing according to [51–53]
τmix =1/h . 22..00
KS
In Fig. 5 we compare the isothermal density depen-
dence of the maximum Lyapunov exponent λ (top 11..00
1
panel), of the smallest positive exponent λ (mid-
2N−3 00
dle panel), and of the Kolmogorov-Sinai (KS) entropy 00 00..22 00..44 00..66 00..88 11 11..22
perparticleh /N (bottompanel)forplanarWCAand
KS
HD fluids. Both systems contain N = 375 particles in
a periodic box with aspect ratio 0.6 and with a temper-
00..66
ature T = 1. As expected, these quantities for the two
HHDD
fluids agree at small densities, but differ significantly for 00..55
large ρ.
For hard disks, van Zon and van Beijeren [54] and de N-3N-3 00..44
Wijn [55] used kinetic theory to obtain expressions for λλ22 00..33
λ and h /N to leading orders of ρ. They agree very
1 KS
well withcomputer simulations ofDellago et al. [56,57]. 00..22
The regime of larger densities, however, is only qualita-
00..11 WWCCAA
tively understood. λHD and hHD/N increase monoton-
1 KS
ically due to the increase of the collision frequency ν. 00
00 00..22 00..44 00..66 00..88 11 11..22
The(first-order) fluid-solidtransitionshows upas astep
between the freezing point of the fluid ( ρHD = 0.88 )
f
andthemeltingpointofthesolid(ρHD =0.91[58])(not 77..00
m HHDD
shown in the figure). These steps disappear, if instead 66..00
of the density the collision frequency is plotted on the
55..00
abscissa. λHD and hHD/N diverge at the close-packed
density due1to the divKerSgence of ν. /N/NSS 44..00
For the WCA fluid, both the maximum exponent and hhKK WWCCAA
33..00
the Kolmogorov-Sinai entropy have a maximum as a
22..00
functionofdensity,andbecomeverysmallwhentheden-
sity of the freezing point is approached, which happens 11..00
for ρWCA = 0.89 [59]. This behavior is not too surpris-
f 00
ing in view of the fact that a harmonic solid is not even 00 00..22 00..44 00..66 00..88 11 11..22
chaotic and λ vanishes. The maximum of λWCA occurs ρρ
1 1
at a density of about 0.9ρWCA and confirms earlier re-
f
sults for Lennard-Jones fluids [60]. At this density chaos FIG.5. IsothermaldensitydependenceofthemaximumLya-
is most pronounced possibly due to the onset of coop- punov exponent, λ (top), of the smallest positive exponent,
1
erative dynamics characteristic of phase transitions. At λ (middle), and of the Kolmogorov-Sinai entropy per
2N−3
about the same density, mixing becomes most effective particle, hKS/N (bottom), for hard and soft-disk systems.
as is demonstrated by the maximum of the KS-entropy. From Ref. [41].
The latter is related to the mixing time in phase space
according to τmix =1/h [51, 52].
KS
It is interesting to compare the Lyapunov time τ = related collisions). For densities ρ > 0.4 the Lyapunov
λ
time is progressively larger than the collision time and
1/λ , which is a measure for the system to “forget” its
1
higher-order corrections such as ring collisions become
past, with the (independently-measured) time between
collisionsofaparticle,τ =1/ν. InFig.6suchacompar- important. Also the lines for τλ and τKS cross even be-
c
fore, at a density 0.1. This is a consequence of the shape
ison is shown for a three-dimensional system of 108 hard
change of the Lyapunov spectrum with density [61]: the
spheresinacubicboxwithperiodicboundaryconditions
[53, 61]. Also included is the behavior of τ ≡N/h . smallpositiveexponentsgrowfasterwithρthanthelarge
KS KS
exponents. ItalsofollowsthatneithertheLyapunovtime
Forsmalldensities,wehaveτ (cid:28)τ ,andsubsequentcol-
λ c
nor the mixing time determine the correlation decay of -
lisions are uncorrelated. This provides the basis for the
say-theparticlevelocities. Forlowerdensities,forwhich
validity of lowest-order kinetic theory (disregarding cor-
8
106 IV. LYAPUNOV MODES
104 τc ρ: 0.1 0.4 1.0 Already the first simulations of the Lyapunov spectra
forharddisksrevealedaninterestingnewstep-likestruc-
τKS
ture of the small positive and negative exponents in the
102 τλ center of the spectrum [6], which is due to degenerate
τ exponents. A similar structure was also found for hard
100 dumbbellsystems[1],towhichwecomebackinSec.VB.
The explanation for this behavior lies in the fact that
the perturbation vectors for these exponents are charac-
10-2
terized by coherent sinusoidal patterns spread out over
the whole physical space (the simulation cell). We have
10-4 referred to these collective patterns as Lyapunov modes.
10-5 10-3 10-1 101 103
The modes are interpreted as a consequence of a sponta-
ν
neousbreakingofcontinuoussymmetriesand,hence,are
intimately connected with the zero modes spanning the
FIG.6. (Coloronline)ComparisonoftheLyapunovtimeτ =
λ centralmanifold[62,63]. Theyappearassinusoidalmod-
1/λ ,thetimeτ ≡N/h ,andthecollisiontimeτ =1/ν,
1 KS KS c ulations of the zero modes in space – with wave number
asafunctionofthecollisionfrequencyν, forasystemof108
k (cid:54)= 0 – due to the spontaneous breaking of the contin-
hard spheres in a cubic box with periodic boundaries. The
vertical lines indicate collision frequencies corresponding to uous symmetries. For k → 0, the modes reduce to a
the densities ρ=0.1,0.4, and 1.0. linear superposition of the zero modes, and their Lya-
punov exponent vanishes. The experimentally observed
wave vectors depend on the size of the system and the
ring collisions are not important, the decay of the veloc- nature of the boundaries (reflecting or periodic).
ity autocorrelation function z(t) is strictly dominated by OurdiscoveryoftheLyapunovmodestriggeredanum-
the collision time. In Ref. [53] we also demonstrate that ber of theoretical approaches. Eckmann and Gat were
forhard-particlesystemsthetimeτ doesnotprovidean the first to provide theoretical arguments for the ex-
λ
upperbound forthe time for which correlation functions istence of the Lyapunov modes in transversal-invariant
may be reliably computed. systems [64]. Their model did not have a dynamics in
For later reference, we show in Fig. 7 the time- phase space but only an evolution matrix in tangent
dependence of the local exponents for (cid:96) = 1 and (cid:96) = space, which was modeled as a product of independent
D =4N =16ofasystemconsistingoffoursmoothhard random matrices. In another approach, McNamara and
disks. The symplectic symmetry as given by Eq. (22) is Mareschal isolated the six hydrodynamic fields related
well obeyed. to the invariants of the binary particle collisions and the
vanishing exponents, anduseda generalizedEnskogthe-
orytoderivehydrodynamicevolutionequationsforthese
250
fields. Their solutions are the Lyapunov modes [65]. In
200
amoredetailedextensionofthisworkrestrictedtosmall
150 ℓ = 1
densities, a generalized Boltzmann equation is used for
100 thederivation[66]. deWijnandvanBeijerenpointedout
50 the analogy to the Goldstone mechanism of construct-
St() 0 ing (infinitesimal) excitations of the zero modes and the
G
Λ Goldstone modes of fluctuating hydrodynamics [67, 68].
-50
They used this analogy to derive the modes within the
-100
framework of kinetic theory [69]. With a few exceptions,
-150 ℓ = 16
a rather good agreement with the simulation results was
-200
achieved, at least for low densities. Finally, Taniguchi,
-250 Dettmann,andMorrissapproachedtheproblemfromthe
1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000
point of view of periodic orbit theory [70] and master
t
equations [71].
Themodeswereobservedforhard-ballsystemsinone,
FIG. 7. (Color online) Test of the symplectic symmetry of two, and three dimensions [48, 49, 61, 62, 72, 73], for
Eq. (22) for a smooth hard-disk system with N = 4. The planar hard dumbbells [18], and also for one- and two-
phase space has 16 dimensions. The local GS exponents for dimensional soft particles [41, 74, 75].
(cid:96) = 1 and (cid:96) = 16 are plotted as a function of time t along a A formal classification for smooth hard-disk systems
short stretch of trajectory.
has been given by Eckmann et al. [73]. If the 4N com-
9
ponents of a tangent vector are arranged according to For the same k, the values of the Lyapunov expo-
nents for the longitudinal and momentum modes coin-
δΓ = (δqx,δqy;δpx,δpy), cide. These modes actually belong to a combined sub-
space LP(n) ≡ L(n) ⊕ P(n). The dimensions of the
the six orthonormal zero modes, which span the central subspaces T(n) and LP(n) determine the multiplicity of
manifold and are the generators of the continuous sym- the exponents associated with the T and LP modes. For
metry transformations, are given by [45, 73] a periodic boundary, say in x direction, the multiplicity
is 2 for the T, and 4 for the LP modes. For details we
1
e1 = √ (px,py; 0,0), refer to Ref. [73].
2K
Thetransversemodesarestationaryinspaceandtime,
1
e = √ (1,0; 0,0), but the L and P modes are not [73, 76, 77]. In the fol-
2
N lowing, we only consider the case of periodic boundaries.
1 Any tangent vector from an LP subspace observed in a
e = √ (0,1; 0,0),
3 N simulation is a combination of a pure L and a pure P
1 mode. The dynamics of such an LP pair may be rep-
e4 = √ (0,0; px,py), resented as a rotation in a two-dimensional space with
2K
a constant angular frequency proportional to the wave
1
e = √ (0,0; 1,0), number of the mode,
5
N
1 ω =vk ,
e = √ (0,0; 0,1). n n
6
N
wheretheLmodeiscontinuouslytransformedintotheP
modeandviceversa. SincethePmodeismodulatedwith
e correspondstoashiftofthetimeorigin,e toachange
1 4 the velocity field of the particles, the experimentally ob-
ofenergy,e ande toauniformtranslationinthexand
2 3 served mode pattern becomes periodically blurred when
y directions, respectively, and e and e to a shift of the
5 6 its character is predominantly P. Since during each half
total momentum in the x and y directions, respectively.
ofaperiod,duringwhichallspacialsineandcosinefunc-
The six vanishing Lyapunov exponents associated with
tionschangesign,themodeisoffsetinthedirectionofthe
these modes have indices 2N −2 ≤ i ≤ 2N +3 in the
wave vector k, which gives it the appearance of a travel-
center of the spectrum. Since, for small-enough k, the
ing wave with an averaged phase velocity v. If reflecting
perturbation components of a particle obey δp = Cδq
boundaries are used, standing waves are obtained. The
for the unstable perturbations (λ > 0), and δq = −Cδp
phase velocity v is about one third of the sound velocity
for the stable perturbations (λ < 0), where C > 0 is a
[76], but otherwise seems to be unrelated to it.
number,onemayrestricttheclassificationtotheδqpart
ThebasisvectorsspanningthesubspacesoftheTand
of the modes and, hence, to the basis vectors e ,e ,e
1 2 3 LP modes may be reconstructed from the experimental
[73]. Now the modes with a wave vector k may be clas-
datawithaleastsquaremethod[37,73]. Asanexample,
sified as follows:
weconsidertheL(1,0)modeofasystemconsistingof198
hard disks at a density 0.7 in a rectangular periodic box
1. Transverse modes (T) are divergence-free vec-
withanaspectratio2/11. TheLP(1,0)subspaceincludes
torfieldsconsistingofasuperpositionofsinusoidal
the tangent vectors for 388 ≤ (cid:96) ≤ 391 with four identi-
modulations of e and e .
2 3
cal Lyapunov exponents. The mutually orthogonal ba-
2. Longitudinal modes (L) are irrotational vector sis vectors spanning the corresponding four-dimensional
fields consisting of a superposition of sinusoidal subspaceareviewedasvectorfieldsinpositionspaceand
modulations of e and e . are given by [73]
2 3
3. Momentum modes (P) are vector fields consist- (cid:18)1(cid:19) (cid:18)2πq (cid:19) (cid:18)1(cid:19) (cid:18)2πq (cid:19)
ing of sinusoidal modulations of e1. Due to the L(1,0): 0 cos L x , 0 sin L x ,
x x
random orientations of the particle velocities, a P
modeisnoteasilyrecognizedasfundamentalmode
in a simulation. However, it may be numerically P(1,0):(cid:18)px (cid:19)cos(cid:18)2πqx(cid:19), (cid:18)px (cid:19)sin(cid:18)2πqx(cid:19),
transformed into an easily recognizable periodic py Lx py Lx
shape [16, 73], as will be shown in Fig. 9 below.
where,forsimplicity,wehaveomittedthenormalization.
The subspaces spanned by these modes are denoted by
IftheGram-Schmidtvectorsareusedforthereconstruc-
T(n), L(n), and P(n), respectively, where the wave
tion, the components corresponding to the cosine are
vector for a periodic box with sides L ,L is k =
x y shownforL(1,0)inFig.8,forP(1,0)inFig.9. Inthelat-
(2πn /L ,2πn /L ), and n ≡ (n ,n ). n and n are
x x y y x y x y tercase, themodestructureisvisibleduetothedivision
integers.
by p and p as indicated. This also explains the larger
x y
scattering of the points. The components proportional
10
0.1 0.1
x
p
δqx 0.0 / x 0.0
q
δ
-0.1 -0.1
-20 -10 0 10 20 -20 -10 0 10 20
q q
x x
0.1 0.1
y
p
δqy 0.0 / y 0.0
q
δ
-0.1 -0.1
-20 -10 0 10 20 -20 -10 0 10 20
q q
x x
FIG. 8. (Color online) Reconstructed L(1,0) mode compo- FIG. 9. (Color online) Reconstructed P(1,0) mode compo-
nents proportional to cos(2πq /L ). The system consists of nentsproportionaltocos(2πq /L )forthesameLPsubspace
x x x x
198smoothharddisksinaperiodicbox. Fordetailswerefer as in Fig. 8.
to the main text.
give analytical expressions for the modes in accordance
to the sine are analogous, but are not shown. A related with the boundary conditions applied. This procedure
analysis may be carried out for the covariant vectors. works for the stationary transverse modes, as well as for
Recently, Morriss and coworkers [77–79] considered in the time-dependent LP modes. For the latter, the par-
detailthetangent-spacedynamicsofthezeromodesand tialdifferentialequationsforthecontinuousperturbation
of the mode-forming perturbations over a time τ, during fields assume the form of a wave equation with traveling
which many collisions occur. In addition, they enforced wave solutions [77]. For quasi-one-dimensional systems
the mutual orthogonality of the subspaces of the modes (withnarrowboxessuchthatparticlesmaynotpasseach
andoftheir conjugate pairs(forwhichthe Lyapunovex- other), the phase velocity is found to depend on the par-
ponentsonlydifferbythesign),withthecentralmanifold ticle density via the mean particle separation and the
byinvokingan(inverse)Gram-Schmidtprocedure,which mean time between collisions (of a particle). For fully
starts with the zero modes and works outward towards two-dimensional models of hard disks, the wave veloc-
modes with larger exponents. With this procedure they ity becomes proportional to the collision frequency of a
were able to construct the modes in the limit of large τ particle, and inversely proportional to the mean squared
andtoobtain(approximate)valuesfortheLyapunovex- distance of a particle to all its neighbors in the first co-
ponentsforthefirstfewLyapunovmodes(countedaway ordination shell [77, 79]. This is an interesting result,
fromthenullspace)[78]. Theagreementwithsimulation since it provides a long-looked-for connection between a
results is of the order of a few percent for the T modes, mode property – the wave velocity – with other density-
and slightly worse for the LP modes. dependent microscopic quantities of a fluid.
For the modes with a wave length large compared
to the inter-particle spacing, a continuum limit for the
perturbation vectors leads to a set of partial differen-
tialequationsfortheperturbationfields,whosesolutions