Thermostats were first introduced as an aid to performing nonequilibrium computer simulations. Only later was it realised that these devices have a fundamental role in the statistical mechanics of many-body systems. The first deterministic method for thermostatting molecular dynamics simulations was proposed simultaneously and independently by Hoover and Evans (Hoover, Ladd and Moran, 1982, and Evans, 1983). Their method employs a damping or friction term in the equations of motion. Initially the use of such damping terms had no theoretical justification. Later it was realised (Evans, Hoover, Failor, Moran and Ladd, 1983) that these equations of motion could be derived using Gauss' principle of least constraint (§3.1). This systematised the extension of the method to other constraint functions.

Using Gauss' Principle (Chapter 3), the isokinetic equations of motion for a system subject to an external field can be written as,

(5.19)

This is the thermostatted generalisation of equation (5.1) where the thermostatting term αp_i has been added. In writing these equations we are assuming:

In order to know that these three conditions are valid, we must know quite a lot about the possible flows induced in the system by the external field. This means that if we are considering shear flow for example, the Reynolds number must be small enough for laminar flow to be stable. Otherwise we cannot specify the streaming component of a particles motion (C_i must contain the local hydrodynamic flow field u(r,t)) and we cannot expect condition (1) to be valid.

The isokinetic expression for the multiplier is easily seen to be ,

(5.20)

It is instructive to compare this result with the corresponding field free multiplier given in (3.32). It is important to keep in mind that the expression for the multiplier depends explicitly on the external field and therefore on time. This is why we define the time and field independent phase variables α₀,α₁ .

It is easy to show that if Gauss' Principle is used to fix the internal energy H_{0 ,} then the equations of motion take on exactly the same form (Evans, 1983), except that the multiplier is,

(5.21)

It may seem odd that the form of the field dependent equations of motion is independent of whether we are constraining the kinetic or the total energy. This occurs because the vector character of the constraint force is the same for both forms of constraint (see §3.1). In the isoenergetic case it is clear that the multiplier vanishes when the external field is zero. This is as expected since in the absence of an external field, Newton's equations conserve the total energy.

Gaussian thermostats remove heat from the system at a rate,

(5.22)

by applying a force of constraint which is parallel to the peculiar velocity of each particle in the system.

We will now discuss the equilibrium properties of Gaussian isokinetic systems in more detail. At equilibrium the Gaussian isokinetic equations become,

(5.23)

with the multiplier given by equation (5.20) with F_e=0. Clearly the average value of the multiplier is zero at equilibrium with fluctuations in its value being precisely those required to keep the kinetic energy constant. Following our assumption that the initial value of the total linear momentum is zero, it is trivial to see that like the kinetic energy, it is a constant of the motion.

The ergodically generated equilibrium distribution function f_T(Γ), can be obtained by solving the Liouville equation for these equations of motion. It is convenient to consider the total time derivative of f. From the Liouville equation (3.34), we see that,

(5.24)

In computing the final derivative in this equation we get 3N identical intensive terms from the 3N derivatives, α(∂∕∂p_i)⋅p_i. We also get 3N terms from p_i⋅∂α∕∂p_i which sum to give -α. Since we are interested in statistical mechanical systems we will ignore terms of relative order 1∕N, in the remaining discussion. It is certainly possible to retain these terms but this would add considerably to the algebraic complexity, without revealing any new physics. This being the case, equation (5.24) above becomes,

(5.25)

From (5.24) it is can be shown that ,

(5.26)

or,

(5.27)

Integrating both sides with respect to time enables us to evaluate the time independent equilibrium distribution function,

(5.28)

where the constant, β=3N∕2K₀. We call this distribution function the isokinetic distribution f_T (Evans and Morriss, 1983). It has a very simple form: the kinetic degrees of freedom are distributed microcanonically, and the configurational degrees of freedom are distributed canonically. The thermodynamic temperatures (∂E∕∂S)_N,V=T of these two sub systems are of course identical.

If one retains terms of order 1∕N in the above derivation, the result is the same except that β=(3N-4)2K₀. Such a result could have been anticipated in advance because in our Gaussian isokinetic system four degrees of freedom are frozen, one by the kinetic energy constraint, and three because the linear momentum is fixed.

One can check that the isokinetic distribution is an equilibrium solution of the equilibrium Liouville equation. Clearly df_T∕dt≠0. As one follows the streaming motion of an evolving point in phase space Γ(t), the streaming derivative of the co-moving local density is,

(5.29)

This is a direct consequence of the fact that for a Gaussian isokinetic system, phase space is compressible. It is clear however, that in the absence of external fields ⟨df_T∕dt⟩=0, because the mean value of must be zero. If we sit at a fixed point in phase space and ask whether, under Gaussian isokinetic dynamics, the isokinetic distribution function changes, then the answer is no. The isokinetic distribution is the equilibrium distribution function. It is preserved by the dynamics. Substitution into the Liouville equation gives,

(5.30)

The proof that the last two terms sum to zero is easily given using the fact that, β=3N∕2K₀ and that K=∑p²∕2m is a constant of the motion.

(5.31)

If the equilibrium isokinetic system is ergodic, a single trajectory in phase space will eventually generate the isokinetic distribution. On the other hand a single isokinetic trajectory cannot ergodically generate a canonical distribution. We can however, ask whether isokinetic dynamics will preserve the canonical distribution. If we integrate the equations of motion for an ensemble of systems which are initially distributed canonically, will that distribution be preserved by isokinetic dynamics? Clearly,

(5.32)

is not identically zero. In this expression K is a phase variable and not a constant, and is only equal to zero on average. K would only be a constant if all members of the ensemble had identical kinetic energies. The mean value of 3N∕2K is of course β.

Consider the time derivative of the canonical average of an arbitrary extensive phase variable, B, where the dynamics is Gaussian isokinetic.

(5.33)

The time derivative of the ensemble average is,

(5.34)

where ∆K≡K-⟨K⟩=K-K₀. Equation (5.34) can be written as the time derivative of a product of three extensive, zero-mean variables.

(5.35)

In deriving these equations we have used the fact that , and that the ensemble average of the product of three extensive, zero mean phase variables is of order N, while K₀=⟨K⟩ is extensive.

The above equation shows that although B is extensive, the change in ⟨B(t)⟩ with time, (as the ensemble changes from canonical at t=0, to whatever for the Gaussian isokinetic equations generate as t→∞) is of order 1 and therefore can be ignored relative to the average of B itself. In the thermodynamic limit the canonical distribution is preserved by Gaussian isokinetic dynamics.

The Gaussian thermostat generates the isokinetic ensemble by a differential feedback mechanism. The kinetic temperature is constrained precisely by setting its time derivative to be zero. Control theory provides a range of alternative feedback processes. After the Gaussian thermostat was developed, Nosé (1984a,b) utilised an integral feedback mechanism. As we will see the Nosé thermostat, especially after a simplifying reformulation by Hoover (1985), provides a simple and direct way of ergodically generating the canonical ensemble.

The original Nosé method considers an extended system with an additional degree of freedom s, which acts like an external reservoir, interacting with the system of interest by scaling all the velocities of the particles, . The new potential energy that Nosé chose to associate with this new degree of freedom was (g+1)k_BT1ns, where g is related to the number of degrees of freedom of the system and T is the desired value of the temperature. It is essentially the choice of the potential for s which leads to dynamics which generate the canonical ensemble.

The equivalent Hoover formulation of the Nosé thermostat uses equations of motion with the same form as the Gaussian equations. The difference being that the thermostatting multiplier α, is determined by a time integral of the difference of the actual kinetic temperature from its desired value. All present applications of the Nosé thermostat use the Hoover reformulation rather than the original, somewhat cumbersome approach.

The Nosé Hamiltonian for the extended system is,

(5.36)

where Q is effectively the mass associated with the heat bath (s is dimensionless so the parameter Q does not have the dimensions of mass). The equations of motion generated by this Hamiltonian are

(5.37)

If we eliminate the variable p_s from the equations of motion obtaining instead of the last two equations a single second order differential equation for s,

(5.38)

If the system is at equilibrium, the average force on the s coordinate must be zero, so that

(5.39)

Suppose we interpret the time appearing in (5.37) to be a non-Galilaean fictitious time, and the real velocities to be v_i=s(dq_i∕dt). The instantaneous temperature is related to ∑_im_iv²_i, and its time averaged value is equal to (g+1)k_BT, where g+1 is the number of degrees of freedom. This is consistent with a non-Galilaean time average being given by

(5.40)

This is an unusual definition of a time average as it implies that equal intervals non-Galilaean time dt, correspond to unequal intervals in real time of dt∕s. Large values of s can be understood as corresponding to a rapid progress of fictitious time t. Division by s in the time averages appearing in (5.40) cancels out the uneven passage of fictitious time restoring the averages to their correct Galilean values.

To calculate the equilibrium distribution function corresponding to the Nosé Hamiltonian we use the fact that for an ergodic system, the equilibrium distribution function for Nosé's extended system is microcanonical. The microcanonical partition function for the extended system is,

(5.41)

(5.42)

where q and p are 3N-dimensional vectors, q≡(q₁,...,q_N) and p≡(p₁,...,p_N). If we change variables from p to p', where p'_i=p_i∕s for all i, then

(5.43)

where H'₀ is the usual N particle Hamiltonian , (the prime indicates that H'₀ is a function of p'). The integral over s can be performed as the only contributions come from the zeros of the argument of the delta function. If , then G has only one zero, that is

(5.44)

Using the identity δ(G(s))-δ(s-s₀)∕G'(s) it is easy to show that performing the integral over s gives

(5.45)

The integral over p_s is the infinite integral of a Gaussian and the result is

(5.46)

If we choose g=3N then this partition function is simply related to the canonical partition function

(5.47)

If the variables q, p, s, p_s are distributed microcanonically then variables p' and q are canonically distributed. The notion of non-Galilaean time makes this formulation of the Nosé thermostat rather cumbersome to use and difficult to interpret.

The next step in the development of this method was made by Hoover (1985) who realised that if one's interest lies solely in computing averages over q, p' in real time then you may as well rewrite the equations of motion in terms of q, p' and real time, t', and eliminate the p, s, p_s, t variables entirely. He used the time transformation

(5.48)

so that dt'=dt∕s, to rewrite the Nosé equations of motion as

(5.49)

where K₀ is the value of the kinetic energy corresponding to the required value of the temperature K₀=(g+1)k_BT∕2, K(p') is the instantaneous value of the kinetic energy, τ is a relaxation time which is related to the mass of the s degree of freedom (τ²=Q∕2K₀) and ζ=p_s∕Q. The motion of the system of interest can now be determined without reference to s. It is an irrelevant variable which can be ignored. The variable dζ∕dt'is a function of p' only, so the complete description of the system can be given in terms of the variables q and p'.

An important result, obtained from this time transformation by Evans and Holian (1985), is that time averages in terms of the variables q , p' and t' take their usual form, that is

(5.50)

To obtain this result we start by considering the Nosé-Hoover phase variable Liouvillean iL_NH(q,p's,p_s;t') and relating it to the Nosé Liouvillean iL_N(q,p's,p_s;t').

(5.51)

Using the results:

(5.52)

and

(5.53)

equation (5.51) becomes

(5.54)

If A is an arbitrary phase variable then the Liouvillean describes the rate of change of A. If we consider A to be a function of q and p then the rate of change of A with respect to time t is

(5.55)

Since iL_N contains no explicit time dependence, integrating with respect to time gives

(5.56)

In a similar fashion we can consider A to be function of q and p'. In that circumstance it is natural to ask for the value of A at t'.

(5.57)

Now A is function of the reduced phase space only, so the dependence on s and p_s can be ignored. These two different representations of the phase variable can be equated. To do this consider the time transformation (5.48). It implies,

(5.58)

So that dτ'=dτ∕s, and

(5.59)

Using (5.59) and the time transformation (5.48) we find that T'=∫^T'₀dt'=∫^T₀dt∕s so that we can rewrite the time average in the usual form,

(5.60)

So the time average over t is equal to the time average over t'. Using the variables q, p' and t' the time average over equal intervals of t' takes the usual form. The time average over q, p and t however, involves the scaling variable s, or equivalently a time average over unequal intervals of t.

One can of course dispense with the original form of Nosé's equations entirely. There is now no need to consider the notion of non-Galilaean time. We simply repeat the derivation we gave for the isokinetic distribution based on the Gaussian isokinetic equations of motion, for the Nosé-Hoover equations. Since there is no need to refer to non-Galilaean time we refer to q, p', t' simply as, q, p, t (dropping the prime). The N particle distribution function f(Γ,ζ) generated by the Nosé-Hoover equations of motion can be obtained by solving the Liouville equation for the equations of motion written in terms of q, p and t. It is convenient to consider the total time derivative of f(Γ,ζ) which from the Liouville equation is

(5.61)

From the equations of motion (5.49), dropping the primes, it is easy to see that dζ∕dt is a function of q and p, and hence independent of ζ. The only nonzero contribution to the right hand side comes from the p dependence of dp∕dt, so that

(5.62)

Consider the time derivative of the quantity H₀+½Qζ²

(5.63)

If we take g=3N-1 then we find that

(5.64)

So that the equilibrium distribution function is the extended canonical distribution f_c,

(5.65)

In the Hoover representation of the equations of motion, the scaling variable s has essentially been eliminated so the number of degrees of freedom of the system, changes from 3N+1 to 3N and g changes from 3N to 3N-1.