Surprisingly often we are interested in the intermediate regime where the Green-Kubo method cannot be applied and where, because of noise, direct NEMD is very inefficient. We have just seen how the TTCF method may be applied to strong fields. It is also the most efficient known method for treating fields of intermediate strength. Before we demonstrate the application of TTCFs to the small field response, we will describe an early method that was used to calculate the intermediate field response.
Prior to the development of Transient Time Correlation Function method, the only way of computing the small but finite field response of many-body systems was to use the Subtraction or Differential response method. The idea behind this method is extremely simple. By considering a sufficiently small field, the systematic response (ie the field induced response) will be swamped by the natural (essentially equilibrium) fluctuations in the system. However it is clear that for short times and small applied fields, there will be a high degree of correlation in the transient response computed with, and without, the external field, (see Fig. 7.6).
Figure 7.6 We depict the systematic nonequilibrium response (the shaded curve) as the difference of the nonequilibrium response from the equilibrium response. By taking this difference we can dramatically reduce the noise in the computed systematic nonequilibrium response. To complete this calculation one averages this differenc over an ensemble of starting states.
If we compute A(t) for two trajectories which start at the same phase, Γ, one with the field on and the other with the field off, we might see what is depicted in Figure 7.6. Ciccotti et. al. (1975, 1976, 1979), realised that, for short times, the noise in A(t) computed for the two trajectories, will be highly correlated. They used this idea to reduce the noise in the response computed at small applied fields.
To use their Subtraction Method one performs an equilibrium simulation (Fe=0), from which one periodically initiates nonequilibrium calculations (Fe≠0). The general idea is shown in Figure 7.7. The phases {Γ i}, are taken as time origins from which one calculates the differences of the response in a phase variable with and without the applied field. The systematic or nonequilibrium response is calculated from the equation,
(7.66)
For many years this was the only method of calculating the small field nonequilibrium response. It suffers from a major problem however. For the method to work, the noise in the the value of A(t) computed with and without the field, must be highly correlated. Otherwise the equilibrium fluctuations will completely swamp the desired response. Now the noise in the two responses will only be correlated if the two systems remain sufficiently close in phase space. The Lyapunov instability (§3.4) will work against this. The Lyapunov instability will try to drive the two systems apart exponentially fast. This can be expected to lead to an exponential growth of noise with respect to time. This is illustrated in Figures 7.8,9 in which the TTCF, denoted (T), and Subtraction techniques, denoted (sub), are compared for the 256 particle WCA system considered in §7.5.
Figure 7.8 shows the shear stress for the three dimensional WCA system at the comparatively small strain rate of γ * = 10-3. At this field strength conventional steady state NEMD is swamped by noise. However the Subtraction technique can be used to substantially improve the statistics. It is important to note that both the Subtraction and TTCF technique are based on an analysis of the transient response of systems. The results compared in Figure 7.8 were computed for exactly the same system using exactly the same data. The only difference between the two sets of results is how the data were analysed. Lyapunov noise is clearly evident in the Subtraction results labelled in Figure 7.8 as Pxy(sub). For longer times, during which we expect the slow nonlinearities to complete the relaxation to the steady state, the Subtraction technique becomes very noisy.
Figure 7.8 Shear stress for the three-dimensional WCA system at a strain rate of γ*=10-3. sub, subtraction technique; T, TTCF.
Figure 7.9 shows the corresponding results for shear dilatancy. Here the Subtraction technique (labelled ‘sub’), is essentially useless. Even the TTCF method becomes somewhat noisy at long times. The TTCF results clearly show the existence of a measurable, intrinsically nonlinear effect even at this small strain rate.
Although the TTCF method allows us to compute the response of systems to fields of arbitrary, even zero, strength, we often require more information about the small field response than it is capable of providing. For example at small fields the response is essentially linear. Nonlinear effects that we may be interested in are completely swamped by the linear response terms. The Differential Transient Time Correlation Function (DTTCF) is an attempt to provide an answer to this problem. It uses a subtraction technique on the TTCFs themselves to formally subtract the linear response.
Figure 7.9 Shear dilatancy for the three-dimensional WCA system at a strain rate of γ*=10-3. sub, subtraction technique; T, TTCF.
In the DTTCF method we consider the difference between B(s) evaluated with and without the external field, starting from the same initial phase point. From the transient correlation function expression this gives
(7.67)
In this equation B(s,γ) is generated from B(0) by the thermostatted field dependent propagator. B(s,0), on the other hand is generated by the zero-field thermostatted propagator. The last term is the integral of an equilibrium time correlation function. This integral is easily recognised as the linear, Green-Kubo response. The first term on the RHS is the integral of a differential transient time correlation function (DTTCF), and is the intrinsically nonlinear response. The LHS is termed the direct differential, or subtraction average.
There are two possible cases; the first in which B has a non-zero linear response term, and the second where the linear response is identically zero. If B is chosen to be Pxy the third term in (7.6.2) is the Green-Kubo expression for the response of the shear stress -η(0)γ, where η(0) is the zero shear rate shear viscosity. The definition of the shear rate dependent viscosity, η(γ)=-⟨Pxy⟩∕γ gives
(7.68)
as the intrinsically nonlinear part of the shear viscosity. As s→∞ the differential transient time correlation function (using the mixing assumption) becomes <(Pxy(s,γ)-Pxy(s,0))> <Pxy> = <Pxy(s,γ)> <Pxy> . This is zero because <Pxy(0)> is zero. On the other hand <(Pxy(s,γ)> is clearly non-zero which means that the use of our trajectory mappings will improve the statistics as s→∞.
To apply the phase space mappings in the differential response method we consider the identity (7.63). We can obtain four different time evolutions of B(Γ) by simply removing the minus signs and parity operators from each of the equivalent forms in equation (7.63). If we use the index α to denote the 4 mappings {I,T,Y,K}, then
(7.69)
This is the direct response of the phase variable B(Γ) from one sampling of Γ, where the mappings are used to generate four starting phase points. To calculate the differential response of B we need to subtract the field free time evolution of B(Γ) from each of these four starting states. The four field free time evolutions are found by setting γ α=0 in equation (7.69). That is
(7.70)
Clearly there are only two different field free time evolutions; the remaining two can be obtained from these by the sign changes of the parity operators. In practice, a single cycle of the numerical evaluation of a differential transient time correlation function will involve the calculation of four field dependent trajectories and two field free trajectories, yielding four starting states.
The use of the symmetry mappings implies some redundancies in the various methods of calculating the response. In particular the direct response of Pxy(t) is exactly equal to the direct differential response for all values of the time. This means that the contribution from the field free time evolutions is exactly equal to zero. This is easy to see from equation (7.69) as there are only two different time evolutions; those corresponding to exp[iLt] and exp[-iLt] respectively, and for Pxy each comes with a positive and negative parity operator. Therefore these two responses exactly cancel for all values of time.
The second redundancy of interest is that the transient response of the pressure p(t) is exactly equal to the differential transient response for all values of time. This implies that the contribution to the equilibrium time correlation function ⟨p(t)Pxy⟩ from a single sampling of Γ is exactly zero. Clearly this equilibrium time correlation is zero when the full ensemble average is taken, but the result we prove here is that the mappings ensure that Σ p(t) Pxy is zero for each starting state Γ for all values of t. The contribution from the field free trajectories is
(7.71)
Again the product of parities ensures that the two field free time evolutions exp[iLt], and exp[-iLt] occur in cancelling pairs. Therefore the field free contribution to the differential time correlation function is exactly zero and the differential transient results are identical the transient correlation function results.
The DTTCF method suffers from the same Lyapunov noise characteristic of all differential or subtraction methods. In spite of this problem Evans and Morriss (1987) were able to show, using the DTTCF method, that the intrinsically nonlinear response of 3-dimensional fluids undergoing shear flow is given by the classical Burnett form (see §9.5). This is at variance with the nonclassical behaviour predicted by mode coupling theory. However, nonclassical behaviour can only be expected in the large system limit. The classical behaviour observed by Morriss and Evans (1987), is presumably the small strain rate, asymptotic response for finite periodic systems.
A much better approach to calculating and analysing the asymptotic nonlinear response will be discussed in §9.5.