The first evidence for the existence of a strange attractor in the phase space of the two-dimensional, two-body planar Couette flow system was obtained by Morriss (1987). He showed numerically that the information dimension of two-body planar Couette flow is a decreasing function of the strain rate, dropping steadily from three towards a value near two, before dropping dramatically at a critical value of the strain rate to become asymptotic to one. These results are for the WCA potential (equation 6.5) at a reduced temperature of 1 and a reduced density of ρ=0.4. The sudden change in dimension, from a little greater than two to near one, is associated with the onset of the string-phase for this system (see §6.4). A change in effective dimensionality for shearing systems of 896 particles, under large shear rates, has been observed. In this case the vector separation between two atoms r_ij=(x_ij,y_ij) has components whose sign is independent of time. This arises because within strings the atoms are ordered, and the strings themselves once formed remain forever intact, (and in the same order). It has been shown that the string phase is an artifact of the definition of the temperature with respect to an assumed streaming velocity profile (§6.4), so it is likely that this decrease in dimensionality is pathological, and not associated with the attractor which is found at intermediate strain rates.

Morriss (1989) has calculated the generalized dimension D_q and the spectrum of singularities f(α) for the steady state phase space distribution function of two dimensional two-body planar Couette flow using the WCA potential at a reduced temperature of 1 and a reduced density of 0.4. This system is identical to that considered in the information dimension calculations referred to above. The maximum resolution of the distribution function was 3x2⁶ bins in each of the three degrees of freedom, leading to more accurate results than the previous information dimension calculations. He found that at equilibrium the discrete probabilities p_i(ε) scale with the dimension of the initial phase space. Away from equilibrium the p_i(ε) scale with a range of indices, extending from the full accessible phase space dimension to a lower limit which is controlled by the value of the shear rate γ.

In Figure 10.8 we present the singularity distribution f(α) for a number of values of the strain rate γ. The results near γ=0 depend significantly on the values of grid size used, and could be improved by considering finer meshes (the minimum grid size is limited by computer memory size). At higher values of γ (say γ=1) the values of f(α) above the shoulder in Figure 10.8, are insensitive to grid size. However, the position of the shoulder does change with grid size. In the limit q→∞, the value of D_q and hence the value of α=α_min for which f(α)→0, is controlled by the scaling of the most probable p_i in the histogram p_max. It is easy to identify p_max and determine its scaling as an independent check on the value of α_min. Just as large positive values of q weight the most probable p_i most strongly, large negative values of q weight the least probable p_i most strongly. The accuracy with which the least probable p_i can be determined limits the minimum value of q for which the calculation of D_q is accurate. This is reflected in poor values of D_q for q≤0.5, and we believe is a contributing factor in obtaining inconsistent values of the fractal dimension D₀.

We interpret the results shown in Figure 10.8 as follows. The value of f(α) is the dimension of the set of points on the attractor which scale as ε^α in the discrete phase space distribution function {p_i}. For this system it implies singularities of the form |Γ-Γ ₀|^α-3 in the underlying (continuous) phase space distribution function f(Γ,γ). At equilibrium most p_i's scale as ε³, with a very narrow spread of lower α values. Indeed with finer grid sizes this distribution may narrow still further. Away from equilibrium two effects are clearly discernible. First the dimension of the set of p_i's which scale as ε³ drops with increasing γ. Second the distribution of values of α increases downwards with the lower limit α_min controlled by the value of γ. This distribution is monotonic with the appearance of a shoulder at an intermediate value of α.

Having calculated the full phase space distribution function on a resolution ε we can investigate the behavior of the various reduced distributions, for example we may consider the coordinate space distribution function f₂(r,Φ), or the distribution of the momentum angle . Each of these reduced distributions is obtained by integrating (or summing) over the redundant coordinates or momenta. Perhaps the most interesting of these reduced distribution functions is the coordinate space distribution f₂(x₁₂,y₁₂), shown in Figure 10.9.

If the underlying continuous distribution function has a singularity of the form |Γ-Γ ₀|^α-3, then f₂ can have singularities of the form |Γ-Γ ₀|^α-2. However, if 2≤α≤3 then these points are no longer singularities, and the reduced distribution f₂ has a different character to the full phase space distribution. If the exponent α-2 is positive, then f₂ is zero at Γ ₀ and the discrete probability p_i(ε) which includes Γ ₀ will scale as ε², whereas if α-2 is negative then f₂ is singular.

In this study all the two variable distribution functions, although being highly structured in many cases, did not show any evidence of singularity sets of non-zero measure. This observation has important ramifications for the Green entropy which we will meet in §10.5.

The complete set of Lyapunov exponents for two and three-dimensional planar Couette flow have been calculated for 2, 4 and 8 particle systems by Morriss (1988,1989). For the two particle system the Lyapunov dimension D_L has been calculated using both the Mori and Kaplan-Yorke conjectures (equations 10.46 and 10.47). This requires the complete set of Lyapunov exponents (that is 3 exponents for N=2) and has the advantage over static dimensions that no subsequent extrapolation procedure is needed. The following table contains the results for the two-body, two-dimensional Couette flow system at the same state point as that used in the information and generalized dimension calculations.

For both the Kaplan-Yorke and Mori forms, the Lyapunov dimension is found to be a decreasing function of the shear rate. This is consistent with the contraction of phase space dimension that we have already seen from the numerical evaluated static dimensions D_q. It confirms that the nonequilibrium distribution function is a fractal attractor whose dimension is less than that of the equilibrium phase space. When the shear rate γ is zero, both methods of calculating the Lyapunov dimension agree. However, as soon as the shear rate changes from zero, differences appear. In the Kaplan-Yorke formula (equation 10.47), the value of n is 2 from γ=0, until the magnitude of λ₂ exceeds that of λ₁ (somewhere between γ=2 and 2.25). This means that 2<D^KY_L<3 in this range. For γ>2, 1<D^KY_L<2 as long as λ₁ remains positive. The value of λ₃ is irrelevant as soon as |λ₂|>λ₁. Then as λ₁ becomes negative the dimension is equal to zero. The Kaplan-Yorke formula can never give fractional values between zero and one. In the Mori formula the value of λ₃ always contributes to the dimension, and its large negative value tends to dominate the denominator, reducing D^M_L. The transition from D^M_L>2 to D^M_L>2 is somewhere between γ=1 and 1.25. Indeed the Mori dimension is systematically less than the Kaplan-Yorke dimension.

In Table 10.2 we compare the values of D_q for 2 particle two-dimensional planar Couette flow for several values of q, with the Kaplan-Yorke Lyapunov dimension for this system obtained from the full spectrum of Lyapunov exponents. Of the two routes to the Lyapunov dimension the Kaplan-Yorke method agrees best with the information dimension results of Table 10.2, whereas the Mori method does not. In particular the Kaplan-Yorke method and the information dimension both give a change from values greater than two, to values less than two at about γ=2.25. There are a number of points to note about the results in this table. First, it can be shown that D₁ is a lower bound for D₀, however the numerical results for D₀ and D₁ are inconsistent with this requirement as D₀<D₁. As we remarked previously, the results for D_q when q<0.5 are poor. It has been argued that the fractal (Hausdorff) dimension and Kaplan-Yorke Lyapunov dimension should yield the same result, at least for homogeneous attractors. In this work we find that D^KY_L is significantly lower than D₁ (which is itself a lower bound on D₀) for all values of the strain rate. Indeed D^KY_L is approximately equal to D_q, where q somewhat greater than 3.

It is possible to calculate the Lyapunov exponents of systems with more than two particles, whereas extending the distribution function histograming algorithms for the information dimension or generalized dimension is much more difficult. The full Lyapunov spectrum has been calculated for 4 and 8 particle planar Couette flow systems in both two and three dimensions.

In Figure 10.10 we show the Lyapunov spectra for the 4 particle system at ρ=0.4 for a range of values of the shear rate. For the equilibrium spectrum (γ=0) one exponent is zero, while the others occur in conjugate pairs {λ_-i,λ_i}, where λ_-i=λ_i. This symmetry is a consequence of the time reversibility of the equations of motion and the conservation of phase space volume from the Liouville theorem. For the two-dimensional system the exponents appear to be essentially linear in exponent number, but a linear fit to the positive exponents is not consistent with an exponent of zero for exponent number zero. As the external field is increased systematic changes in the Lyapunov spectrum occur.

The positive branch decreases, with the smallest positive exponent decreasing most. The largest positive exponent seems almost independent of the external field. We expect that the most vigorous mixing in phase space, which is controlled by the positive exponents, is first a function of the curvature of the particles themselves (the higher the curvature, the more defocusing is each collision), and second depends on the collision frequency (and hence the density). It could be argued that the insensitivity of the largest exponent is associated with only a small change in collision frequency with strain rate, at this density. The zero exponent becomes more negative with increasing field, as does the negative branch of the Lyapunov spectrum. The change in the negative branch is larger than the change in the positive branch. The change in the sum of each exponent pair is the same, that is λ_-i+λ_i=c, where c is constant independent of i and related directly to the dissipation. The change in the exponent which is zero at equilibrium is ½c.

The idea of being able to characterize the Lyapunov spectrum without having to calculate all of the exponents is very attractive, as the computation time for the Gaussian constraint method depends on the fourth power of the number of particles N. We decided to compare the Lyapunov spectra as a function of system size, at the same state point. It is well known that the bulk properties will have some system size dependence, but the trends as a function of density and temperature should be reliable. In Figure 10.11 we present the Lyapunov spectra for an equilibrium system at ρ=0.4, for a range of system sizes N=2,4 and 8. Each spectra is scaled so that the largest positive and negative exponents have the same exponent number regardless of system size. These results look very encouraging as the spectra of all three systems are very similar. The linear fit to the positive branch for N=4 and N=8 have slightly different slopes but the qualitative features are the same.

In Figure 10.12 we present the Lyapunov spectra for a strain rate of γ=1.0 at ρ=0.4, for system sizes of N=2,4 and 8. This shows that there is also a close correspondence between the results at different system sizes away from equilibrium.

In Figure 10.13 we show the Lyapunov dimension of the planar Couette flow system at ρ=0.4 as a function of strain rate, for a range of system sizes. For each system size the Lyapunov dimension is scaled by the equilibrium value, so that the plotted results represent the proportional reduction in dimension. The qualitative trends are the same. There is a decrease in dimension with increasing strain rate. The proportional change in dimension is greatest for the two particle system and smallest for the eight particle system, whereas the absolute changes are in the opposite order.

In summary, the results confirm the dimensional contraction observed previously in two body, two-dimensional planar Couette flow simulations. The initial phase space dimension of D=2dN-2d-1, contracts with increasing external field, and the distribution function is only nonzero on a fractal attractor of dimension less than 2dN-2d-1. Although the results for these systems differ in detail from the generalized dimension results, the observation of significant dimensional contraction is universal. An approach which may help reduce the magnitude of the numerical calculations is the observation that the qualitative features of the spectra are essentially independent of system size.

If we consider the volume element V_2dN where 2dN is the phase space dimension of the initial system (d is the spatial dimension and N is the number of particles), then we have that the phase space compression factor gives the rate of change of phase space volume (see equation 3.78), so that the average of the divergence is equal to the sum of the Lyapunov exponents. A careful calculation of the divergence for the SLLOD algorithm, taking into account the precise number of degrees of freedom gives

(10.51)

where ⟨P^K_xy⟩ is the kinetic contribution to the shear stress and V is the volume. The term involving ⟨P^K_xy⟩ is order one whereas the first term is order N, so for many particle systems the second term can be ignored. For the systems considered here both terms must be included. This is a valuable consistency check on the accuracy of the numerical calculation of Lyapunov exponents.

We have now identified two effects associated with the phase space distribution functions of nonequilibrium systems; the first was dimensional contraction, and the second is a range of sets of fractional power law singularities. The two results are consistent in the sense that as each distribution function is normalized, the loss of probability due to dimensional contraction, is compensated for by the appearance of singularities in the distribution function.

Studies of two and three-dimensional colour diffusion systems by Posch and Hoover (1988) have produced an impressive calculation - the full Lyapunov spectrum for a three dimensional system of 32 repulsive Lennard-Jones atoms (185 Lyapunov exponents) - as well as results for the same system with 8 atoms. Lyapunov spectra are insensitive to ensemble, both at and away from equilibrium. All indications are that nonequilibrium systems are also insensitive to the details of the ensemble or thermostatting mechanism. On the other hand boundary effects do have a significant influence on the shape of spectra for small system. In particular, the homogeneous algorithms for shear flow (such as SLLOD) give different Lyapunov exponents to boundary reservoir methods (Posch and Hoover, 1989).

As small NEMD simulations of planar Couette flow and colour diffusion are dominated by a fractal attractor whose dimension is determined by the strength of the applied field, this behaviour can be expected for all nonequilibrium steady state simulations. The existence of a fractal attractor is a vital clue to understanding the nonequilibrium entropy, but as yet we only have information concerning the rate of approach of a trajectory to the attractor, and measures of its effective dimension. We know a good deal about the structure of the attractor, and the singularities of the nonequilibrium distribution function. Some recent work in the study of dynamical systems (Takahashi and Oono, 1984) shows that modeling chaotic behaviour with statistical mechanical analogues is a useful approach however, but to date the approach parallels irreversible thermodynamics with a continuous production of entropy. For a theory of nonequilibrium steady states, we need to be able to calculate an entropy shift from equilibrium to the steady state which is finite. The appearance of an attractor, and the relative stability of entropy producing trajectories provides a plausible explanation for the observation of irreversibility and a mechanism for the resolution of Löschmidt's paradox (Holian, Hoover and Posch, 1987).

It is interesting to make a connection between the results given here and the numerical calculations based on the Kawasaki distribution function. In §7.7 we described some very recent numerical studies of the Kawasaki form for the full nonlinear response of an equilibrium system subject to the sudden application of a fixed shear rate. From a theoretical point of view there are two points of interest in this stress growth experiment. First, is the renormalized Kawasaki shear stress equal to that observed directly? Second, how does the Kawasaki normalization behave as a function of time? The results show that the renormalized Kawasaki shear stress is in quite good agreement with the direct result, and that the Kawasaki normalization which is one initially, decreases with time. The results obtained here for the 2-body system suggest that the Kawasaki distribution function may have singularities which compensate for the observed decrease in both the individual probabilities and the normalization, and that these singularities are not adequately represented in the phase space sampling used.

Equation (10.4) implies that if we consider a comoving phase space volume element containing a fixed number of trajectories, then the local density of phase space increases indefinitely because the associated Lagrangian volume is constantly decreasing (because the sum of the Lyapunov exponents is negative). Since the contraction of the accessible phase space is continuous there is in a sense, no possibility of generating a steady state distribution function. Computed from the ostensible phase space the volume of accessible phase space shrinks at a constant rate becoming zero at infinite time. A steady state is in fact characterized by a constant rate of decrease in the relative volume occupied by accessible phase space. This is in spite of the fact that in a steady state averages of phase variables are constant. This apparent contradiction can be understood by considering the following example.

Suppose we consider a system which at t=0 occupies a 2-dimensional phase space 0<x,y<L. Suppose that by some means this thermostatted system is subject to a dissipative external field which, after initial transients, causes the distribution function, f(x,y), to collapse towards a one dimensional attractor, x²+y²=r². At some time t, the distribution function is given by the equation,

(10.52)

Further, we suppose that the width of the annulus which forms the distribution function satisfies an equation of motion,

(10.53)

for some positive constant value of α. It is easy to see that in the steady state , df∕dt = αf, which is the analog of (10.4). The phase space distribution function diverges at a constant rate, α. In spite of this, if we compute the phase average of a nonsingular phase variable B(x,y), time averages will converge exponentially fast towards their steady state values, <B(t)> - <B(∞)> ~ e^-αt. This example points out that although the distribution function, as computed from the ostensible phase space, may be diverging at a constant rate, steady state phase averages may still be well defined and convergent. The distribution function computed from within the accessible phase space has no singularities, f_acc(x,y,t) ≡ f(x,y,t)∕(2πrΔ(t)) = 1, ∀t, provided, r² < x²+y² < (r+Δ(t))². In our example it is always uniform and constant in time. Phase averages are fundamentally functions of phase space distances not of volumes. Indeed the notion of a phase space volume is somewhat arbitrary.