In the first few sections of this chapter we have given a description of the mechanics of individual N-particle systems. The development which follows describes an ensemble of such systems; that is an essentially infinite number of systems characterised by identical dynamics and identical state variables (N, V, E or T etc.) but different initial conditions, (Γ(0)). We wish to consider the average behaviour of a collection of macroscopically identical systems distributed over a range of initial states (microstates). In generating the ensemble we make the usual assumptions of classical mechanics. We assume that it is possible to know all the positions and momenta of an N particle system to arbitrary precision at some initial time, and that the motion can be calculated exactly from the equations of motion.
The ensemble contains an infinite number of individual systems so that the number of systems in a particular state may be considered to change continuously as we pass to neighbouring states. This assumption allows us to define a density function f(Γ,t), which assigns a probability to points in phase space. Implicit in this assumption is the requirement that f(Γ,t), has continuous partial derivatives with respect to all its variables, otherwise the phase density will not change continuously as we move to neighbouring states. If the system is Hamiltonian and all trajectories are confined to the energy surface then f(Γ,t) will not have a continuous partial derivatives with respect to energy. Problems associated with this particular source of discontinuity can obviously be avoided by eliminating the energy as a variable, and considering f(Γ,t) to be a density function defined on a surface of constant energy (effectively reducing the dimensionality of the system). However it is worth pointing out that other sources of discontinuity in the phase space density, may not be so easily removed.
To define a distribution function for a particular system we consider an ensemble of identical systems whose initial conditions span the phase space specified by the macroscopic constraints. We consider an infinitesimal element of phase space located at Γ ≡ (q,p). The fraction of systems δN, which at time t have coordinates and momenta within δq, δp of q, p is used to define the phase space distribution function f(q,p,t), by
(3.36)
The total number of systems in the ensemble is fixed, so integrating over the whole phase space we can normalise the distribution function,
(3.37)
If we consider a small volume element of phase space, the number of trajectories entering the rectangular volume element δq δp through some face will in general be different from the number which leave through an opposite face. For the faces normal to the q1-axis, located at q1, and q1+δq1, the fraction of ensemble members entering the first face is
Similarly the fraction of points leaving through the second face is
Combining these expressions gives the change in δN due to fluxes in the q1 direction
(3.38)
Summing over all coordinate (and momentum) directions gives the total fractional change δN as
(3.39)
Dividing through by the phase space volume element δq δp we obtain the rate of change in density f(q,p), at the point (q,p),
(3.40)
Using the notation, Γ = (q,p) = (q1,q2,...,q3N,p1,p2,...,p3N) for the 6N-dimensional phase point, this may be written as
(3.41)
This is the Liouville equation for the phase space distribution function. Using the streaming or total time derivative of the distribution function, we can rewrite the Liouville equation in an equivalent form as,
(3.42)
This equation has been obtained without reference to the equations of motion. Its correctness does not require the existence of a Hamiltonian to generate the equations of motion. The equation rests on two conditions: that ensemble members cannot be created or destroyed and that the distribution function is sufficiently smooth that the appropriate derivatives exist. Λ(Γ) is called the phase space compression factor since it is equal to the negative time derivative of the logarithm of the phase space distribution function.
(3.43)
The Liouville equation is usually written in a slightly simpler form. If the equations of motion can be generated from a Hamiltonian, then it is a simple matter to show that Λ(Γ)=0. This is so even in the presence of external fields which may be driving the system away from equilibrium by performing work on the system.
(3.44)
The existence of a Hamiltonian is a sufficient, but not necessary condition for the phase space compression factor to vanish. If phase space is incompressible then the Liouville equation takes on its simplest form,
. (3.45)
The following sections will be devoted to developing a formal operator algebra for manipulating the distribution function and averages of mechanical phase variables. This development is an extension of the treatment given by Berne (1977) which is applicable to Hamiltonian systems only. We will use the compact operator notation
(3.46)
for the Liouville equation, equation (3.41). The operator iL is called the distribution function (or f-) Liouvillean. Both the distribution function f, and the f-Liouvillean are functions of the initial phase Γ. We assume that there is no explicit time dependence in the equations of motion (time varying external fields will be treated in Chapter 8). Using this notation we can write the formal solution of the Liouville equation for the time dependent N-particle distribution function f(Γ,t) as
(3.47)
where f(Γ,0), is the initial distribution function. This representation for the distribution function contains the exponential of an operator, which is a symbolic representation for the infinite series of operators. The f -propagator is defined as,
(3.48)
The formal solution given above can therefore be written as
(3.49)
This form makes it clear that the formal solution derived above is the Taylor series expansion of the explicit time dependence of f(Γ,t), about f(Γ,0).
We will need to consider the time evolution of functions of the phase of the system. Such functions are called phase variables. An example would be the phase variable for the internal energy of a system, H0(Γ)= ∑ip2i∕2mi +Φ(q1,...qn). Phase variables by definition, do not depend on time explicitly, their time dependence comes solely from the time dependence of the phase Γ. Using the chain rule, the equation of motion for an arbitrary phase variable B(Γ) can be written as
(3.50)
The operator associated with the time derivative of a phase variable iL(Γ) is referred to as the phase variable (or p-) Liouvillean. The formal solution of this equation can be written in terms of the p -propagator, eiLt. This gives the value of the phase variable as a function of time
(3.51)
This expression is very similar in form to that for the distribution function. It is the Taylor series expansion of the total time dependence of B(t), expanded about B(0). If the phase space compression factor Λ(Γ) is identically zero then the p-Liouvillean is equal to the f-Liouvillean, and the p-propagator is simply the adjoint or Hermitian conjugate of the f-propagator. In general this is not the case.
In this section we will derive some of the more important properties of the Liouville operators. These will lead us naturally to a discussion of various representations of the properties of classical systems. The first property we shall discuss relates the p-Liouvillean to the f-Liouvillean as follows,
(3.52)
This is true for an arbitrary distribution function f(0). To prove this identity the LHS can be written as
(3.53)
The boundary term (or surface integral) is zero because f(S)→0 as any component of the momentum goes to infinity, and f can be taken to be periodic in all coordinates. If the coordinate space for the system is bounded then the surface S is the system boundary, and the surface integral is again zero as there can be no flow through the boundary.
Equations (3.52 & 53) show that L, L are adjoint operators. If the equations of motion are such that the phase space compression factor, (3.43), is identically zero, then obviously L=L and the Liouville operator is self-adjoint, or Hermitian.
We can calculate the value of a phase variable B(t) at time t by following B as it changes along a single trajectory in phase space. The average ⟨B(Γ(t))⟩ can then be calculated by summing the values of B(t) with a weighting factor determined by the probability of starting from each initial phase Γ. These probabilities are chosen from an initial distribution function f(Γ,0). This is the Heisenberg picture of phase space averages.
(3.54)
The Heisenberg picture is exactly analogous to the Lagrangian formulation of fluid mechanics; we can imagine that the phase space mass point has a differential box dΓ surrounding it which changes shape (and volume for a compressible fluid) with time as the phase point follows its trajectory. The probability of the differential element, or mass f(Γ)dΓ remains constant, but the value of the observable changes implicitly in time.
The second view is the Schrödinger, or distribution based picture (Fig. 3.3).
(a) The Heisenberg picture: 
(b) The Schrödinger picture: 
; 
the local density of ensemble representatives.
In this case we note that ⟨B(t)⟩ can be calculated by sitting at a particular point in phase space and calculating the density of ensemble points as a function of time. This will give us the time dependent N-particle distribution function f(Γ,t). The average of B can now be calculated by summing the values of B(Γ) but weighting these values by the current value of the distribution function at that place in phase space. Just as in the Eulerian formulation of fluid mechanics, the observable takes on a fixed value B(Γ) for all time, while mass points with different probability flow through the box.
(3.55)
The average value of B changes with time as the distribution function changes. The average of B is computed by multiplying the value of B(Γ), by the probability of find the phase point Γ at time t, that is f(Γ,t).
As we have just seen these two pictures are of course equivalent. One can also prove their equivalence using the Liouville equation. This proof is obtained by successive integrations by parts, or equivalently by repeated applications of equation (3.52). Consider
(3.56)
One can unroll each p-Liouvillean in turn from the phase variable to the distribution function (for the first transfer we consider (iL)n-1B to be a composite phase variable) so that equation (3.56) becomes,
This is essentially the property of phase and distribution function Liouvilleans which we have already proved, applied to nth Liouvillean. Repeated application of this result leads to
So finally we have the result,
(3.57)
The derivation we have used assumes that the Liouvillean for the system has no explicit time dependence. (In Chapter 8 we will extend the derivation of these and other results to the time dependent case.) Our present derivation make no other references to the details of either the initial distribution function, or the equations of motion for the system. This means that these results are valid for systems subject to time independent external fields, whether or not those equations are derivable from an Hamiltonian. These results are also independent of whether or not the phase space compression factor vanishes identically.
A final point that can be made concerning the Schrödinger and Heisenberg pictures is that these two ways of computing phase averages by no means exhaust the range of possibilities. The Schrödinger and Heisenberg pictures differ in terms of the time chosen to calculate the distribution function, f(Γ,t). In the Heisenberg picture that time is zero while in the Schrödinger picture the time is t. One can of course develop intermediate representations corresponding any time between 0 and t (eg. the interaction representation).