Having explored some of the fundamentals of nonlinear response theory, we are now in a better position to comment on one of the early criticisms of linear response theory. In an oft-cited paper van Kampen (1971), criticised linear response theory on the basis that microscopic linearity which is assumed in linear response theory, is quite different from the macroscopic linearity manifest in linear constitutive relations. Van Kampen correctly noted that to observe linear microscopic response (ie of individual particle trajectories) over macroscopic time (seconds, minutes or even hours), requires external fields which are orders of magnitude smaller than those for which linear macroscopic behaviour is actually observed. Therefore, so the argument goes, the theoretical justification of, the Green-Kubo relations for linear transport coefficients, is suspect.

In order to explain his assertion that linearity of microscopic motion is entirely different from macroscopic linearity, van Kampen considered a system composed of electrons which move, apart from occasional collisions with impurities, freely through a conductor. An imposed external electric field, F_e, accelerates the particles between collisions. The distance an electron moves in a time t, under the influence of the field, is 1/2t²(eF_e∕m). In order for the induced current to be linear one requires that t²(eF_e∕2m) << d, the mean spacing of the impurities. Taking d~ 100Å and t to be a macroscopic time, say 1 second, we see that the field must be less than ~10^-18Volts/cm!

As a criticism of the derivation of linear response theory, this calculation implies that for linear response theory to be valid, trajectories must be subject to a linear perturbation over macroscopic times - the time taken for experimentalists to make sensible measurements of the conductivity. This however, is incorrect.

The linear response theory expression for the conductivity, σ (≡J∕F_e) is,

(7.76)

Now it happens that in three dimensional systems the integral of the equilibrium current autocorrelation function converges rapidly. (In two dimensional systems this is expected not to be so.) The integral in fact converges in microscopic time, a few collision times in the above example. Indeed if this were not so one could never use equilibrium molecular dynamics to compute transport coefficients from the Green-Kubo formulae. Molecular dynamics is based on the assumption that transport coefficients for simple fluids can be computed from simulations which only follow the evolution systems for ~10^-10 seconds. These times are sufficient to ensure convergence of the Green-Kubo correlation functions for all the Navier-Stokes transport coefficients. If we require microscopic linearity over 10^-10 seconds (rather than van Kampen's 1 second) then we see that the microscopic response will be linear for fields less than about 100Volts/cm, not an unreasonable number. It simply does not matter that for times longer than those characterising the relaxation of the relevant GK correlation function, the motion is perturbed in a nonlinear fashion. In order for linear response theory to yield correct results for linear transport coefficients, linearity is only required for times characteristic of the decay of the relevant correlation functions. These times are microscopic.

We used nonequilibrium molecular dynamics simulation of shear flow in an atomic system to explore the matter in more detail (Morriss et. al., 1989). We performed a series of simulations with and without an imposed strain rate, γ (≡∂u_x∕∂y), to measure the actual separation d, of phase space trajectories as a function of the imposed strain rate. The phase space separation is defined to be,

(7.77)

where Γ ≡ (q ₁,q ₂,... q _N,p ₁,p ₂, .. ,p _N) is the 6N-dimensional phase space position for the system. In measuring the separation of phase space trajectories we imposed the initial condition that at time zero the equilibrium and nonequilibrium trajectories start from exactly the same point in phase space, d(0,γ)=0, ∀γ. We used the 'infinite checker board' convention for defining the Cartesian coordinates of a particle in a periodic system. This eliminates trivial discontinuities in these coordinates. We also reported the ensemble average of the phase space separation, averaged over an equilibrium ensemble of initial phases, Γ(0,0).

The equations of motion employed were the SLLOD equations. As we have seen the linear response computed from these equations is given precisely, by the Green-Kubo expression for the shear viscosity. The system studied in these simulations was the Lennard-Jones fluid at its triple point (ρ^*=ρσ³=0.8442,T^*=k_BT∕ε=0.722, t^*=t(ε∕m)^1/2σ^-1). A Lees-Edwards periodic system of 256 particles with a potential truncated at, r^*=r∕σ=2.5, was employed.

Before we begin to analyse the phase separation data we need to review some of the relevant features of Lennard-Jones triple point rheology. Firstly, as we have seen (§6.3) this fluid is shear thinning. The strain rate dependent shear viscosities of the Lennard-Jones triple point fluid are set out in the table below.

‡ NEMD estimated.

The most important relevant fact that should be noted from these results is that for reduced strain rates, γ ^*<~10^-2, the fluid is effectively Newtonian with a viscosity which varies at most, by less than ~4% of its zero shear value. (Because of the uncertainty surrounding the zero shear viscosity, we cannot be more certain about the degree of nonlinearity present at γ ^*=0.01.)

The second relevant fact that we should remember is that the GK equilibrium time correlation function whose integral gives the shear viscosity, has decayed to less than 1% of its zero time value at a reduced time t^*=2.0. Values are shown below.

Of course the viscosity which is the time integral of this correlation function converges relatively slowly due to the presence of the slowly decaying t^-3/2 long time tail. Here again there is some uncertainty. If one believes that enhanced long time tail phenomena (§6.3), are truly asymptotic and persist indefinitely then one finds that the viscosity converges to within ~13% of its asymptotic value at t^*=1.0 and to within ~5% of the asymptotic value at t^*=10.0. (If we map our simulation onto the standard Lennard-Jones representation of argon, t^*=1.0 corresponds to a time of 21.6 picoseconds.) If enhanced long time tails are not asymptotic then the GK integrand for the shear viscosity converges to within ~5% of its infinite time value by t^*=2.

The only important observation that concerns us here is that the GK estimate for the shear viscosity is determined in microscopic time, a few hundreds of picoseconds at the very most, for argon. This observation was omitted from van Kampen's argument. We call the range of times required to ensure say 5%, convergence of the GK expression for the viscosity, the GK time window.

Figure 7.12 shows the common logarithm of the ensemble average of the phase space separation plotted as a function of reduced time for various values of the imposed shear rate. The shear rates employed were: γ*=1.0, 10^-1, 10^-2, 10^-3, 10^-5, 10^-7. Note that for the standard Lennard-Jones argon representation, these strain rates correspond to shear rates of 4.6*10¹¹ to 4.6*10⁵ Hz. It will be clear from the present results that no new phenomena would be revealed at strain rates less than γ*~10^-4.

One can see from the figure that at a shear rate of 10^-7, the phase space separation increases very rapidly initially and then slows to an exponential increase with time. The same pattern is followed at a strain rate of 10^-5 except that the initial rise is even more rapid than for a strain rate of 10^-7. Remember that at t=0 the phase space separations start from zero, and therefore the logarithm of the t=0 separations is -∞, for all strain rates.

For strain rates >10^-5, we notice that at long times the phase separation is a constant independent of time. We see an extremely rapid initial rise, followed by an exponential increase with a slope which is independent of strain rate, followed at later times by a plateau. The plateau is easily understood.

The simulations shown in Figures 7.12,13 are carried out at constant peculiar kinetic energy Σp_i ²∕2m = 3Nk_BT. The 3N components of the phase space momenta therefore lie on the surface of a 3N-dimensional sphere of radius, r_T = √(3Nmk_BT). Once the phase space separation exceeds this radius, the curved nature of accessible momentum space will be apparent in our phase space separation plots. The arrow marked on Figure 7.12 shows when the logarithm of the separation is equal to this radius. The maximum separation of phase points within the momentum sub-space is of course 2r. It is clear therefore that the exponential separation must end at approximately, d(t,Γ) = r_T. This is exactly what is observed in Figure 7.12.

Between the plateau and the initial (almost vertical) rise is an exponential region. As can be see from the graph the slope of this rise is virtually independent of strain rate. The slope is related to the largest positive Lyapunov exponent for the system at equilibrium. The Lyapunov exponent measures the rate of separation of phase trajectories that start a small distance apart, but which are governed by identical dynamics. After initially being separated by the external field, the rate of phase space separation thereafter is governed by the usual Lyapunov instability. The fact that the two trajectories employ slightly different dynamics is a second order consideration. The Lyapunov exponents are known to be insensitive to the magnitude of the perturbing external field for field strengths less than 10^-2.

This conjecture regarding the role played by the Lyapunov exponent in the separation of equilibrium and nonequilibrium trajectories which start from a common phase origin is easily verified numerically. Instead of measuring the separation, d, induced by the strain rate, we ran a simulation in which two trajectories started at slightly different phases and which evolved under (identical) zero strain rate equations of motion. The resulting displacement is shown in Figure 7.12 and labelled as ‘lyap’ in the legend. One can see that the slope of this Lyapunov curve is essentially identical to the exponential portions of the strain rate driven curves. The time constants for the exponential portions of the curves are given in Table 7.5.

At this stage we see that even at our smallest strain rate, the trajectory separation is exponential in time. It may be thought that this exponential separation in time supports van Kampen's objection to linear response theory. Surely exponentially diverging trajectories imply nonlinearity? The assertion turns out to be false.

In Figure 7.13 we examine the field dependence of the phase separations in more detail. In this figure we plot the ratio of the separations to the separation observed for a field, γ*=10^-7.

If the ensemble averaged trajectory response is linear then each of the curves in Figure 7.13 will be equispaced horizontal lines. The curves denoted 'av' refer to the ensemble averaged separations shown in Figure 7.12. One can see immediately that within the GK time window, t^* < 2.0, all the separations are linear in the field except for the largest two strain rates γ*=1.0,0.1. We should expect that all strain rates exhibiting linearity within the GK time window should correspond to those systems exhibiting macroscopic linear behaviour (ie. those which are Newtonian). Those exhibiting microscopic nonlinearity within the GK time window should display non-Newtonian macroscopic behaviour. Comparing table 7.3 with Figure 7.12, this is exactly what is seen.

Although systems at a shear rate γ*=10^-2 & 10^-4, do exhibit a nonlinear growth in the phase space separation, it occurs at times which are so late, that it cannot possibly effect the numerical values of the shear viscosity. These nonlinearities occur outside the GK time window.

A possible objection to these conclusions might be: since we are computing ensemble averages of the phase space separations, it might be the averaging process which ensures the observed microscopic linearity. Individual trajectories might still be perturbed nonlinearly with respect to the strain rate. This however, is not the case. In Figure 7.13 the symbols plotted represent the phase space separation induced in single trajectories. For all strain rates a common phase origin is used. We did not average over the time zero phase origins of the systems.

What we see is a slightly noisier version of the ensemble averaged results. Detailed analysis of the un-averaged results reveals that:

The exponential Lyapunov separation is of course only expected to be observed 'on average' either by employing time or ensemble averages. The main point we make here is that even for individual trajectories where phase separation is not exactly exponential in time, trajectory separation is to 4 significant figure accuracy, linear in the field. The linearity of the response is not produced by ensemble averaging.

We conclude from these studies that within the GK time window, macroscopic and microscopic linearity are observed for identical ranges of strain rates. For times shorter than those required for convergence of the linear response theory expressions for transport coefficients, the individual phase space trajectories are perturbed linearly with respect to the strain rate for those values of the strain rate for which the fluid exhibits linear macroscopic behaviour. This is in spite of the fact that within this domain the strain rate induces an exponential separation of trajectories with respect to time. We believe that many people have assumed an exponential trajectory separation in time implies an exponential separation with respect to the magnitude of the external field. This work shows that within the GK time window, the dominant microscopic behaviour in fluids which exhibit linear macroscopic behaviour, is linear in the external field but exponential in time.

We have seen in Figure 7.12 that for intermediate times the phase separation takes the form,

(7.78)

where the Lyapunov time, τ_L, is the inverse of the largest Lyapunov exponent for the system at equilibrium. We can explain why the phase separation exhibits this functional form and moreover, we can make a rough calculation of the absolute magnitude of the coefficient, A. We know that the exponential separation of trajectories only begins after a time which is roughly the Maxwell relaxation time τ_M, for the fluid. Before the particles sense their mutual interactions, the particles are freely streaming with trajectories determined by the initial values of . After this initial motion the particles will have coordinates and momenta as follows,

(7.79)

When this approximation breaks down, approximately at the Maxwell relaxation time, τ_M≡η∕G, the phase separation d(τ_M,γ) will be,

(7.80)

For our system this distance is,

(7.81)

We have used the fact that the reduced Maxwell time is 0.137. After this time the phase separation can be expected to grow as,

(7.82)

where, as before τ_L is the inverse of the largest zero-strain rate Lyapunov exponent. For fields less than γ ^*=10^-2, the equilibrium Lyapunov time dominates the denominator of the above expression. This explains why the slopes of the curves in Figure 7.12 are independent of strain rate. Furthermore by combining equations (7.70), (7.81) and (7.82) we see that in the regime where the strain rate corrections to the Lyapunov exponents are small, the phase separation takes the form given by equation (7.78) with the coefficient, A~8.7. Equation (7.82) is plotted, for a reduced strain rate of 10^-7, as a dashed line in Figure 7.12. It is in reasonable agreement with the results. The results for other strain rates are similar. The greatest uncertainty in the prediction is the estimation of the precise time at with Lyapunov behaviour begins.