In calculations of transient time correlation functions it is convenient to generate the initial ensemble of starting states from a single field free, Gaussian isokinetic trajectory. As Gaussian isokinetic dynamics ergodically generates the isokinetic ensemble, a single field free trajectory is sufficient to sample the ensemble. At equally spaced intervals along this single field free trajectory (every N_e timesteps), field dependent simulations are started and followed for N_n timesteps. The number N_n should be greater than the characteristic time required for the system to relax to a steady state and N_e should be large enough to ensure that the initial phases are uncorrelated. Each of these cycles gives one initial phase Γ, for the transient correlation function. This process can be made more efficient if we use this single equilibrium starting state to provide more than one initial phase for the nonequilibrium trajectories. To do this we use a group of phase space mappings.

In this section we develop mappings of the phase, Γ which have useful properties, for the theoretical interpretation and practical implementation of nonlinear response theory. For convenience we shall write the phase, Γ, as (q,p)=(x,y,z,p_x,p_y,p_z) where each of the components x,y,z,p_x,p_y,p_z is itself an N-dimensional vector. The time evolution of an arbitrary phase variable B(Γ) is governed by the phase variable propagator exp[iLt], so that B(t)=B(Γ(t))=exp[iLt]B(Γ). Note that the propagator is an operator which acts on the initial phases Γ, so in order to calculate the action of the propagator on a phase variable at a time other than zero, B(t) has to be expressed as a function of the initial phases Γ and not the current phases Γ(t). We assume that the equations of motion have no explicit time dependence (by way of a time dependent external field). The propagator is therefore a shift operator. In the time dependent case, the propagator is not a simple shift operator and the results which follow will need to be generalised. We leave this generalisation until Chapter 8.

The phase variable B at time t, B(t) can be traced back to time zero by applying the negative-time phase variable propagator exp[-iLt],

(7.40)

Reversing the sign of the time in the propagator retraces the original trajectory. It is possible to return to the original phase point Γ(0) without changing the sign of the time. This is achieved by mapping the phase point Γ(t) so that a combination of positive time evolution and mapping takes Γ(t)⇒Γ(0). This mapping is called the time reversal mapping M^T. For field free equations of motion, this is straightforward as the mapping simply consists of reversing the signs of all the momenta.

(7.41)

It is important to realise that this process does not lead to a retracing of the original trajectory, as everywhere along the return path the momenta are the opposite sign to those of the forward path. Noting that e^iLt=e^iL(Γ)t, this can be summarised: M^Te^iLt M^Te^iLtΓ(0)=M^Te^iLt M^TΓ(t)=e^-iLtΓ(t)=Γ(0). These results will be derived in more detail later.

Given an initial starting phase Γ=(x,y,z,p_x,p_y,p_z) then four starting phases, which occur within the equilibrium distribution with the same probability as Γ, can be obtained using the mappings M^I, M^T, M^Y and M^K ;

(7.42)

Here M^I is the identity mapping; M^T is the time reversal mapping introduced above; M^Y is termed the y-reflection mapping; and M^K is called the Kawasaki mapping (it is the combined effect of time reversal and y-reflection mapping M^K=M^TM^Y). For shear flow these four configurations give four different starting states, and lead to four different field dependent trajectories from the single equilibrium phase point Γ. Each of the mappings consists of a pair of reflections in a coordinate or momentum axis. In total there are 2⁴ states that can be obtained using the reflections of a 2-dimensional phase space however, only 2³ of these states will result in at most a sign change in the instantaneous shear stress P_xy(Γ). Only 2² of the remaining mappings lead to different shearing trajectories. The shear stress obtained from trajectories starting from Γ₁ and -Γ₁ for example, are identical. The probability of each of these states occurring within the equilibrium distribution, is identical because the Hamiltonian H₀ is invariant under these mappings.

There is a second, more important, advantage of this procedure. If we examine the transient response formula (7.32), we see that at long time the correlation function ⟨B(t)P_xy(0)⟩ approaches ⟨B(∞)⟩⟨P_xy(0)⟩ . The steady state average of B is usually non-zero (in contrast to equilibrium time correlation functions). To minimise the statistical uncertainties in calculating the transient correlation integral, it is necessary to choose equilibrium starting states Γ in such a way that ⟨P_xy(0)⟩≡0. The phase mapping procedure described above achieves this. If the shear stress computed from the original starting phase is P_xy, then the shear stress from Γ^T is also equal to P_xy, but the shear stresses from both Γ^Y and Γ^K are equal to -P_xy. This means that the sum of the shear stresses from these four starting phases is exactly zero, so if each chosen Γ is mapped in this way the average shear stress is exactly zero regardless of the number of samplings of Γ. The statistical difficulties at long time, associated with a small non-zero value the average of P_xy(0), are eliminated.

There are two further operations on these mappings which we need to complete the development of the mapping algebra. First we need to know how each of the mappings affect phase variables. Second we must understand the effect of the mapping on the phase variable Liouvillean iL(Γ), as it is also a function of the initial phase point Γ. To do this we need to know how the equations of motion transform. First we will discuss the transformation of Hamiltonian equations of motion under the mapping, and then consider the transformation of the field dependent dynamics. This will require an extension of the mappings to include the field itself.

To illustrate the development we will consider the time reversal mapping M^T in detail, and simply state the results for other mappings. In what follows the mapping operator M^T operates on all functions and operators (which depend upon Γ and γ) to its right. A particular example is useful at this stage, so consider the shear stress P_xy

(7.43)

Here P_xy is mapped to the same value. For thermodynamically interesting phase variables the operation of the mappings involve simple parity changes

(7.44)

where p^X_B=±1. In the following table we list the values of the parity operators for shear stress, pressure and energy for each of the mappings.

The operation of the mapping M^T on the Hamiltonian equations of motion is

(7.45)

(7.46)

where the transformed coordinate and momenta are denoted by the superscript (T). The vector character of the force F is determined by the coordinate vector q, so that under this mapping the force is invariant. Because and p change sign under the mapping M^T, the phase variable Liouvillean becomes

(7.47)

It is straightforward to use this result and the series expansion of the propagator to show that

(7.48)

To see exactly how this combination of the M^T mapping,and forward time propagation combine to give time reversal we consider the time evolution of Γ itself,

(7.49)

This implies that

. (7.50)

If we start with Γ(0), propagate forward to time t, map with M^T (changing the signs of the momenta), propagate forward to time t, and map with M^T (changing the signs of the momenta again), we return to Γ(0). An analogous identity can be constructed by considering Γ(0)=exp[iL(Γ)t]Γ(-t), that is

(7.51)

This says that we can complete a similar cycle using the backward time propagator exp[iLt] first. These to results demonstrate the various uses of this time reversal mapping.

When the equations of motion for the system involve an external field the time reversal mapping can be generalised to include the field. This is necessary if we wish to determine whether a particular mapping leads to different field dependent dynamics. Here we limit consideration to the isothermal SLLOD algorithm for shear flow. It is clear that all the momenta must change sign so a suitable choice for the mapping is

(7.52)

As the field has units of inverse time the field changes sign together with the momenta. The equations of motion for the mapped variables become

(7.53)

and

(7.54)

Notice also that for the thermostatting variable α

(7.55)

as the numerator changes sign and the denominator is invariant under the time reversal mapping. The mapping of the Liouvillean is similar to the field free case and it can be shown that

(7.56)

In the field dependent case the two operators, equations (7.50, 7.51) generalise to

(7.57)

(7.58)

As a phase variable by definition is not a function of the field, the parity operators associated with mapping phase variables are unchanged.

The second mapping we consider is the y-reflection mapping M^Y, as it acts to change the sign of the shear rate but not the time or the Liouvillean. This mapping is defined by

(7.59)

This mapping consists of a coordinate reflection in the x,z-plane, and momenta reflection in the p_x,p_z-plane. Substituting this mapping into the SLLOD equations of motion shows that the time derivatives of both y and p_y change sign, while the thermostatting variable remains unchanged. The y-reflection Liouvillean is related to the standard Liouvillean by

(7.60)

We now define the combination Kawasaki mapping M^{K ,} which consists of the time reversal mapping followed by the y reflection mapping, so that

(7.61)

Under the Kawasaki mapping the Liouvillean is transformed to

(7.62)

Using the results obtained in this section it easy to show that the following four time evolutions of the phase variable B yield identical values. That is

(7.63)

Notice that these four time evolutions involve changing the sign of the time and/or the sign of the field. If we consider the phase variable P_xy(Γ), the time evolution leads to a negative average value at long time, and where a single sign change is made in the propagator, the parity operator is -1. The third equality has been used to interpret the propagation of the dissipative flux in the Kawasaki exponent; negative time evolution with a positive external field from Γ, is equivalent to positive time evolution with a positive field from Γ^K. As each of the time evolutions in equation (7.63) represent different mathematical forms for the same trajectory, the stabilities are also the same.

The Kawasaki mapping is useful as an aid to understanding the formal expression for the Kawasaki distribution function. The particular term we consider is the time integral of the dissipative flux in the Kawasaki exponent

using the Kawasaki mapping the negative time evolution can be transformed to an equivalent positive time evolution. To do this consider

(7.64)

The last equality follows from the fact that γ^K=γ. So we may think of P_xy(-s,γ,Γ) as equivalent (apart from the sign of the parity operator) to the propagation of P_xy forward in time, with the same γ, but starting from a different phase point Γ^K. The probability of this new phase point Γ^K in the canonical (or isothermal) distributionis the same as the original Γ, as the equilibrium Hamiltonian H₀, is invariant under time reversal and reflection. Therefore the Kawasaki distribution function can be written as

(7.65)

In this form the sign of the exponent itself changes as well as the sign of the time evolution. At sufficiently large time P_xy(s,Γ^K,γ) approaches the steady state value ⟨P_xy(s,γ)⟩, regardless of the initial phase point Γ^K.