1. Introduction

The pioneering work of Bernal and Fowler (1933) and Pauling (1935) showed that ice is sufficiently disordered to present interesting and difficult configurational problems but with a sufficiently simple kind of disorder that there was reasonable hope that definitive quantitative conclusions could be drawn, in particular, the mathematical problem of the residual entropy is sufficiently simple so as to be communicable to mathematicians and elegant enough to appeal to them. But even though the formulation of the mathematical problem is simple and elegant, the only exact calculations are for hypothetical square ice (Lieb, 1967) which is the simplest two-dimensional analogue of three-dimensional hexagonal or cubic ice. However, the two-dimensional results are extremely valuable as a testing ground for the approximate calculation methods which are used to study the three-dimensional ices. This development is reviewed in Section 4.

The theory of the dielectric constant of ice also involves a configurational problem even more difficult than the residual entropy. In the past this difficult configurational problem has tended to obscure the theoretical foundations. In particular, two of the giants in statistical mechanics, namely Onsager and Kirkwood, had quite different fundamental equations, as is reviewed in Section 2. This fundamental difference cannot be resolved by computation. Since it appears that theoreticians are about equally divided between the two theories, an experiment is proposed in Section 2 to resolve the issue.

The theory of the dielectric constant has also been confused by the existence of two quantities g and G called correlation or polarization factors. At the time of the last Ice Symposium it was believed by most researchers that there was only one correlation factor g named after Kirkwood. However, Stillinger and Cotter (1973) showed that things are not so simple. Their ideas were developed further by Nagle (1974[a]) with conclusions quite different from those of Stillinger and Cotter. This development, which is quite subtle but of great importance for the interpretation of and the feasibility of doing approximate calculations, is discussed in Section 3.

2. The Basic Dielectric-Constant Equation

The Kirkwood-Fröhlich theory of dielectric fluids (Kirkwood, 1939; Fröhlich, 1949, ch. 2) is often extended to ice and yields

The dipole moment μ is related to the dipole moment μ _v of the water molecule in the vapor and the high-frequency dielectric constant ε _∞ by μ = (ε _∞ + 2) μ _v/3. The Kirkwood correlation factor g, the calculation of which requires the use of configurational statistics, accounts for pair correlations due to short-range forces. Equation (1) was an extension of Onsager’s theory of dielectric liquids which ignored short-range forces and therefore had no g factor. However, Onsager (1936) in his original liquid-dielectric paper noted that μ in condensed hydrogen-bonded systems should not be expected to be simply related to μ _v but should be even more enhanced. For ice Onsager (Onsager and Dupuis, 1960, 1962; Onsager and Runnels, 1969) consistently used a different formula, which is analogous to one developed by Slater (1941), for KDP, namely,

to estimate μ in ice from measured values of the dielectric constant.

In the next section we will discuss the relation between g and G in Equations (1) and (2). Anticipating that discussion gives g = G for real ice with a non-zero concentration of Bjerrum defects. Thus, the Equations (1) and (2) differ by a factor . The physical models corresponding to Equations (1) and (2) both have the Bernai Fowler–Pauling ice rules, with a few Bjerrum faults, and the statistics associated with these rules determine g and G. The model in Equation (2) assumes that all interactions, including dipolar interactions, are subsumed in the ice rules, and that the effective electric field seen by a dipole is just the applied field E. In contrast, the Kirkwood–Fröhlich model modifies the E field seen by a central dipole due to the interaction with the ice, with dielectric constant ε ₀, that is outside a large imaginary sphere containing the central dipole. The basis of this modification is entirely macroscopic, but, as will be discussed in more detail in the next section, it may he thought of as the effect of long-range dipolar interactions which might be supposed to exist in addition to the ice-rule interactions.

Hollins (1964) has discussed the difference between these two approaches and suggests that the f factor should be absent. He suggests that the ice crystal consists of constant tubes of polarization along the bonds which terminate only at the surface. Thus, ?·P = 0 and the field seen by a dipole would be the applied field E with no modification, thus supporting the model in Equation (2). The omission, or at least the reduction, of long-range dipolar interactions also follows from this picture as well as from the discussion given by Slater (1941).

The reader will perhaps have surmised that this reviewer prefers Equation (2) to (1). However, too many competent researchers have believed in Equation (1) to dismiss it lightly. The best way to decide which, if either, theory is correct is by experiment. The experimental problem essentially is to determine the dipole moment μ at low frequencies. If one could apply a large enough field to orient all the protons and measure the polarization, then one would have μ = (V/N) P. Unfortunately, such a large field would enhance μ by molecular polarization, if it did not first cause dielectric breakdown. Suppose, instead, that neutron diffraction were used to measure the proton occupational probabilities α ₁ and α ₂ in the two positions on each bond in a small field E. For cubic ice with the field along a crystal axis this would give

For hexagonal ice the equation is only a little more complicated. To distinguish between the two theories would require a measurement of α ₂—α ₁ to about 10%. An experimental resolution of this problem would also be important in understanding cooperative phenomena, such as ferroelectricity, in other hydrogen-bonded crystals.

3. The Correlation Factors G and g

The Kirkwood correlation factor g is defined for isotropic systems (Stillinger and Cotter, 1973) as

where ω denotes a large sphere containing the central dipole labelled I. The order of the limits means that the entire system V is always much larger than the large sphere ω, so the correlations are never summed over all molecules j in the system. This contrasts with the definition of G given for isotropic systems by

where N is the total number of dipoles and the correlations are summed over all molecules in the system.

In this paragraph the derivation of Equation (5) will be given to help avoid any misunderstandings of its meaning and origin. Letting E be the internal electric field, the partition function for the model described by Equation (2) is

and the polarization P is

where 1_E is the unit vector in the direction of the field. Taking another derivative with respect to E, proceeding to the zero-field limit to obtain the electric susceptibility, and using Equation (2) which defines G, one obtains

The simplification of Equation (8) to (5) follows easily for cubic ice by symmetry since the sum in Equation (8) is the same when 1 _E is taken in any of the three directions parallel to the three orthogonal dipole pair types, and Equation (5) is just the average of the three sums. For lattices with lower symmetry, such as hexagonal ice, the dielectric constant and hence G (and also g) will depend on the field direction and so no such simple reduction of Equation (8) to Equation (5) is possible. In Kirkwood’s (1939) paper one would not take the step from his Equation (7) to Equation (8) for anisotropic lattices and this leaves one an equation for g similar to our Equation (8) for G. However, the calculations to be discussed in the next section show that G _╧ is very close to for hexagonal ice and so no real harm is done by thinking about Equations (4) and (5). In the case of square ice only four of the six H₂O configurations have dipole moments. Using the square lattice symmetry reduces Equation (8) to Equation (5) but with a factor of to prevent undercounting due to the configurations with zero moment. However, it must be emphasized that this difference is only a local geometrical complication similar to the anisotropy complication in hexagonal ice and does not affect the use of square ice as a testing ground for approximate methods of calculation.

In this paragraph an alternative but completely equivalent form for G is given which is more useful for some of the computations. Let n — n ₊ — n _– be the net number of dipoles pointed in the direction of the field. Then the partition function is

where W(n) is the number of states with polarization P = nμ/N. Using the maximum-term theorem

gives

Taking the derivative of Equation (11) with respect to n gives

From Equations (12) and (2) one therefore has

This formula relates G to the second derivative of the entropy with respect to the “configurational” polarization n/N. In this form it is apparent that G/3 is the same as the Γ defined by Hubmann (unpublished) and his analysis of NH₃-doped ice gives the result Γ = 0.90 ±5% or G = 2.70±0.15.

Until a few years ago it was widely supposed by most researchers that the Kirkwood factor g in Equation (4) and G in Equation (5) were the same. However, spurred by Monte Carlo calculations (Rahman and Stillinger, 1972) that give g ≈ 2.1 whereas series approximations to he discussed later gave G ≈ 3, Stillinger and Cotter (1973) convincingly argued that g and G should not be equal for perfect ice which obeys the Bernal–Fowler–Pauling ice rules completely. The reason for the possible difference between g and G comes from the fact that rigid ice rules cause very long-range correlations which fall off slowly with distance with a power-law dependence (Sutherland, 1970)

rather than the usual exponential decay

Because of Equation (14a) there may be a finite contribution outside the sphere in Equation (4) (no matter how large) which contributes to G but not to g for perfect ice. The slow decay of correlations in Equation (14a) is caused by a conservation rule in perfect ice which demands that the net orientation of the dipoles in two parallel planes be exactly equal with no fluctuations. For example, in hexagonal ice the conservation rule requires that the bond polarization in the bonds parallel to the c-axis in one layer be exactly the same for all layers.

Another reason for correlations to decay algebraically as in Equation (14a) rather than exponentially as in Equation (14b) is the long-range dipolar interactions implicitly assumed in the Kirkwood–Fröhlich theory of Equation (1) but absent in the Onsager–Slater theory of Equation (2). In fact, if Equation (14b) held for the Kirkwood–Fröhlich theory, then any large volume ω in Equation (4) would give the same value ?f g but the f factor in Equation (1) would depend upon the shape, which would be a fatal inconsistency. However, with the exception of an early calculation by Powles (1952) none of the approximate calculations which use Equations (1) and (4) actually include an electrostatic dipolar interaction in the calculation of g. This calculational procedure can be justified because the dipolar interaction over a sphere, and only over a sphere, averages to zero. While this does not guarantee that the value of g will remain unchanged when dipolar interactions are added, it does suggest that such a change will be small for a sphere. This is the basis for asserting that the long-range dipolar forces are included in the f factor in Equation (1) and that the calculation of g involves only the same ice-rule interactions that determine G in Equation (2). Thus, an answer to Hollins’ (1964) question: “Why choose a sphere?” in the Kirkwood–Fröhlich theory is that only for a sphere are the effects of the dipolar interactions and the ice-rule interactions likely to be computationally separable into the f and g factors respectively.

The reader may have noticed that the product of g = 2.1 and the factor f = 3ε ₀/(2ε ₀ + 1) in Equation (1) give nearly the same result for ε ₀ as does Equation (2) with G = 3 and so it might appear that the two theories are reconciled. Unfortunately, this is not the case, as we now show. In real ice there is a small but non-zero concentration of Bjerrum faults which relaxes the perfect ice rules and which breaks the conservation rule mentioned at the end of the last paragraph. Thus, for large enough r > r _s the correlations decay exponentially as Equation (14b) and not as Equation (14a). From the concentration of Bjerrum faults the screening distance r _s can be estimated as about 200 Å and is certainly not macroscopic, although it is large enough so that neutron scattering experiments to test the ice rules would probably not be able to differentiate real ice from perfect ice. Thus, for real ice the correlations outside a large Kirkwood sphere decay exponentially and there is only a vanishing extra contribution to Equation (5). Therefore, G = g for real ice. It should be emphasized that long-range dipolar interactions do not appreciably affect this conclusion because they do not apply to the G calculation at all and, from the discussion in the preceding paragraph, can be largely excluded from the computation of g.

Let us follow Stillinger and Cotter (1973) by introducing the parameter λ = tanh βw where 2w is the energy required to create a pair of Bjerrum faults. Then, the perfect ice limit is represented by λ = 1 and real ice is represented by λ close to but smaller than 1. Now, one can consider G and g as functions of λ, and indeed this is natural for the series-expansion calculations to be discussed in the next section and ought to be a relevant variable to any calculation. The conclusion of the last two paragraphs is that

but that it is possible that

Assuming Inequality (15b), then at least one of the functions g or G is discontinuous as a function of λ at λ = 1. This has two implications. First, if the discontinuity is significantly large, then an evaluation at λ = 1 is not useful for real ice with λ < 1. Second, most computations using approximate methods will give a smooth curve as a function of λ. For example, the Monte Carlo method requires a finite lattice and only a finite number of terms can be obtained using the series method, so in both cases the very long-range correlations will not be fully taken into account. Thus, approximation methods will necessarily be erroneous when one is calculating discontinuous functions near the discontinuity. One is likely to have much higher accuracy in approximating a continuous function, even if the function has a radius of convergence given by |λ| = 1. In the next section computational evidence is given that supports the continuity of G at λ = 1 and that therefore suggests that G is the best quantity to compute for real ice even if one wishes to know g.

4. Computational Methods

Mathematically exact computations are rare in the statistical mechanics of cooperative phenomena and ice researchers are fortunate in having some exact results for the perfect ice model for which Equation (2) applies. The residual entropy of square ice was calculated exactly to be by Lieb (1967) using the Bethe ansatz method developed by Yang and Yang (1966) for quantum magnetic problems in one dimension. From ensuing work by Sutherland and others (1967) and Lieb (1967) one can deduce (Nagle, 1973[b], 1974[a]) that

The numbers for S and G for the square lattice are expected to deviate more from the Pauling type of dendritic approximation, which gives S = R ln (3/2) and G = 3, than the numbers for the cubic or hexagonal lattices because of the existence of cycles with only four bonds in the square lattice compared to the shortest cycles of six bonds in the three-dimensional lattices, and this is borne out by the approximate calculations to be described. However, there is no reason to believe that the configurational problems are qualitatively different in square ice and in hexagonal or cubic ice. In particular, the zero-moment dipole configurations in square ice are similar to the dipole configurations in three-dimensional ices whose moments are perpendicular to the electric field. Therefore, any approximation should be put to this square-ice test.

Rahman and Stillinger (1972) performed Monte Carlo calculations for g using finite cubic ice lattices with 4 096 molecules and hexagonal ice lattices with 2 048 molecules. These calculations gave the surprising result that g ≈ 2.1 whereas the earlier work had strongly suggested g = 3. Unfortunately, no similar calculations have yet been performed on the square lattice or for the quantity G, but it is expected that they will be made before the next Ice Symposium.

By far the most successful approximate computational method for lattice statistical problems in ice has been the series-expansion method. With this technique one first develops a perturbation formalism in which each successive correction is a term in a power series. Then, one evaluates exactly as many of the terms of the series as one can. The approximation part of the method involves guessing the answer from the finite number of terms that one has computed. The simplest way to do this is to add up only the available terms. Refined approximations attempt to extrapolate the behavior of the remaining terms from the first n terms. The initial step in obtaining a useful series expansion for the residual entropy of ice was made by DiMarzio and Stillinger (1964) and the method was developed by the present author (Nagle, 1966, 1968, 1974[a]). The residual entropy was calculated for the three-dimensional ices and also for square ice (Nagle, 1966). The subsequent computation of the exact answer for square ice (Lieb, 1967) confirmed the reliability of the method and the calculated value for real ice is in good agreement with experiment. In addition to its applicability to the ice problem this series expansion method has been generalized in a form called the weak graph-series method (Nagle, 1968, 1974[b]). As such it includes older series methods used for the Ising model and has been applied to a number of other problem areas (Descamps and Coulon, 1977[a]; Wu, 1969; Nagle, 1968). Since the developments in this paragraph have been reviewed elsewhere (Nagle, 1974[b]), let us proceed to the more recent work on the dielectric constant.

At the time of the last Ice Symposium there were three groups completing dielectric calculations using the weak graph-series method (Gobush and Hoeve, GH, 1972; Stillinger and Cotter, SC, 1973; Nagle, 1973[b], 1974[a]). GH and SC were primarily interested in computing g and SC also computed the first three individual correlation functions. Both groups straight-forwardly extended the expansion method for residual entropy to obtain expansions for the correlation functions. This requires the counting of numerous graphs with two vertices of odd degree, which are all zero in the residual-entropy calculation. Even the recovery of the dendritic approximation for g or G requires a summation over an infinite subset of all such graphs. Under these circumstances it is not surprising that SC and GH restricted their attention to cubic ice only. In contrast 1 was interested primarily in G. Also, in the course of developing the series method I had noticed a trick that allowed the G series to be renormalized so that the dendritic approximation was the first term in the series and higher terms had no vertices of degree one and so were much less numerous. Therefore, it was not too difficult to treat square ice and hexagonal ice with field parallel and perpendicular to the c-axis as well as cubic ice (Nagle, 1973[b], 1974[a]).

The series are most naturally expressed in powers of the variable introduced by Stillinger and Cotter (1973), namely λ = tanh βw where w is the energy of a Bjerrum fault. For the individual correlation functions SC obtained series such as

where x = λ/3, for the nearest-neighbor correlation function of cubic ice. Gobush and Hoeve (1972) did not obtain series for individual correlations but instead added them up algebraically to obtain the first twelve terms of the following series for cubic ice,

Both GH and SC call Equation (18) the series for g. My series for G added the 13th term to Equation (18) (Nagle, 1974[a]). That Equation (18) represents both G and g is expected from Equation (15a) since they must be the same function for λ < 1.

The evaluation by SC of the correlation-function series in Equation (17) at λ = 1 gave good agreement with the Monte Carlo results of Rahman and Stillinger (1972), thereby supporting the result g ≈ 2. 1. Those correlation series, such as Equation (17), also indicated quite strongly that the radius of convergence is at λ = 1, consistent with the long-ranged nature of the correlations for the perfect-ice model. In contrast the series in Equation (18) gives a value g(1) = G(1) ≈ 3. The discrepancy here is due to the fact that the series in Equation (18) cannot possibly represent g(1). SC point out in their study of the mth-neighbor correlation functions that all coefficients of x ⁿ with n > m are negative. Thus, for a Kirkwood sphere including up to rth neighbors of a central dipole, the g series will coincide with the G series up to x ^r but for further terms the g series coefficients will change character and become increasingly negative. SC describe this as an “avalanche” of negative terms. For λ < 1, as the value of r becomes large, one will obtain g = G because all the series are inside their radii of convergence and terms beyond r become negligible. In contrast, at λ = 1 the correlation-function series are at the radii of convergence. In this case one should first sum the correlation-function series for each mth-neighbor correlation and then sum these values over m to obtain g in accordance with the Kirkwood formulation as SC did. This is a different order of summation than the one given by Equation (18). Therefore, one has a simple mathematical explanation of the difference between g(1) and G(1) because it is well known that series which are only conditionally convergent can give different values when the order of summation is changed (Knopp, 1947). I have given a simplified mathematical example that shows this explicitly (Nagle, 1974[a]).

There are two G series for hexagonal ice, one for the field parallel to the c-axis and one for the field perpendicular to the c-axis. Neither series has been written down in powers of λ; instead a less natural grouping of terms was used (Nagle, 1973[b]). Essentially all terms to order λ ¹¹ are included plus a subset of higher-order terms. It is clear that the two series are formally different, but evaluation at λ = 1 gives nearly the same values, G ₁ ≈ G _╧ ≈ 3. Thus, configurational statistics contributes negligible anisotropy to the dielectric constant. At the last Ice Symposium I stated in discussions that the 15% anisotropy measured by some groups (Humbel and others, 1953); but not by others (Wörz and Cole, 1969), requires additional intermolecular interactions in the basic ice-rule model. There is an alternative suggestion by Jaccard (1964) that the bond dipole moments may be different for the bonds parallel to the c-axis compared to the others, and if true this would avoid the need for further computations. Jaccard’s hypothesis requires the configurational polarization factor G to be. very nearly the same in both directions whereas the additional interaction hypothesis requires the c-axis G factor to be larger. It seems that neutron diffraction experiments could in principle measure the configuration polarization in the two directions and so resolve the question.

In the case of square ice the series calculation has been carried out to 11th order (Nagle, 1974[a]). The series for square ice appears more erratic than the one for cubic ice; it is also typical in magnetic series expansion work that the series for three-dimensional lattices are smoother than the series for two-dimensional lattices (Fisher, 1965). Nevertheless, various ways of analysing the series lead to values of G(1) that are well within 0.05 of the exact result (Nagle, 1974[a]). This leads one to the conclusion that G(λ) is probably continuous at λ = 1 for square ice. Of course, this is no guarantee that G(λ) is continuous at λ = 1 for hexagonal ice nor is there any direct evidence that g(1) < G(1) for square ice. But in the absence of any reason why the two- and three-dimensional lattices should not behave qualitatively the same, it is most reasonable to conclude that g(λ) is probably discontinuous at λ = 1 and that for real ice one should use g ≈ G ≈ 3 in Equations (1) and (2).

Using g = G = 3, ε ₀ = 91 at T = 273 K, ρ = 3.07 × 10²² molecules/cm³, and ε _∞ = 3.1 yields μ = 2.42 D from Equation (1) and μ = 2.93 D from Equation (2). In comparison the formula μ = (ε _∞ + 2) μ _v/3 where the vapor dipole moment μ _v = 1.83 D yields μ = 3.1 D, although, as discussed by Hollins (1964), one might think of using n ² = 1.77 in place of ε _∞ and this gives μ = 2.3 D. However, in this reviewer’s opinion the theory behind this kind of formula is not so strong for a complicated hydrogen-bonded solid such as ice that these latter dipole moments should be taken too seriously. Of much greater significance is the 21% difference between the value 2.42 D given by Equation (1) and the value 2.93 D given by Equation (2).

5. Other Problems in Configurational Statistics of Ice

This review has been primarily concerned with the theory of the dielectric constant and lo a lesser extent the entropy in the low-pressure phases of ice I. However, there are other problems for which configurational statistics are relevant.

At the last Ice Symposium, Villain and Schneider (1973) initiated the new development in configurational statistics of calculating directly the Fourier transform χ(q) of the correlation functions so as to compare theory directly with elastic neutron scattering experiments. They used a random-walk approximation which required the alteration of some of the graphical contributions. Very recently, Descamps and Coulon (1977[b]) have applied the weak-graph-series method directly to χ(q) and have obtained series for the hexagonal lattice up to 8th order with no graphical alteration. This expansion requires the counting of graphs with two vertices of odd degree, similar to the correlation-function calculations. It would appear to this reviewer to be similar in spirit to the calculation of the dielectric G in that it automatically sums the series for all correlations up to a certain order rather than summing over as many orders as possible for each correlation and then summing over correlations. This approach seems sensible and might avoid some of the problems associated with singularities in the correlation functions. The results seem to be in qualitative agreement with the earlier work of Villain and Schneider (1973) and the experiments (Axe and Hamilton, unpublished).

Also, high-pressure forms of ice are a challenging problem area, first to construct reasonable models and then to try to do appropriate configurational statistics. Since ice III and ice VII order protonically at low temperatures (Kamb, 1973), and even ice VI has a temperature-dependent g factor (Johari and Whalley, 1976), it is clear that the models will require additional interactions compared with the ice I models and not just different geometries. Since the nature of these interactions is not obvious, Kamb (1973) proposed that quantities such as the configurational entropy be calculated directly from the measured entropy probabilities of the proton configuration, thus by-passing a basic theory, and such a formula of the same approximate type as Pauling’s (1935) dendritic formula was produced (Nagle, 1973[a]). However, progress with other hydrogen-bonded systems which order at low temperatures such as copper formate tetrahydrate (Allen, 1974) and tin chloride dihydrate (Salinas and Nagle, 1974) make one optimistic about the prospects of dealing with the high-pressure forms of ice at a more basic theoretical level.