Since the dimension of the accessible phase space decreases to less than the ostensible 2dN dimensions, the volume of the accessible phase space, as measured from the ostensible space is zero. The entropy of a system is proportional to the logarithm of the accessible phase volume. Since that volume as determined from the ostensible phase space, is zero, the entropy will diverge to negative infinity. These simple observations explain the divergence of entropy as computed in the ostensible space. Presumably the thermodynamic entropy should be arrived at by integrating over the accessible phase space only. This would remove the apparent divergence. However the determination of the topology of the phase space which is accessible to nonequilibrium steady states is exceedingly complex. Even the dimension of the accessible space is only known approximately. Such a program for the calculation of the nonequilibrium entropy would therefore appear quite hopeless.

The fine grained entropy as computed from the ostensible phase space dimension has a number of further difficulties. From a quantum mechanical point of view, if a system such as the one depicted in Figure 10.9 is meant to represent argon, it is in violation of the Heisenberg uncertainty principle. The uncertainty principle puts an absolute limit on the degree to which a distribution function can be fractal. There is a lower limit imposed by Planck's constant, to the scale of features in that can be found in phase space. The extreme fineness of the filaments depicted in Figure 10.9 implies extreme sensitivity to external perturbations. The finer the length scale of the phase space structures, the more sensitive those structures will be to external perturbations. If the distribution function is fractal, there is no limit to the smallness of the phase space structures and therefore no limit to the sensitivity of the full distribution function to uncontrolled external perturbations. In an experiment, averaging over an ensemble of possible external fluctuations would of course wash out the fine structure below a critical length scale. The precise cut-off value would be determined by the amplitude and spectrum of the external fluctuations. This washing out of fine structure provides an ansatz for the computation of the entropy of nonequilibrium steady states.

Evans (1989) described a systematic method for computing the coarse grained entropy of nonequilibrium steady states. The coarse graining is introduced by decomposing the Gibbs (1902) entropy, into terms arising from the partial distribution functions involving correlations of successive numbers of particles. If the expansion is carried out to order N, the total number of particles in the system, the results will of course be identical to the fine-grained Gibbs entropy. The expansion has been tested at equilibrium and it has been found that for densities less than ~75% of the freezing density, the singlet and pair contributions to the entropy appear to be accurate to more than ~90%. At equilibrium, the expansion therefore appears to converge rapidly. Away from equilibrium the expansion will consist of a series of finite terms until the dimension of the partial distribution function exceeds the dimension of the accessible phase space. Once this occurs all succeeding terms will be infinite. The method yields finite terms below this dimension because all the lower dimensional integrals are carried out in the accessible phase space.

Green (1952) used Kirkwood's factorization of the N-particle distribution function to write an expansion for the entropy. If we define z-functions in an infinite hierarchy, as

(10.54)

where the various f-functions are the partial 1,2,3, .. -body distribution functions, then Green showed that Gibbs' fine grained entropy (equation 10.1.1) can be written as an infinite series,

(10.55)

Using equation (10.5.1) one can easily show that the entropy per particle is given by the following series.

(10.56)

In deriving this equation we have assumed that the fluid is homogeneous. This enables a spatial integration to be performed in the first term. This equation is valid away from equilibrium. Using the fact that at equilibrium the two body distribution function factors into a product of kinetic and configurational parts equation (10.5.3) for two dimensional fluids, reduces to,

(10.57)

where g(r₁₂) is the equilibrium radial distribution function. Equation (10.57) has been tested using experimental radial distribution function data by Mountain and Raveché (1971) and by Wallace (1987). They found that the Green expansion for the entropy, terminated at the pair level, gives a surprisingly accurate estimate of the entropy from the dilute gas to the freezing density. As far as we know prior to Evans’ work in 1989, the Green expansion had never been used in computer simulations. This was because, in the canonical ensemble, Green’s entropy expansion is non-local. In Evans’ calculations the entropy was calculated by integrating the relevant distribution functions over the entire simulation volume. A recent reformulation of (10.5.4) by Baranyai and Evans, (1989), succeeds in developing a local expression for the entropy of a canonical ensemble of systems. Furthermore the Baranyai-Evans expression for the entropy is ensemble independent.

Evans (1989) used a simulation of 32 soft discs (Φ(r) = ε(σ∕r)¹² truncated at r∕σ=1.5) to test equation (10.57) truncated at the pair level. All units were expressed in dimensionless form by expressing all quantities in terms of the potential parameters σ,ε and the particle mass m. Table 10.3, below shows some of the equilibrium data gathered for the soft disc fluid. All units are expressed in reduced form. Each state point was generated from a ten million timestep simulation run using a reduced timestep of 0.002. The energy per particle is denoted e, and the total 1 and 2-body entropy per particle is denoted by s. The entropy was calculated by forming histograms for both g(r) and f(p). These numerical approximations to the distribution functions were then integrated numerically. The radial distribution function was calculated over the minimum image cell to include exactly the long ranged, nonlocal, contributions arising from the fact that at long range, g(r) = (N-1)∕N. The equipartition, or kinetic, temperature corrected for O(1∕N) factors, is denoted by T_k. The thermodynamic temperature T_th was calculated from equation (10.57) using the thermodynamic relation, T_th=∂e∕∂s)_V. For each density the three state points were used to form a simple finite difference approximation for the derivative.

The analytical expression for the kinetic contribution to the entropy was not used, but rather this contribution was calculated from simulation data by histograming the observed particle velocities and numerically integrating the single particle contribution. The numerical estimate for the kinetic contribution to the entropy was then compared to the theoretical expression (basically the Boltzmann H-function) and agreement was observed within the estimated statistical uncertainties.

By using the entropies calculated at ρ = 0.6, 0.7 to form a finite difference approximation to the derivative ∂s∕∂ρ^-1 one can compare the pressure calculated from the relation p=T∂S∕∂V)_E, with the virial expression calculated directly from the simulation. The virial pressure at e=2.134, ρ=0.65, is 3.85 whereas the pressure calculated exclusively by numerical differentiation of the entropy is 3.72 ± 0.15. The largest source of error in these calculations is likely to be in the finite difference approximation for the various partial derivatives.

‡ The uncertainties in the entropies are ±0.005.

Away from equilibrium the main difficulty in using even the first two terms in equation (10.5.3) is the dimensionality of the required histograms. The nonequilibrium pair distribution function does not factorize into a product of kinetic and configurational parts. One has to deal with the full function of 6 variables for translationally invariant two dimensional fluid. In his work, Evans reduced the density to ρ~0.1 where the configurational contributions to the entropy should be unimportant. He evaluated the entropy of the same system of 32 soft discs, but now the system was subject to isoenergetic planar Couette flow, using the SLLOD equations of motion. In this simulation a constant thermodynamic internal energy H₀ ≡ Σp²∕2m + Φ was maintained. The thermostatting multiplier α, takes the form (see equation 5.2.3),

(10.58)

where P_xy is the xy-element of the pressure tensor.

To check the validity of our assumption that at these low densities, the configurational parts of the entropy may be ignored, he performed some checks on the equilibrium thermodynamic properties of this system. Table 10.4 shows the thermodynamic temperature computed using a finite difference approximation to the derivative, ∂e∕∂s, (e=<H₀>∕N, s=S∕N). It also shows the kinetic temperature computed using the equipartition expression. At equilibrium, the data at a reduced density of 0.1 predicts a thermodynamic temperature which is in statistical agreement with the kinetic temperature, 2.12±0.04 as against 2.17, respectively. The equilibrium data at e=2.134, ρ=0.1, gives a thermodynamic pressure of 0.22, in reasonably good agreement with the virial pressure (including both kinetic and configurational components) of 0.24. The disagreement between the thermodynamic and the kinetic expressions for both the temperature and the pressure arise from two causes; the absence of the configurational contributions, and the finite difference approximations for the partial derivatives.

Figure 10.14 shows the analogue of Figure 10.9 for a 32 particle system under shear. The nonequilibrium pair distribution function is free of the singularities apparent in the 2-particle system. The reason why it is smooth is that for 1 and 2-particle distributions in systems of many particles, one averages over all possible positions and momenta for the other N-2 particles. This averaging washes out the fine structure. These distributions even at very high strain rates, are not fractal. If the Green expansion converges rapidly we will clearly arrive at a finite value for the entropy.

Table 10.4 gives the computed kinetic contribution to the entropy as a function of energy, density and strain rate. At low densities the increased mean free paths of particles relative to the corresponding situation in dense fluids means that considerably longer simulation runs are required to achieve an accuracy comparable to that for dense fluids. The data given in table 10.4 is taken from 15 million timestep simulation runs. Away from equilibrium the strain rate tends to increase the mixing of trajectories in phase space so that the errors actually decrease as the strain rate is increased.

For a given energy and density, the entropy is observed to be a monotonically decreasing function of the strain rate. As expected from thermodynamics, the equilibrium state has the maximum entropy. Although there is no generally agreed upon framework for thermodynamics far from equilibrium, it is clear that the entropy can be written as a function, S = S(N,V,E,γ). Defining T_th as ∂E∕∂S)_V,γ , p_th as T∂S∕∂V)_E,γ and ζ_th as -T∂S∕∂γ)_E,V, we can write,

(10.59)

Some years ago Evans and Hanley (1980) proposed equation (10.59) as a generalized Gibbs relation, however, at that time there was no way of directly computing the entropy or any of the free energies. This forced Evans and Hanley to postulate that the thermodynamic temperature was equal to the equipartition or kinetic temperature, T_k ≡ 2K∕(dNk_B), for systems in d dimensions. Evans and Hanley observed that away from equilibrium, although the pressure tensor is anisotropic, the thermodynamic pressure must be independent of the manner in which a virtual volume change is performed. The thermodynamic pressure must therefore be a scalar. They assumed that the thermodynamic pressure would be equal to the simplest scalar invariant of the pressure tensor that was also consistent with equilibrium thermodynamics. In two dimensional systems they assumed that p=(P_xx+P_yy)∕2.

Since we can now calculate the coarse grained Gibbs entropy directly, we can check the correctness of these postulates. We assume, that the internal energy is given by the sum of the peculiar kinetic energy and the potential energy, that we know the system volume and strain rate and that the thermodynamic entropy is equal to the coarse grained Gibbs entropy which at low densities can be approximated by the first term of equation (10.56). Table 10.4 below shows a comparison of kinetic and thermodynamic temperatures for the 32-particle soft-disc system.

As has been known for some time (Evans, 1983), ∂T_k∕∂γ)_V,E is negative leading to a decrease in the kinetic temperature with increasing strain rate. For this low density system the effect is far smaller than has been seen for moderately dense systems. At a density of 0.1 the kinetic temperature drops by 0.3% as the shear rate is increased to unity. The precision of the kinetic temperature for these runs is about 0.01%. The thermodynamic temperature also decreases as the strain rate is increased but in a far more dramatic fashion. It decreases by 10% over the same range of strain rates. The results clearly show that away from equilibrium the thermodynamic temperature is smaller than the kinetic or equipartition temperature. As the strain rate increases the discrepancy grows larger.

Using the simulation data at e=2.134, one can estimate the thermodynamic pressure as a function of strain rate. Table 10.5 shows the finite difference approximation for the thermodynamic pressure, p_th, the hydrostatic pressure, p_tr= (P_xx+P_yy)∕2 and the largest and smallest eigenvalues of the pressure tensor p₁,p₂ respectively. As expected the hydrostatic pressure increases with shear rate. This effect, known as shear dilatancy, is very slight at these low densities. The thermodynamic pressure shows a much larger effect but it decreases as the strain rate is increased. In an effort to give a mechanical interpretation to the thermodynamic pressure we calculated the two eigenvalues of the pressure tensor. Away from equilibrium, the diagonal elements of the pressure tensor differ from one another and from their equilibrium values, these are termed normal stress effects. The eigenvalues are influenced by all the elements of the pressure tensor including the shear stress. One of the eigenvalues increases with strain rate while the other decreases and within statistical uncertainties the latter is equal to the thermodynamic pressure.

‡ Away from equilibrium the uncertainties in the entropy are ±0.005.

Evans (1989) conjectured that the thermodynamic pressure is equal to the minimum eigenvalue of the pressure tensor, that is p_th= p₂. This relation is exact at equilibrium and is in accord with our numerical results. It is also clear that if the entropy is related to the minimum reversible work required to accomplish a virtual volume change in a nonequilibrium steady state system, then p₂dV is the minimum pV work that is possible. If one imagines carrying out a virtual volume change by moving walls inclined at arbitrary angles with respect to the shear plane then the minimum virtual pV work (minimized over all possible inclinations of the walls) will be p₂dV.

Figure 10.15 shows the kinetic contribution to the entropy as a function of strain rate for the 32-particle system at an energy e=2.134 and a density ρ=0.1. The entropy seems to be a linear function of strain rate for the range of strain rates covered by the simulations. Combining these results with those from Table 10.4 allows us to compute ζ_th as a function of strain rate. For γ=0.0, 0.5, 1.0 we find that ζ_th∕N = 1.22, 1.08, and 0.91 respectively. Most of the decrease in ζ is due to the decrease in the thermodynamic temperature with increasing strain rate. We have assumed that asymptotically s is linear in strain rate as the strain rate tends to zero. It is always possible that at strain rates which are too small for us to simulate, that this linear dependence gives way to a quadratic variation.

Although these calculations are restricted to the low density gas regime, the results suggest that a sensible definition for the nonequilibrium entropy can be given. A definition, based on equation (10.56), avoids the divergences inherent in the fine grained entropy due to the contraction of the nonequilibrium phase space. At low densities this entropy reduces to the Boltzmann entropy implicit in the Boltzmann H-function. Our entropy is, for states of a specified energy and density, a maximum at equilibrium.

Defining a temperature on the basis of this entropy, indicates that far from equilibrium there is no reason to expect that the equipartition, or kinetic temperature is equal to the thermodynamic temperature. Similarly there seems to be no reason to expect that the average of the diagonal elements of the pressure tensor will be equal to the thermodynamic pressure far from equilibrium. The concept of minimum reversible virtual work, together with our numerical results suggests that the thermodynamic pressure is instead equal to the minimum eigenvalue of the pressure tensor.

It remains to be seen whether the entropy so defined, is a local maximum in nonequilibrium steady states. If this can be satisfactorily demonstrated then we will have for the first time a fundamental basis for a generalized thermodynamics of steady states far from equilibrium.