# The 6-vertex model of hydrogen-bonded crystals with bond defects

## 1 Introduction

The 6-vertex model on a square lattice  describes hydrogen-bonded crystals in two dimensions. Historically, it was Slater  who first considered the evaluation of the residual entropy of ice, a hydrogen-bonded crystal, under the assumptions that (i) there is one hydrogen atom on each lattice edge, and (ii) there are always two hydrogen atoms near, and away from, each lattice site (the ice rule). Under these assumptions there are 6 possible hydrogen configurations at each site, and one is led to a 6-vertex model. The exact residual entropy of the “square” ice, i.e., ice on the square lattice, was obtained by Lieb , which gives rise to a numerical number surprisingly close to the experimental residual entropy of real ice (in three dimensions). The 6-vertex model is therefore an accurate description of hydrogen-bonded crystals.

In real hydrogen-bonded crystals, however, there exist bonding defects . One way through which bond defects can occur is caused by the double-well potential seen by hydrogen atoms between two lattice sites. When two hydrogens occupy the two off-center potential wells along a given lattice edge, the assumption (i) above is broken, albeit the ice rule (ii) is still intact. This leads to the fact that one, two or zero hydrogen atoms can be present on a lattice edge. Indeed, one of us  has considered this possibility in a percolation model of supercooled water. The same model has later been considered by Attard and Batchelor  who analyzed it using series analyses. More recently, Attard  reformulated the problem as a 14-vertex model with Bjerrum bond defects, and analyzed the 14-vertex model using an independent bond approximation. Here, using a somewhat different mapping, we establish the exact equivalence of the 6-vertex model with bond defects with an 8-vertex model in an external field. As a result, we are able to analyse the exact solution in a particular parameter subspace. We also discuss the general solution of the 6-vertex model with bond defects on the Bethe lattice.

## 2 Equivalence with an 8-vertex model

Consider a square lattice $``$ of $`N`$ sites under periodic boundary conditions so that there are $`2N`$ lattice edges. The lattice is hydrogen-bonded with defects such that there can be one, two or zero hydrogen atoms on each lattice edge. As the hydrogen atoms are placed off-center on the edges, we place two Ising spins $`\sigma ,\sigma ^{}`$ on each lattice edge such that $`\sigma =1`$ denotes that the site is occupied by a hydrogen and $`\sigma =1`$ the site is empty. However, the ice rule dictates that the sum of the 4 Ising spins $`\sigma _1,\sigma _2,\sigma _3,\sigma _4`$ surrounding a square lattice sites must vanish. There are altogether six ice-rule configurations as shown in Fig. 1.

More generally, we consider the ice-rule model with weights

$$\omega _1=\omega _2=a,\omega _3=\omega _4=b,\omega _5=\omega _6=c,$$
(1)   
where $`\omega _i`$ is the weight of the $`i`$th configuration shown in Fig. 1. Denote the weights (1) by $`\omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)`$ where the subscripts are indexed as shown and $`\omega _1=\omega (1,1,1,1)`$, etc. The weight $`\omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)`$ satisfies the “spin reversal” symmetry

$$\omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)=\omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)$$
(2)   
and vanishes except for the six weights given in (1). To each lattice edge containing two spins $`\sigma `$ and $`\sigma ^{}`$, introduce an edge factor $`E(\sigma ,\sigma ^{})`$ to reflect the effect of bond defects. Then, the partition function of interest is

$$Z=\underset{\sigma =\pm 1}{}\underset{\mathrm{vertices}}{}\omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)\underset{\mathrm{edges}}{}E(\sigma ,\sigma ^{}).$$
(3)   
In the 6-vertex model without bond defects we have $`E(\sigma ,\sigma ^{})=(1\sigma \sigma ^{})/2`$ so that there is precisely one hydrogen on each lattice edge. The model considered by Attard  is described by

$`E(\sigma ,\sigma ^{})`$ $`=`$ $`w_+,\sigma =\sigma ^{}=1`$ (4)   
$`=`$ $`w_{},\sigma =\sigma ^{}=1`$   
$`=`$ $`1,\sigma =\sigma ^{}.`$   
The percolation model of  is equivalent to a special case of (4) with $`w_+=w_{}=e^{2K}`$ and

$$E(\sigma ,\sigma ^{})=e^{K(\sigma \sigma ^{}+1)}=e^K(\mathrm{cosh}K)(1+z\sigma \sigma ^{}),$$
(5)   
where $`z=\mathrm{tanh}K`$. For our purposes, we shall restrict our considerations to the percolation model (5).

Attard and Batchelor  adopted an arrow representation for the hydrogen configurations which, due to the occurrence of defects, led to a 14-vertex model with Bjerrum defects. A weak-graph transformation  is then carried out for the 14-vertex model. Here, we expand the partition function directly. Substituting (5) into (3) and expanding the second product over the edges of $``$, we obtain an expansion of $`2^{2N}`$ terms. To each term in the expansion we associate a bond graph by drawing bonds on those edges corresponding to the $`z`$ factors contained in the term. This leads to a 16-vertex model on $``$. Besides an overall Boltzman factor $`(e^K\mathrm{cosh}K)^{2N}`$, the 16-vertex model has vertex weights

$$W=\underset{\sigma _1\sigma _2\sigma _3\sigma _4}{}\left(\omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)(\sqrt{z}\sigma _i)\right),$$
(6)   
where the product is taken over those incident edges with bonds. The symmetry relation (2) now implies that $`W=0`$ whenever there are an odd number of incident bonds, and the 16-vertex model becomes an 8-vertex model.

Using the bond configurations of the 8-vertex model shown in Fig. 2, it is straightforward to deduce using (6) the following vertex weights:

$`W_1`$ $`=`$ $`{\displaystyle \omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)}=2(a+b+c)`$   
$`W_2`$ $`=`$ $`z^2{\displaystyle \sigma _1\sigma _2\sigma _3\sigma _4\omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)}=2z^2(a+b+c)`$   
$`W_3`$ $`=`$ $`W_4=z{\displaystyle \sigma _2\sigma _4\omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)}=2z(ab+c)`$   
$`W_5`$ $`=`$ $`W_6=z{\displaystyle \sigma _1\sigma _2\omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)}=2z(a+bc)`$   
$`W_7`$ $`=`$ $`W_8=z{\displaystyle \sigma _2\sigma _3\omega (\sigma _1,\sigma _2,\sigma _3,\sigma _4)}=2z(abc),`$ (7)   
where $`W_i`$ is the vertex weight of the $`i`$th vertex and we have used the fact that in the non-vanishing $`\omega `$’s we have $`\sigma _1\sigma _2\sigma _3\sigma _4=1`$. Now, in an 8-vertex model configurations, vertices 3 and 4, 5 and 6, and 7 and 8, always occur in pairs and/or in even numbers. Therefore we can conveniently replace the relevant weights by their absolute values and arrive at, after dividing all weights by a common factor $`z`$,

$`W_1`$ $`=`$ $`z^1(a+b+c)`$   
$`W_2`$ $`=`$ $`z(a+b+c)`$   
$`W_3`$ $`=`$ $`W_4=|ab+c|`$   
$`W_5`$ $`=`$ $`W_6=|a+bc|`$   
$`W_7`$ $`=`$ $`W_8=|abc|.`$ (8)   
The vertex weights (8) describes an 8-vertex model in an external electric field $`h=(\mathrm{ln}z)/2`$ in both the vertical and horizontal directions .

More generally, if we allow different values of $`w_+=w_{}=w_i,i=1,2`$ in (4) and (5) for the horizontal and vertical edges respectively, and write $`h=(\mathrm{ln}z_1)/2`$ and $`v=(\mathrm{ln}z_2)/2`$, where $`z_i=(w_i1)/(w_i+1),i=1,2`$, then one arrives at an 8-vertex model with weights

$`W_1`$ $`=`$ $`e^{h+v}(a+b+c)`$   
$`W_2`$ $`=`$ $`e^{hv}(a+b+c)`$   
$`W_3`$ $`=`$ $`e^{hv}|ab+c|`$   
$`W_4`$ $`=`$ $`e^{vh}|ab+c|`$   
$`W_5`$ $`=`$ $`W_6=|a+bc|`$   
$`W_7`$ $`=`$ $`W_8=|abc|.`$ (9)   
In ensuing discussions we shall consider the general model (9).

## 3 The free-fermion solution

The free-fermion model  is defined as a particular case of the eight-vertex model in which the vertex weights satisfy the relation

$$W_1W_2+W_3W_4=W_5W_6+W_7W_8,$$
(10)   
a condition equivalent to the consideration of a noninteracting many-fermion system in an $`S`$-matrix formulation of the 8-vertex model . In the present case the free-fermion condition (10) is satisfied when either $`a=0`$ or $`b=0`$ . Without the loss of generality, we consider $`a=0`$ and $`b>c`$.

The closed expression for the free energy of the free-fermion model is well-known  and after some algebraic manipulation using results of , we obtain

$$\beta f=\underset{N\mathrm{}}{lim}\frac{1}{N}\mathrm{ln}Z=\mathrm{ln}(b+c)+\frac{1}{4\pi }_0^{2\pi }𝑑\varphi \mathrm{ln}(A+Q^{1/2}),$$
(11)   
where

$`Q`$ $`=`$ $`[\mathrm{sinh}(2v+2h)+k^2\mathrm{sinh}(2v2h)+2k\mathrm{cosh}2v\mathrm{cos}\varphi ]^2+4k^2\mathrm{sin}\varphi ^2,`$   
$`A`$ $`=`$ $`\mathrm{cosh}(2v+2h)+k^2\mathrm{cosh}(2v2h)+2k\mathrm{sinh}2v\mathrm{cos}\varphi `$ (12)   
with $`k=(bc)(b+c)`$.

The critical condition of the free-fermion model is given by 

$$W_1+W_2+W_3+W_4=2\mathrm{max}\{W_1,W_2,W_3,W_4\},$$
(13)   
where

$`W_1`$ $`=`$ $`e^{h+v}(b+c)`$   
$`W_2`$ $`=`$ $`e^{hv}(b+c)`$   
$`W_3`$ $`=`$ $`e^{hv}(bc)`$   
$`W_4`$ $`=`$ $`e^{vh}(bc).`$ (14)   
Thus, the $`h`$-$`v`$ plane is divided into four regions depending on which vertex $`1`$, $`2`$, $`3`$ or $`4`$ has the largest weight. Denoting the four regions by $`I`$, $`II`$, $`III`$ and $`IV`$ respectively, as shown in Fig. 3, the critical condition (13) can be rewritten as

$`{\displaystyle \frac{b}{c}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{tanh}v+e^{2h}}{1e^{2h}\mathrm{tanh}v}},\text{region}I,`$   
$`{\displaystyle \frac{b}{c}}`$ $`=`$ $`{\displaystyle \frac{1e^{2h}\mathrm{tanh}v}{\mathrm{tanh}v+e^{2h}}},\text{region}II,`$   
$`{\displaystyle \frac{b}{c}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{tanh}ve^{2h}}{1+e^{2h}\mathrm{tanh}v}},\text{region}III,`$   
$`{\displaystyle \frac{b}{c}}`$ $`=`$ $`{\displaystyle \frac{1+e^{2h}\mathrm{tanh}v}{\mathrm{tanh}ve^{2h}}},\text{region}IV.`$ (15)   
The critical condition (15) is plotted in Fig. $`4a4d`$ for four specific values of $`b/c`$. Generally, the free energy exhibits a logarithmic singularity at the phase boundaries with exponents $`\alpha =\alpha ^{}=0`$ . When $`b=c`$ for which $`Q`$ is a complete square, however, we have

$`\beta f`$ $`=`$ $`\mathrm{max}\{\mathrm{ln}W_1,\mathrm{ln}W_2\}`$ (16)   
$`=`$ $`\mathrm{ln}(2b)+|h+v|`$   
and the phase boundary $`h+v=0`$ separates the two frozen states $`W_1`$ and $`W_2`$.

## 4 The Bethe-like lattice

Now, we consider the 6-vertex model with bond defects on a Bethe-like lattice with plaquettes as shown in Fig. 5. The study of systems in the Bethe-like lattices is an alternative approach to the usual mean-field theory. The main features of the model under investigation will be obtained by studying the properties of the free energy.

A free energy in a region deep inside a Bethe-like lattice must be carefully defined. It cannot be obtained by directly evaluating the logarithm of the partition function in which the contribution from the outside of this region is not negligible and as result the system exhibits an unusual type of phase transitions without long-range order . Recently, a method for the surface independent free energy calculation is presented . The free energy $`f_{\mathrm{}}`$ per plaquette of our model is expressed as

$$\beta f_{\mathrm{}}=\underset{n\mathrm{}}{lim}\frac{1}{2}\left(\mathrm{ln}Z_n3\mathrm{ln}Z_{n1}\right)$$
(17)   
where $`Z_n`$ and $`Z_{n1}`$ is the partition functions of the 6-vertex model with bond defect on the Bethe-like lattice consist of $`n`$ and $`n1`$ generations respectively.

The calculation on the Bethe-like lattice is based on a recursion method . When the tree is cut at the central plaquette, it is separate into 4 branches, each of which contains 3 branches. Then the partition function of interest (3) can be written as follows:

$$Z_n=\left(\frac{\omega +1}{2}\right)^{N_b^{(n)}}\underset{\{\sigma _{0i}\}}{}\omega (\sigma _{01},\sigma _{02},\sigma _{03},\sigma _{04})g_n(\sigma _{01})g_n(\sigma _{02})g_n(\sigma _{03})g_n(\sigma _{04}),$$
(18)   
where $`N_b^{(n)}=2(3^n1)`$ is the number of bonds, $`n`$ is the number of generations and $`g_n(\sigma _{0i})`$ is in fact the partition function of a branch nearest to the $`0i`$ site.

Each branch, in turn, can be cut along any site of the first generation. then the expressions for $`g_n(\sigma _{0i})`$ can therefore be written in the form:

$$g_{n+1}(\sigma _{01})=\underset{\{\sigma _{1i}\}}{}(1+z\sigma _{01}\sigma _{13})\omega (\sigma _{11},\sigma _{12},\sigma _{13},\sigma _{14})g_n(\sigma _{11})g_n(\sigma _{12})g_n(\sigma _{14})$$
(19)   
After dividing $`g_n()`$ by $`g_n(+)`$, we obtain a recursion relation for $`x_n=g_n()/g_n(+)`$. Let us consider the case when series of solution of recursion relation converge to a stable point at $`n\mathrm{}`$, namely,

$$\underset{n\mathrm{}}{lim}x_n=x.$$

We obtain the following equation:

$$x=\frac{(1z)x^2+(1+z)x}{(1+z)x^2+(1z)x}.$$
(20)   
We are now in a position to compute the free energy per plaquette of our model. Using Eqs. (17), (18), (19) and (20), the expression for the free energy functional can be written as

$`\beta f_{\mathrm{}}`$ $`=`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \frac{1}{2}}\left(N_b^{(n)}3N_b^{(n1)}\right)\mathrm{ln}{\displaystyle \frac{\omega +1}{2}}`$   
$`+`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \frac{1}{2}}\left[\mathrm{ln}\mathrm{\Phi }(x_{(n1)})+\mathrm{ln}\mathrm{\Psi }(x_{(n)})3\mathrm{ln}\mathrm{\Psi }(x_{(n1)})\right]`$   
where

$$\mathrm{\Psi }(x)=2(a+b+c)x^2\text{and}\mathrm{\Phi }(x)=(a+b+c)^4x^4[(1+z)x+1z]^4$$

It easy to see that $`N_b^{(n)}3N_b^{(n1)}=4`$ for all $`n`$. Thus the free energy per plaquette can be finally written as

$$\beta f_{\mathrm{}}=\mathrm{ln}\frac{(\omega +1)^2}{2}(a+b+c)+\mathrm{ln}\frac{(1+z)x+1z}{2}.$$
(22)   
Together with the expression (20) for $`x`$ it gives the free energy per plaquette of the 6-vertex model with bond defects.

The equation of state (20) always has $`x=1`$ as fixed point solution. In this case the free energy per plaquette is

$$\beta f_{\mathrm{}}=\mathrm{ln}\frac{(w+1)^2}{2}(a+b+c)$$
(23)   
In the case of $`a=b=c=1`$, we recover the result obtained by Attard and Batchelor  for the six-vertex model with bond defect in the mean-field (independent vertex) framework,

$$\beta f_{\mathrm{}}=\mathrm{ln}\frac{3}{2}(w+1)^2,$$

which reduces to Pauling’s estimate when $`w=0`$ .

For the specific heat, we obtain

$`C`$ $`=`$ $`{\displaystyle \frac{d}{dT}}\left[T^2{\displaystyle \frac{d}{dT}}(\beta f_{\mathrm{}})\right]`$ (24)   
$`=`$ $`{\displaystyle \frac{2w}{(w+1)^2}}(\mathrm{ln}w)^2+{\displaystyle \frac{ab(\mathrm{ln}a/b)^2+ac(\mathrm{ln}a/c)^2+bc(\mathrm{ln}b/c)^2}{(a+b+c)^2}},`$   
or, explicitly,

$$C=\frac{1}{T^2}\frac{2K^2}{(\mathrm{cosh}K/T)^2}+\frac{\epsilon _{ab}^2e^{\epsilon _{ab}/T}+\epsilon _{ac}^2e^{\epsilon _{ac}/T}+(\epsilon _{ab}\epsilon _{ac})^2e^{(\epsilon _{ab}+\epsilon _{ac})/T}}{T^2(1+e^{\epsilon _{ab}/T}+e^{\epsilon _{ac}/T})^2}$$
(25)   
where $`\epsilon _{ab}=\epsilon _a\epsilon _b,\epsilon _{ac}=\epsilon _a\epsilon _c`$ and

$$a=e^{\epsilon _a/T},b=e^{\epsilon _b/T},c=e^{\epsilon _c/T}.$$

The specific heat versus $`T/K`$ for $`a=b`$ is plotted in Fig. 6.

## 5 Summary and discussion

In this paper we have considered a 6-vertex model of hydrogen-bonded crystals with bond defects. We have established the exact equivalence of the 6-vertex model with bond defects with an 8-vertex model in an external electric field. Using this equivalence we solve exactly our model in the free-fermion subspace. We also obtain the exact solution of the 6-vertex model with bond defects on a Bethe-like lattice.

In , one of us has used the percolation representation of the hydrogen-bonded model to argue that the specific heat of the system will increase as the temperature decreases from very high temperature. Figure 6 indeed shows such behavior. However, the specific heat in Fig. 6 does not diverge. It is of interest to calculate the specific heat for the square or diamond lattices by Monte Carlo method. We expect to find singular behavior of specific heat in such systems.

ACKNOWLEDGMENTS

This work was partly supported by the National Science Council of the Republic of China (Taiwan) under grant number NSC 89-2112-M-001-005. The work of FYW is supported in part by the National Science Foundation grant DMR-9614170. FYW thanks C. K. Hu for the hospitality at the Academia Sinica and T. K. Lee for the hospitality at the National Center for Theoretical Sciences where this work is completed.