The number dependence of the maximum Lyapunov exponent

PHYSICA ELSEVIER

Physica A 240 (1997) 96-104

The number dependence of the maximum Lyapunov exponent D e b r a J. Searles a, Denis J. Evans a'*, Dennis J. Isbister b "Research School of Chemistry, Australian National University, GPO Box 414, Canberra, ACT 2601, Australia bDepartment of Physics, University College, University of NSW, ADFA, Canberra, ACT 2600, Australia

Abstract

We examine the number dependence of the largest Lyapunov exponent for nonlinear dynamical systems in one, two and three Cartesian dimensions. Our results suggest that the largest Lyapunov exponent diverges logarithmically with system size, independently of the number of Cartesian dimensions and the interaction potential.

The maximum Lyapunov exponent characterises the maximum dynamical instability present in a dynamical system. It can be defined as the maximum rate of exponential separation of two trajectories in phase space which are generated by the same equations of motion but which originate from two infinitesimally close initial conditions [1]. For an ergodic system, this rate of separation is (except for a set of measure zero) independent of the initial conditions. A dynamically unstable system has exponentially diverging trajectories and therefore has a positive maximal Lyapunov exponent as the first exponent of the set or spectrum of Lyapunov exponents. In this paper we consider the dependence of the maximum Lyapunov exponent, 2 . . . . on the number of particles for systems in one, two and three Cartesian dimensions. The dynamical systems which we examine include the classical Fermi-Pasta-Ulam coupled oscillator problem (in its fl form), as well as the traditional molecular dynamics of N particles acting under short ranged forces in two and three dimensions. In the 1960s, Lebowitz and Percus [2] showed that corrections to the microscopic distribution functions and macroscopic thermodynamic functions for systems at equilibrium were of order l/N, where N is the number of particles in the system under

*Corresponding author. 0378-4371/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PII S 0 3 7 8 - 4 3 7 1 ( 9 7 ) 0 0 1 3 3 - 7

D.J. Searles et al./Physica A 240 (1997) 96-104

97

study. This convergence of the thermodynamic properties to the so-called thermodynamic limit (in which the volume V ~ oo, and number of particles N ~ oo but the ratio, the number density, is fixed n = N / V ) was later verified using simulations in different ensembles. For an incompressible viscous fluid obeying the Navier Stokes momentum equation, Ruelle [3] found that in two Cartesian dimensions the sum of the positive Lyapunov exponents is bounded above by a fixed multiple of the total viscous dissipation. It has subsequently been assumed for systems in which the thermodynamic properties are well behaved in the thermodynamic limit, that the largest Lyapunov exponent is finite in the same limit. Recent exact results linking Lyapunov exponents to the thermophysical properties of the system have only reinforced this assumption. In particular, the conjugate pairing rule of Evans et al. [4] states that a transport coefficient such as electrical conductivity is proportional to the sum of the maximal Lyapunov exponents divided by the square of the electric field which drives the electric current. If the linear (small field) transport coefficient is finite in the thermodynamic limit, then the sum of the maximal Lyapunov exponents divided by the square of the driving thermodynamic force must also be finite in the thermodynamic limit. It has been widely assumed that not only is the sum of the maximal Lyapunov exponents finite but the maximal Lyapunov exponents are individually finite in the thermodynamic limit. The escape rate formalism of Gaspard and coworkers [5,6] gives an exact relation between transport coefficients and Lyapunov exponents: again the results seems to endorse the assumption that if the transport coefficient is finite in the thermodynamic limit, the corresponding Lyapunov exponents themselves must also be convergent. It has been suggested that due to the so-called long time tail effect, the Navier Stokes transport coefficients diverge in two dimensions [7]. This observation has lead some investigators to study the number dependence of the largest Lyapunov exponent in two dimensional systems. Stoddard and Ford [8] used arguments based on kinetic theory to predict a logarithmic number dependence of "~maxfor two dimensional systems. Molecular dynamics simulations by Posch and Hoover [9] suggested that )-max increases logarithmically with N for systems confined by corrugated boundaries both at equilibrium and undergoing boundary-driven shear flow. Subsequently Hoover and Posch [10] found no definite evidence of a logarithmic size dependence when Lees Edwards shearing periodic boundary conditions were used. In their studies of the one dimensional Fermi-Pasta Ulam system Livi et al. [11] concluded in agreement with previous studies [12] that, "it must be registered that the exponential divergence of nearby trajectories is independent of N (the number of coupled oscillators) in the limit of large N ( ~ 40)". We are unaware of any previous studies of the N dependence of the largest Lyapunov exponent for systems in three or more Cartesian dimensions. In the present paper we provide computer simulation results which suggest that for systems in one, two or three Cartesian dimensions, the largest Lyapunov exponent diverges logarithmically with the number of particles in the system. Each of the systems we study employs periodic boundary conditions.

D.J. Searles et al./Physica A 240 (1997) 96 104

98

The one dimensional system investigated in this work is the Fermi-Pasta-Ulam fl (FPU) model [11], which represents the dynamics of a one-dimensional chain of N non-linearly coupled oscillators. The Hamiltonian for this system is 2 p { ~-"~(X i - - X I + I ) 2 +

H =

(X i --

X/+1)4

(1)

i=1

where x~ is the displacement relative to the equilibrium position of oscillator i and Pi is the momentum conjugate to xi. Livi et al. [11] integrated the resulting equations of motion Xi -~" Pi , Pi = F i = (Xi+ l - 2 x i + x i - 1 ) -~- I [ ( X i + 1 - - Xi) 3 -1- (Xi -- X I - 1 ) 3]

(2)

with the Verlet leap-frog algorithm. Throughout this work we use the Gear predictor-corrector (4th order) algorithm with small time steps (0.001 ~
VwcA(rij)

{4(~ ----

~ ) + 1,

O,

rij < 21/6 ,

(3) rij > 21/6 ,

where rij is the particle separation. In both the two- and three-dimensional systems the temperature was kept constant using a Gaussian thermostat [14]. Both WCA systems were studied at equilibrium. In addition, the three-dimensional WCA system was investigated away from equilibrium under the influence of shear while being thermostatted using standard nonequilibrium molecular dynamics techniques. The equations of motion of the two- and three-dimensional systems employed the thermostatted SLLOD equations of motion which give an exact description of adiabatic shear flow arbitrarily far from equilibrium [14], ¢1i -~ Pl + iTyi, Pi = F i - iTpyi -- ePi •

(4)

The qi are the particle Cartesian position vectors, Pi are the conjugate momenta, F~ are the forces due to the interaction potential of Eq. (3), 7 is the strain rate and ~ is the thermostat multiplier which maintains constant temperature. The thermostat multiplier is given by = ~V=IFi'PIN --2YPxiPyi ~'~ i = lPi

(5)

D.,L Searles et a l . / P h y s i c a A 240 (1997) 96 104

99

The equilibrium equations of motion and thermostat multiplier are obtained from Eq. (4) and (5) with the strain rate, 7, set to zero, 7 = 0. The simulations for the systems in two and three Cartesian dimensions were carried out at a temperature T = y~ p2/gN = 1.0 where 9 the number of Cartesian dimensions. The number density, n = N / V , was 0.8. "~maxis given by the rate of exponential growth in the separation of two infinitesimally displaced trajectories known as the mother (the fiducial) and daughter (tangent} trajectories. In a simulation, the separation is small but finite and therefore must be rescaled periodically to ensure the separation remains small throughout the simulation [15]. In this work, a continuous rescaling method was used [16]. The mother trajectory, F(t), is determined using the equations of motion above whereas the equations of motion of the daughter trajectory, F1 (t), are modified to ensure that the separation of the mother and daughter trajectories remains constant. The equation of motion of the constrained tangent vector, 6 C= F1 - F where F1 is the phase-space vector of the constrained daughter trajectory, is given by 6c = T ( F ) - 6 - f6

(6)

where T is the stability matrix for the equations of motion (0['/OF = T(F)). The constraint multiplier f maintaining the separation 6 c is given by 6,.T.6C ~ - i ] ~, ~ , , ii

(7)

The Lyapunov exponent can be obtained from the time average of ~ (derived in [16]), t

;m.,

,4~,

,im'f ds f(s).

~I1~(0) 11 - t ~ t

(8)

0

Fig. 1 shows the dependence of ~max o n toglo N for the F P U system with smaller energy density E / N = 1.0. It is clear that 2m,~ is consistent with a logarithmic variation with N for sufficiently large N. The function )~max= a + b / N + c log~o N (a = 0.05005( + 0.00045), b = - 0.2134( + 0.0097), c = 0.00239( + 0.00012) is shown in the figure) gives a good fit to the data. A logarithmic dependence of 2m~x was also found for the same system at a higher energy density E / N = 26.4. The data plotted in Fig. 1 was also examined to show the convergence of the time average shown in (8). Generally the data approached its long time limit with a lit time dependence. The rate of convergence and the length of the simulation runs is consistent with the estimated statistical uncertainties denoted with the error bars in Fig. 1. At the expense of introducing an extra variable, the data shown in Fig. 1 can be fitted almost as well by the function "{max= a' + b'/N + c'N d' (a' = - 0.21798 + 5.4, b' = - 0.21442 ± 0.017, c' = 0.32214 + 5.4, d' -- 0.0031325 ___0.05). However, the extremely large uncertainities in the fitted coefficients indicates that there is insufficient structure in the data to determine the coefficients. For large N, the two functions have

100

D.J. Searles et al./Physica A 240 (1997) 96-104

0.065

'

'

'

' ' ' ' ' 1

'

'

'

' ' ' ' ' 1

'

'

'

'

' ' ' ' ' 1

'

'

'

''"1

.

.

...~"

0.060

0.055

1 0

max

/

0.050

O

/ / /

0.045

¢ ,

0.040 10

100

.

.

,

, , , . I

1000 N

,

,

,

,

. . . .

I

10000

.

.

.

.

.

,I

100000

Fig. 1. The dependence Of2maxof log10 N for the FPU systemwith E/N = 1.0 determined using molecular dynamics simulations. The crosses are data calculated in this work and the circles were obtained by Livi et al. [11]. The curve is a best fit of the function described in the text, L.ax= 0.05005(+0.00045) - 0.2134( +0.0097)/N + 0.00239( +0.00012)loglo N.

a similar variation with N since d', which is positive, is very close to zero and d2maJdN ~ c ' d ' N n ' - i ~ (c/ln 10)/N and c'd' ~ c / l n 10. There are large correlations in the uncertainties of the coefficients partly because the data only determines the product, c'd', and does not individually determine c' and d'. While the two functions are essentially identical at N = 1000, for Avogadro's number of particles the difference between the two functions is less than 25%. Using computer simulations alone it will be impossible to distinguish between the two functions. Since d' is positive, this power law dependence is still consistent with the divergence of the largest L y a p u n o v exponent in the thermodynamic limit. Fig. 1 also compares the 2maxdata obtained in this work with that calculated by Livi et al. [11] (LPR). The two sets of calculations are consistent but L P R only studied systems where N ~< 320. L P R calculated 2max for systems of various energy densities and concluded that the L y a p u n o v spectrum converges in the thermodynamic limit (N/> O(102)) and that j-max is finite in the thermodynamic limit. They also suggest that the convergence is more rapid for the larger energy densities. Replotting the data of L P R in the manner shown in Fig. 1 indicates that their systems were simply too small to observe the logarithmic divergence. In Fig. 2 we examine the dependence of J'max on loglo N for the two-dimensional system. Again the data are consistent with a logarithmic divergence with N. The smooth curve through the data is a fit to the function, )-max= a + b / N + c loglo N, where a, b, c are parameters. Fig. 2 also shows the corresponding dependence of the

D.J. Searles et al./Physica A 240 (1997) 96-104

3.65

'

•

'

6.32

'''''1

3.64 3.63 p

6.31

A

i

max

3.62

101

i

6.3 max

,

/

i i

3.61

6.29

A

6.28

3.6 6.27

3.59

/

3.58

"~ ......

6.26

'~- -'~-- --o ..........

Z,

3.57 10

100

6.25 1000

l04

105

N

Fig. 2. The dependence of "~max(crosses) and pressure (squares) on loglo N for WCA particles in two Cartesian dimensions determined using molecular dynamics simulations. The simulation was carried out at T = 1.0 and n = 0 . 8 . The curves are the best fits to the functions described in the text: "~-max = 3.5286(+0.00062) - 0.38 ( +0.24)/N + 0.0278( +0.0018)logto N for the Lyapunov exponent (crosses), and, p = 6.2565(+0.0001) + 3.18( _+0.19)/N for the pressure (squares).

pressure on log10 N. The curve through the pressure data is a fit to the function, p = d + e/N, where d, e are parameters. Comparison of the two curves shows that for N ~< ~ 900, the pressure decreases until it has relaxed to its asymptotic value. For N >I ~ 900, the pressure is essentially independent of N. The variation of the largest Lyapunov exponent is completely different. For N / > ~ 50, "~max increases logarithmically. It continues to increase logarithmically even after the pressure has relaxed to its asymptotic value. Fig. 3 shows the dependence of '~max on log10 N for the three-dimensional system at equilibrium. A logarithmic variation is observed. Fig. 3 also shows the dependence of the pressure on log10 N for the system. The smooth curves through the data are fits to the same functional forms as used in the two-dimensional case. We have also examined the number dependence of the largest Lyapunov exponent for this system under shear but we do not display these results here. From all our results a consistent picture emerges: after thermophysical properties such as the pressure and the stress have relaxed to their asymptotic, large system values, the largest Lyapunov exponent continues to increase - apparently logarithmically. This logarithmic divergence appears to be independent of either the number of Cartesian dimensions or the form of the particle interaction potential. The results indicate that the behaviour of the thermodynamic properties have no direct bearing on the convergence of the maximum Lyapunov exponent in the thermodynamic limit,

D.J. Searles et al./Physica A 240 (1997) 96-104

102

4.512

&

6.63 i

4.51

6.625 ~

4.508

P

•,

4.506

/ ~ax/

6.62

~'max

6.615 4.504 6.61 4.502 6.605

4.5

6,6

4.498 J i

4.496 100

i

i

i

I

1000

6.595 104

N Fig. 3. The dependence of 2~ax and p on loglo N for W C A particles in three Cartesian dimensions determined using molecular dynamics simulations at equilibrium at T = 1.0 and n = 0.8. The curves are the best fits to the functions described in the text: 2 m a x = 4 . 4 7 4 9 ( + 0 . 0 0 1 9 ) - - 0 . 3 1 0 ( + 0 . 0 8 9 ) / N +0.00922( ___0.00057)logloN for the Lyapunov exponent (crosses), and, p = 6.6007(+0.0004) + 3.10( ___0.12)/N for the pressure (squares).

as has been previously assumed. Even if the variation of the largest exponent with N is not logarithmic, our data shows that its rate of "convergence" is vastly slower than is the convergence of thermodynamic properties in the thermodynamic limit. The divergence of the maximum Lyapunov exponent implies that the Kolmogorov-Sinai entropy per degree of freedom [1], must also diverge in the thermodynamic limit. In a nonequilibrium steady state the sum of all the Lyapunov exponents is proportional to the transport coefficient [4]. In order to obtain a finite transport coefficient in the linear (small field) thermodynamic limit (the order of the limits is important: the field must be taken to zero before we take the thermodynamic limit) the divergence of the maximum Lyapunov exponent must be cancelled by a corresponding divergence of the minimum Lyapunov exponent. Similar results are to be expected for each pair of Lyapunov exponents in systems which satisfy the conjugate pairing

rule. We shall now make some remarks about why this logarithmic divergence is difficult to observe in actual calculations. Firstly, the divergence is extremely weak. If we increase N from say 1000 to Avogadro's number, the largest Lyapunov exponent for the system shown in Fig. 1 2max merely doubles. Secondly, it is clear that the combination of two systems without coupling will result in a system with a Lyapunov spectrum which is simply the union of the two sets of Lyapunov exponents of the isolated systems. When the two systems are coupled, this degeneracy is broken and the present work suggests that the maximum Lyapunov

D.J. Searles et al./Physica A 240 (1997) 96 104

103

exponent increases as the logarithm of the number of degrees of freedom. There are two simple arguments that one can make from this observation. A spurious effect can occur if one does not allow the system to relax over a "sufficiently long time" before measuring the asymptotic rate of separation of the trajectories: this "sufficiently long time" must be greater than the time it takes for the distant parts of a large system to become "aware" of their mutual existence. The sound speed is the fastest communication that is possible in the systems studied in this work. Thus Lyapunov exponents should be measured after times t, =- ( N / n ) l / d / c , where c is the sound speed, n the number density and d the number of dimensions. If the system is simulated for times shorter than tc, the system cannot "know" how large it really is; distant parts of the system will be effectively uncoupled, resulting in a degenerate Lyapunov spectrum with a suppressed logarithmic divergence. Likewise for inaccurate calculations, the presence of noise or error can effectively s c r e e n long range couplings. This could then also result in spurious degeneracy and suppression of the logarithmic divergence in the Lyapunov spectrum. There is a relation between the accuracy with which you attempt to determine the Lyapunov exponent and the time over which the system must be studied. To resolve the largest Lyapunov exponent to a relative accuracy, re, requires using our method [17], a simulation over a time ts >>. - ln(re)/(2m,x re). In Fig. 1 our largest claimed error is about 1%. This implies a minimum simulation time, ts ~ 7500. With the exception of N = 100000, all the simulations shown in Fig. 1 were for times much longer than this. For example for N = 2560 ts ~ 250000. Even for N = 100000 the simulation time was 4200. As argued above, if the simulations were carried out for too short a time period one would expect the divergence to be suppressed. Finally we would like to make two further remarks. Firstly in a time averaged sense the largest Lyapunov exponent measures the greatest dynamical instability that is available to a dynamical system. The other exponents measure lesser instabilities. The eigenvector associated with the largest exponent seeks out that direction in tangent space for which dynamical instability is greatest. This exponent does not measure the a v e r a g e instability. As one increases the number of degrees of freedom in a system one increases the number of exponents and one should not be surprised that when one adds ever more particles to a system, one is presented with ever more opportunities for s u s t a i n i n g ever more unstable processes thereby leading to an ever increasing value for the largest Lyapunov exponent. A consequence of this possible divergence of the largest Lyapunov exponent is that the time over which dynamical trajectories are reversible (to any specified level of accuracy) must decrease to zero in the thermodynamic limit! This is because the reversal time, t . . . . varies as [18] trey ~--ln(6F)/)~min, where 2 m i n is the smallest Lyapunov exponent of the system and 6F is the trajectory error for a typical timestep. At equilibrium ~tml n = - - )~max and trey ~ + ln(6F)/2 .... Therefore although the Universe is governed by reversible equations of motion our work suggests that if the Universe is regarded as being infinite, actual reversal of Universal motion is impossible for any finite time no matter how short.

104

D.J. Searles et al. / P h y s i c a A 240 (1997) 9 6 - 1 0 4

Acknowledgements D J S w o u l d like to t h a n k t h e A u s t r a l i a n R e s e a r c h C o u n c i l for t h e a w a r d o f a n A R C F e l l o w s h i p . T h e a u t h o r s w o u l d like to t h a n k Dr. G . P . M o r r i s s for useful c o m m e n t s a b o u t this w o r k .

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [-10] [11] [12]

[13] [14] [15] [16] [17]

[18]

J.-P. Eckmann, D. Ruelle, Rev. Mod. Phys. 57 (1985) 617. J.L. Lebowitz, J.K. Percus, Phys. Rev. 124 (1961) 1673. D. Ruelle, Comm. Math. Phys. 87 (1982) 287. D.J. Evans, E.G.D. Cohen, G.P. Morriss, Phys. Rev. A 42 (1990) 5990; C.P. Dettmann, G.P. Morriss, Phys. Rev. E 53 (1996) 5545. P. Gaspard, G. Nicolis, Phys. Rev. Lett. 65 (1990) 1693. P. Gaspard, F. Baras, Phys. Rev. E 51 (1995) 5332. W.G. Hoover, W.T. Ashurst, R.J. Olness, J. Chem. Phys. 60 (1974) 4043; D.J. Evans, G.P. Morriss, Phys. Rev. Lett. 51 (1983) 1776; G.P. Morriss, D.J. Evans, Phys. Rev. A 39 (1989) 6335. S.D. Stoddard, J. Ford, Phys. Rev. A 8 (1973) 1504. H.A. Posch, W.G. Hoover, Phys. Rev. A 39 (1989) 2175. W.G. Hoover, H.A. Posch, Phys. Rev. E 51 (1995) 273. R. Livi, A. Politi, S. Ruffo, J. Phys. A 19 (1986) 2033; R. Livi, A. Politi, S. Ruffo, A. Vulpiani, J. Stat. Phys. 46 (1986) 147. G. Benettin, G.L. Vecchio, A. Tenenbaum, Phys. Rev. A 22 (1980) 1709; M. Casartelli, E. Diana, L. Galgani, A. Scotti, Phys. Rev. A 13 (1976) 1921; Y. Pomeau, A. Pumir, P. Pelce, J. Stat. Phys. 37 (1984) 39. J.D. Weeks, D. Chandler, H.C. Andersen, J. Chem. Phys. 54 (1971) 5237. D.J. Evans, G.P. Morriss, Statistical Mechanics of Nonequilibrium Liquids, Academic Press, London, 1990. G. Benettin, L. Galgani, J.-M. Strelcyn, Phys. Rev. A 14 (1976) 2338. W.G. Hoover, H.A. Posch, Phys. Lett. A 113 (1985) 82; I. Goldhirsch, P.L. Sulem, S.A. Orszag, Physica D 27 (1987) 311; S. Sarman, D.J. Evans, G.P. Morriss, Phys. Rev. A 45 (1992) 2233. We note that with an accuracy, re, of ~ 1%, we cannot hope to resolve differences between individual Lyapunov exponents (except when N < 100). This is because the spacing between the exponents of O(1/N). For N = 100000 the relative separation between the largest and the second largest exponent is obviously ~ 0.001%. To carry out a simulation to this accuracy would require formidable resources since in this case ts ~ 107. D.J. Searles, D.J. Evans, Aust J. Phys. 49 (1996) 39.