For a system undergoing planar Couette flow the transient correlation function expression for the canonical ensemble average of a phase variable B is,

(9.1)

This expression relates the nonequilibrium value of a phase variable B at time t, to the integral of a transient time correlation function (the correlation between P_xy in the equilibrium starting state, P_xy(0), and B at time s after the field is turned on). The temperature implied by the β is the temperature of the initial ensemble. The steady state is tied to the initial ensemble by the constraint of constant peculiar kinetic energy. For systems that exhibit mixing, equation (9.1) can therefore be rewritten as,

(9.2)

where the difference variable ΔB(s) is defined as the difference between the phase variable at s and its average value at s,

(9.3)

Systems which are not expected to exhibit mixing are turbulent systems or systems which execute quasi-periodic oscillations.

An important property of the Gaussian thermostat is that although it fixes the kinetic energy of a system, the Gaussian isokinetic Liouville operator is independent of the temperature of the initial distribution. For each member of the ensemble, the Gaussian thermostat simply constrains the peculiar kinetic energy to be constant. As the Liouvillian, and the propagator in (9.2), are independent of the value of the temperature we can calculate the temperature derivative very easily. The result is,

(9.4)

The first term on the right hand side of (9.4) is the equilibrium contribution. This is easily seen by setting t=0. The second and third terms are nonequilibrium terms. In deriving the second term on the right hand side of (9.4) we use equation (9.3) to simplify a number of terms. It is worth noting that equation (9.4) is not only valid in the steady state limit t→∞, but is also correct for all intermediate times t, which correspond to the transients which take the system from the initial equilibrium state to the final nonequilibrium steady state.

If we choose to evaluate the temperature derivative of the internal energy H₀, we can calculate the specific heat at constant volume and external field, C_v,Fe. The result is (Evans and Morriss, 1987),

(9.5)

Again the first term on the right hand side is easily recognised as the equilibrium specific heat. The second and third terms are nonlinear nonequilibrium terms. They signal the breakdown of local thermodynamic equilibrium. In the linear regime for which linear response theory is valid, they are of course both zero. The third term takes the form of a transient time correlation function. It measures the correlations of equilibrium energy fluctuations, ΔH₀(0), with the transient fluctuations in the composite-time variable, Δ( H₀(s) J(0) ). The second term can of course be rewritten as the integral of a transient time correlation function using (9.1).

Consider the Schrödinger form,

(9.6)

The thermostatted Kawasaki form for the N-particle distribution function is,

(9.7)

Since f(t) is a distribution function it must be normalised. We guarantee this by dividing the right hand side of equation (9.7) by its phase integral. If we take the initial ensemble to be canonical, we find,

(9.8)

The exponents contains a divergences due to the fact that the time average of J(-s) is nonzero. This secular divergence can be removed by multiplying the numerator and the denominator of the explicitly normalised form by exp[+βF_{e 0}∫^t ds <J(-s)>]. This has the effect of changing the dissipative flux that normally appears in the Kawasaki exponent from J(-s) to ΔJ(-s), in both the numerator and denominator. The removal of the secular divergence has no effect on the results computed in this chapter and is included here for largely aesthetic reasons.

(9.9)

The average of an arbitrary phase variable B(Γ) in the renormalized Kawasaki representation is,

(9.10)

To obtain the temperature derivative of equation (9.10) we differentiate with respect to β. This gives

(9.11)

Using the Schrödinger-Heisenberg equivalence we transfer the time dependence from the distribution function to the phase variable in each of the terms in equation (9.11). This gives

(9.12)

Substituting the internal energy for B in equation (9.12) and making a trivial change of variable in the differentiation (β→T) and integration (t-s→s), we find that the specific heat can be written as,

(9.13)

The first term gives the steady state energy fluctuations and the second term is a steady state time correlation function. As t → ∞, the only times s, which contribute to the integral are times within a relaxation time of t, so that in this limit the time correlation function has no memory of the time at which the field was turned on.

These theoretical results for the specific heat of nonequilibrium steady states have been tested in nonequilibrium molecular dynamics simulations of isothermal planar Couette flow (Evans, 1986 and Evans and Morriss, 1987). The system studied was the Lennard-Jones fluid at its triple point, (k_BT∕ε=0.722, ρσ³=0.8442). 108 particles were employed with a cutoff of 2.5σ.

Table 9.1. Lennard-Jones Specific Heat Data.^‡

Potential: Φ(r) = 4ε[(r∕σ)^-12 - (r∕σ)^-6 ].

State point: T* = 0.722, ρ* = 0.8442, γ* = 1.0, N = 108, r_c * = 2.5.

‡ Reduced units are denoted by *. Units are reduced to dimensionless form in terms of the Lennard-Jones parameters, m,σ,ε. γ ^*=∂u_x∕∂y σ(m∕ε)^1/2. ∆t*=0.004. < >_ss denotes nonequilibrium steady state average.

The steady state specific heat was calculated in three ways: from the transient correlation function expression equation (9.5), from the Kawasaki expression equation (9.13) and by direct numerical differentiation of the internal energy with respect to the initial temperature. The results are shown in the Table 9.1 below. Although we have been unable to prove the result theoretically, the numerical results suggest that the integral appearing on the right hand side of (9.5) is zero. All of our simulation results, within error bars, are consistent with this. As can be seen in the Table 9.1 the transient correlation expression for the specific heat predicts that it decreases as we go away from equilibrium. The predicted specific heat at a reduced strain rate (γσ(m∕ε)^1/2) = 1 is some 11% smaller than the equilibrium value. This behaviour of the specific heat was first observed in 1983 by Evans (Evans, 1983).

Table 9.2. Comparison of Soft Sphere Specific Heats as a function of Strain Rate^‡

Potential: Φ(r) = ε(r∕σ)^-12 . State point: T* = 1.0877, ρ* = 0.7, N = 108, r_c * = 1.5.

‡ Note: In these calculations, the transient time correlation function integral, (9.5), was assumed to be zero. Data from (Evans, 1983, Evans, 1986, Evans and Morriss, 1987).

The results obtained from the Kawasaki formula show that although the internal energy fluctuations are greater than at equilibrium, the specific heat decreases as the strain rate is increased. The integral of the steady state energy-stress fluctuations more than compensates for increase in internal energy fluctuations. The Kawasaki prediction for the specific heat is in statistical agreement with the transient correlation results. Both sets of results also agree with the specific heat obtained by direct numerical differentiation of the internal energy. Table 9.2 shows a similar set of comparisons based on published data (Evans, 1983). Once again there is good agreement between results predicted on the basis of the transient correlation formalism and the direct NEMD method.

As a final comment of this section we should stress that the specific heat as we have defined it, refers only to the derivative of the internal energy with respect to the temperature of the initial ensemble (or equivalently, with respect to the nonequilibrium kinetic temperature). Thus far, our derivations say nothing about the thermodynamic temperature ( ≡ ∂E∕∂S ) of the steady state. We will return to this subject in Chapter 10.