In this chapter we will discuss how an external field Fe, perturbs an N-particle system. We assume that the field is sufficiently weak that only the linear response of the system need be considered. These considerations will lead us to equilibrium fluctuation expressions for mechanical transport coefficients such as the electrical conductivity. These expressions are formally identical to the Green-Kubo formulae that were derived in the last chapter. The difference is that the Green-Kubo formulae pertain to thermal transport processes where boundary conditions perturb the system away from equilibrium - all Navier-Stokes processes fall into this category. Mechanical transport coefficients on the other hand, refer to systems where mechanical fields which appear explicitly in the equations of motion for the system, drive the system away from equilibrium.
As we will see it is no coincidence that there is such a close similarity between the fluctuation expressions for thermal and mechanical transport coefficients. In fact one can often mathematically transform the nonequilibrium boundary conditions for a thermal transport process into a mechanical field. The two representations of the system are then said to be congruent.
A major difference between the derivations of the equilibrium fluctuation expressions for the two representations is that in the mechanical case one does not need to invoke Onsager's regression hypothesis. The linear mechanical response of a nonequilibrium system is analysed mathematically with resultant expressions for the response that involve equilibrium time correlation functions. In the thermal case - Chapter 4 - equilibrium fluctuations were studied and after invoking Onsager's hypothesis, the connection with nonequilibrium transport coefficients was made. Given a congruent mechanical representation of a thermal transport process, one can in fact prove the validity of Onsager's hypothesis.
The mechanical field Fe, performs work on the system, preventing relaxation to equilibrium. This work is converted into heat. It is easy to show that the rate at which the field performs work on the system is, for small fields, proportional to F2e. As such this is, at least formally, a nonlinear effect. This is why, in the complete absence of any thermostatting mechanism, Kubo (1957) was able to derive correct expressions for the linear response. However in spite of heating being a nonlinear effect, a thermostatted treatment of linear response theory leads to a considerably more satisfying discussion. We will therefore include in this chapter a description of thermostats and isothermal linear response theory.
Consider a system of N atoms suddenly subject, at t=0, to a time dependent external field, Fe. The generalisation of our discussion to vector or tensor fields is straightforward. For simplicity we will assume that the particles move in a three dimensional Cartesian space. For times greater than zero the system is assumed to obey the dynamics given in the equations below,
(5.1)
The phase variables Ci(Γ) and Di(Γ) describe the coupling of the field to the system. We assume that the equations have been written in such a way that at equilibrium in the absence of the external field the canonical kinetic energy K, satisfies the equipartition relation,
(5.2)
This implies that the canonical momenta give the peculiar velocities of each of the particles and that therefore,
(5.3)
In this case H0,
(5.4)
is the instantaneous expression for the internal energy. We do not assume that a Hamiltonian exists which will generate the field-dependent equations of motion. In the absence of the external field and the thermostat, H0 is the total energy, and is therefore a constant of the motion. The rate of change of internal energy due to the field is
(5.5)
where J(Γ), is called the dissipative flux.
The response of the system to the external field can be assessed by monitoring the average response of an arbitrary phase variable B(Γ) at some later time t. The average response of the system is the response that is obtained by perturbing an ensemble of initial phases. It is usual to select the starting states from the equilibrium canonical ensemble, thus
(5.6)
The average response ⟨B(t)⟩ can be calculated from the expression,
(5.7)
This is the Schrödinger representation for the response of B. The problem of determining the response then reduces to determining the perturbed distribution function f(t). The rate of change in the perturbed distribution function is given by the Liouville equation
(5.8)
The Γ(t) in these equations is given by the first order form of the equations of motion with the external field evaluated at the current time, t.
If the equations of motion are derivable from a Hamiltonian it is easy to show that 
, (§3.3). We will assume that even in the case where no Hamiltonian
exists which can generate the equations of motion (5.1), that 
. We refer to this condition as the Adiabatic Incompressibility of
Phase Space (AIΓ). A sufficient, but not necessary, condition for this to hold is
that the unthermostatted or adiabatic equations of motion are derivable from a
Hamiltonian. It is of course possible to pursue the theory without this condition but in
practise it is rarely necessary to do so (the only known exception is discussed: Evans
and MacGowan, 1986).
Thus in the adiabatic case if AIΓ holds, we know that the Liouville operator is Hermitian (see §3.3 & §3.5) and therefore,
(5.9)
If we denote the Liouvillean for the field free equations of motion as iL0, and we break up the total Liouvillean into its field free and field dependent parts, equation (5.8) becomes,
(5.10)
where the distribution function f(Γ,t), is written as fc+∆f(Γ,t). Since H0 is a constant of the motion for the field free adiabatic equations of motion, iL0 therefore preserves the canonical ensemble,
(5.11)
Substituting (5.11) into equation (5.10) we see,
(5.12)
In (5.12) we are ignoring perturbations to the distribution function which are second order in the field. (The Schrödinger-Heisenberg equivalence (§3.3), proves that these second order terms for the distribution are identical to the second order trajectory perturbations.) In §7.8 we discuss the nature of this linearisation procedure in some detail. To linear order, the solution of equation (5.12) is,
(5.13)
The correctness of this solution can easily be checked by noting that at t=0, (5.13) has the correct initial condition, (∆f(Γ,t=0)=0) and that the solution for ∆f(Γ,t) given in (5.13) satisfies (5.12) for all subsequent times.
We will now operate on the canonical distribution function with the operator, iL(t). We again use the fact that iL0 preserves the canonical distribution function.
(5.14)
The adiabatic time derivative of H0 is given by the dissipative flux (5.5), so,
(5.15)
The time argument associated with i∆L(s) is the time argument of the external field.
Substituting (5.15) into (5.13) and in turn into (5.7), the linear response of the phase variable B is given by
(5.16)
In deriving the third line of this equation from the second we have unrolled the propagator from the dissipative flux onto the response variable B. Note that the propagator has no effect on either the canonical distribution function (which is preserved by it), or on the external field Fe(t) which is not a phase variable.
It is usual to express the result in terms of a linear susceptibility χBJ, which is defined in terms of the equilibrium time correlation function of B and J,
(5.17)
To linear order, the canonical ensemble averaged linear response for B(t) is,
(5.18)
This equation is very similar to the response functions we met in Chapter 4 when we discussed viscoelasticity and constitutive relations for thermal transport coefficients. The equation shows that the linear response is non-Markovian. All systems have memory. All N-body systems remember the field history over the decay time of the relevant time correlation function, ⟨B(t)J(0)⟩. Markovian behaviour is only an idealisation brought about a lack of sensitivity in our measurements of the time resolved many-body response.
There are, a number of deficiencies in the derivation we have just given. Suppose that by monitoring ⟨B(t)⟩ for a family of external fields Fe, we wish to deduce the susceptibility χ(t). One cannot blindly use equation (5.18). This is because as the system heats up through the conversion of work into heat, the system temperature will change in time. This effect is quadratic with respect to the magnitude of the external field. If χ increases with temperature, the long time limiting value of ⟨B(t)⟩ will be infinite. If χ decreases with increasing temperature the limiting value of ⟨B(t)⟩ could well be zero. This is simply a reflection of the fact that in the absence of a thermostat there is no steady state. The linear steady state value for the response can only be obtained if we take the field strength to zero before we let time go to infinity. This procedure will inevitably lead to difficulties in both the experimental and numerical determination of the linear susceptibilities.
Another difficulty with the derivation is that if adiabatic linear response theory is applied to computer simulation, one would prefer not to use canonical averaging. This is because a single Newtonian equilibrium trajectory cannot generate or span the canonical ensemble. A single Newtonian trajectory can at most span a microcanonical subset of the canonical ensemble of states. A canonical evaluation of the susceptibility therefore requires an ensemble of trajectories if one is using Newtonian dynamics. This is inconvenient and very expensive in terms of computer time.
One cannot simply extend this adiabatic theory to the microcanonical ensemble. Kubo (1982) recently showed that if one subjects a cumulative microcanonical ensemble (all states less than a specified energy have the same probability) to a mechanical perturbation, then the linear susceptibility is given by the equilibrium correlation of the test variable B and the dissipative flux J, averaged over the delta microcanonical ensemble (all states with a precisely specified energy have the same probability). When the equilibrium ensemble of starting states is not identical to the equilibrium ensemble used to compute the susceptibilities, we say that the theory is ergodically inconsistent. We will now show how both of these difficulties can be resolved.