Exact residual entropies of hexagonal and cubic ice: largest eigenvalues of the layer transfer matrices $\mathbf{M}$, $\mathbf{M}\mathbf{M}^T$, and $\mathbf{M}'$

The number of proton arrangements compatible with the ice rules — the residual entropy — is known exactly only for two-dimensional models, not for real three-dimensional ice. Li, Cen, Wang, and Yang [Phys. Rev. B 110, 195414 (2024)] reformulated the residual entropies of hexagonal ice (Ih) and cubic ice (Ic) as largest eigenvalues of explicitly constructible layer transfer matrices: $\mathbf{M}$ for ice Ic, and $\mathbf{M}\mathbf{M}^T$ or an equivalent rearranged matrix $\mathbf{M}'$ for ice Ih, with $w_{\mathrm{Ic}} = \lim_{m\to\infty} \lambda_{\max,\mathbf{M}}^{1/m}$ and $w_{\mathrm{Ih}} = \lim_{m\to\infty} \lambda_{\max,\mathbf{M}\mathbf{M}^T}^{1/2m} = \lim_{m\to\infty} \lambda_{\max,\mathbf{M}'}^{1/m}$, and recovered Onsager's inequality $w_{\mathrm{Ih}} \ge w_{\mathrm{Ic}}$. The open problem is the exact evaluation of these quantities: provably exact values of $w_{\mathrm{Ic}}$ and $w_{\mathrm{Ih}}$, and the decision whether $w_{\mathrm{Ih}} = w_{\mathrm{Ic}}$ or the inequality is strict. Linked exact-solvability questions in the same framework are the frustrated two-dimensional bilayer ice models associated with $\mathbf{M}\mathbf{M}^T$ (equivalently, Ising models on stacked Kagomé bilayer lattices) and the classification of staggered vertex models that preserve thermodynamic-limit free-energy equality with their unstaggered counterparts. Later work supplies high-precision tensor-network values ($w \approx 1.5074585$ for both phases) and rigorous bounds $1.503360 \le w \le 1.540196$, but no exact result: the problem remains open.

Background

When water ice is cooled toward absolute zero, its protons do not order. Each oxygen–oxygen bond carries exactly one proton, and every oxygen keeps exactly two protons close to it, so that along the four bonds of each oxygen two protons point in and two point out — the ice rules, also called the two-in/two-out rule. Configurations of proton positions obeying these rules are called ice-ruled, and their number $W$ fixes the residual entropy $S = \ln W$ of the crystal; it is conventional to quote the configuration constant per water molecule, $w = W^{1/N_{\mathrm{H_2O}}}$, so that $S = N_{\mathrm{H_2O}} \ln w$. Pauling's 1935 counting argument gave $w = 3/2$; Onsager later proved that this value is a lower bound for any four-coordinated ice (as recounted in the works cited below); Nagle's 1966 series expansions refined the estimate to $w \approx 1.50685$ for real ice. For half a century the only exact result in the family was Lieb's 1967 solution of the two-dimensional square-ice model, $w = (4/3)^{3/2}$, obtained by transfer-matrix methods. In the transfer-matrix method, the number of allowed configurations of a system built from $n$ repeated identical units is written as $\mathrm{Tr}\, T^n$, where the matrix $T$ adds one unit and carries the boundary information between adjacent units; the thermodynamic free energy per unit is governed by the largest eigenvalue $\lambda_{\max}$ of $T$, and the method is exact when $T$ can be diagonalized in the limit of infinitely wide units.

The two common ambient phases of ice — hexagonal ice (ice Ih, the ordinary form) and cubic ice (ice Ic) — have oxygen sublattices that are different stackings of the same puckered hexagonal monolayer: cubic ABC-type stacking for Ic and alternating mirror stacking for Ih. The exact residual entropies $w_{\mathrm{Ic}}$ and $w_{\mathrm{Ih}}$ of these three-dimensional phases are not known; determining them exactly is the subject of this problem.

Li, Cen, Wang, and Yang (2024) recast the problem in a transfer-matrix form suggested by Nagle's thesis. Take a hexagonal monolayer containing $m$ water molecules; each molecule has one vertical bond leaving the layer, so there are $m/2$ lower and $m/2$ upper vertical bonds whose proton positions (one of two per bond) label the boundary conditions. Define $\mathbf{M}$ as the $2^{m/2} \times 2^{m/2}$ matrix whose element $M_{A,B}$ is the number of ice-ruled configurations of the monolayer when the lower and upper vertical-bond configurations are $A$ and $B$. Stacking $n$ such monolayers in the cubic sequence gives $W = \mathrm{Tr}\, \mathbf{M}^n$, so

$$w_{\mathrm{Ic}} = \lim_{m\to\infty} \lambda_{\max,\mathbf{M}}^{1/m},$$

while a bilayer — a monolayer together with its mirror image — is the two-dimensional periodic unit of ice Ih, so $\mathbf{M}\mathbf{M}^T$ is a transfer matrix for ice Ih:

$$w_{\mathrm{Ih}} = \lim_{m\to\infty} \lambda_{\max,\mathbf{M}\mathbf{M}^T}^{1/2m} = \lim_{m\to\infty} \lambda_{\max,\mathbf{M}'}^{1/m},$$

where $\mathbf{M}'$ is an alternative monolayer-based transfer matrix for ice Ih obtained by rearranging the elements of $\mathbf{M}$. Because the largest eigenvalue of $\mathbf{M}\mathbf{M}^T$ is the square of the spectral norm (the largest singular value) of $\mathbf{M}$, one has $\lambda_{\max,\mathbf{M}\mathbf{M}^T} \ge \lambda_{\max,\mathbf{M}}^2$, hence $w_{\mathrm{Ih}} \ge w_{\mathrm{Ic}}$ — Onsager's inequality, rediscovered by Li et al. Li et al. could not diagonalize these matrices; they obtained instead a lattice-specific lower bound $w(\mathrm{Ih}) \ge 1.501060$ and Wang–Landau estimates (a flat-histogram Monte Carlo method that counts configurations by random walks in energy) for two-dimensional models associated with $\mathbf{M}\mathbf{M}^T$.

Several reductions of these matrices are exactly solvable and anchor the problem. The sum of all elements of $\mathbf{M}$ counts the ice-ruled configurations of the free hexagonal monolayer and was solved exactly by the same group through a mapping to the antiferromagnetic Ising model on the Kagomé lattice (the lattice of corner-sharing triangles). The trace of $\mathbf{M}'$ counts the configurations of another two-dimensional ice model and was solved exactly by mapping it to a staggered six-vertex model — a vertex model, meaning arrows on edges with a weight assigned to each allowed arrow pattern at a vertex, in which the vertex weights alternate between two sublattices — with the weights staggered row by row; Li et al. proved that for this row-by-row staggering the partition function equals that of the unstaggered model in the thermodynamic limit. By contrast, the two bilayer ice models defined by the element sum and the trace of $\mathbf{M}\mathbf{M}^T$ (bilayers with open and with periodic boundary conditions on the vertical bonds) are equivalent to Ising models on stacked Kagomé bilayer lattices; these frustrated two-dimensional models are not solved — only Wang–Landau estimates exist, $w \approx 1.7878$ (open) and $w \approx 1.5311$ (periodic). The paper further poses the generalization of its row-by-row equivalence to more general staggered vertex models, in particular the staggered eight-vertex model (a vertex model allowing all arrow patterns with an even number of arrows pointing in) with row-by-row weights, and notes that applying tensor-network methods — in which the ice rule is encoded in local tensors so that the degeneracy becomes the dominant eigenvalue of a transfer operator built from them — to these matrices was in progress.

That program has since been carried out numerically. Building on tensor-network residual-entropy computations for three-dimensional frustrated magnets, Xu, Lin, and Zhang (2026) computed both entropies to six-digit accuracy — $S_{\mathrm{Ih}} = 0.4104251(6)$ and $S_{\mathrm{Ic}} = 0.4104248(14)$ per molecule, i.e. $w \approx 1.5074585$ for both phases — using a projected-entangled-pair-operator representation of the transfer operator (the transfer operator built from the local tensors described above), variational infinite projected entangled-pair states (a variational description of an infinite two-dimensional system by a single repeated local tensor), and a split corner-transfer-matrix renormalization group (an iterative algorithm that contracts such an infinite network by growing a finite environment), and showed that if the ice-Ic transfer operator is exactly normal, meaning $\hat{M}\hat{M}^T = \hat{M}^T\hat{M}$, then $S_{\mathrm{Ih}} = S_{\mathrm{Ic}}$ rigorously; numerically the operator is close to normal but with a small violation (fidelity $F_\infty = 0.99990325(28)$). Chen and Ran (2026) proved that hexagonal stacking maximizes the entropy among all periodic ice-I polytypes with uniform stacking registry and certified the rigorous bounds $1.503360 \le w \le 1.540196$ with $w_{\mathrm{Ic}} \le 1.527699$. No closed-form value, no proof or refutation of $w_{\mathrm{Ih}} = w_{\mathrm{Ic}}$, and no exact solution of the bilayer or staggered-vertex questions exist.

Problem Statement

The problem is the exact evaluation of the residual entropies of ice Ih and ice Ic in the layer-transfer-matrix formulation, together with the exact solvability of the reductions and generalizations posed in the same source paper. All components share one objective: replacing numerical estimates of the ice-rule counting problem by exact results.

(1) Primary objective — exact spectra. Determine the exact thermodynamic-limit values

$$w_{\mathrm{Ic}} = \lim_{m\to\infty} \lambda_{\max,\mathbf{M}}^{1/m}, \qquad w_{\mathrm{Ih}} = \lim_{m\to\infty} \lambda_{\max,\mathbf{M}\mathbf{M}^T}^{1/2m} = \lim_{m\to\infty} \lambda_{\max,\mathbf{M}'}^{1/m},$$

that is, closed-form (or otherwise provably exact) values of the largest eigenvalues of the ice-Ic monolayer transfer matrix $\mathbf{M}$, of the ice-Ih bilayer matrix $\mathbf{M}\mathbf{M}^T$, and of the alternative ice-Ih matrix $\mathbf{M}'$, in the sense established by the construction of Li et al. (matrix dimension $2^{m/2}$ for a monolayer of $m$ sites, $W = \mathrm{Tr}\, \mathbf{M}^n$).

(2) Equality or strict separation. Decide whether $w_{\mathrm{Ih}} = w_{\mathrm{Ic}}$ or $w_{\mathrm{Ih}} > w_{\mathrm{Ic}}$; Onsager's inequality gives only $w_{\mathrm{Ih}} \ge w_{\mathrm{Ic}}$. One sufficient route, established by Xu, Lin, and Zhang, is a proof that the ice-Ic transfer operator is exactly normal, $\hat{M}\hat{M}^T = \hat{M}^T\hat{M}$, which implies $w_{\mathrm{Ih}} = w_{\mathrm{Ic}}$; since normality is sufficient but not necessary, refuting it would not by itself settle the question.

(3) Linked exact-solvability question — bilayer reductions. Solve exactly the two frustrated two-dimensional bilayer ice models obtained from $\sum_{ij} (\mathbf{M}\mathbf{M}^T)_{ij}$ (open boundary conditions on the vertical bonds) and $\mathrm{Tr}\, \mathbf{M}\mathbf{M}^T$ (periodic boundary conditions), equivalently the Ising models on stacked Kagomé bilayer lattices; only Wang–Landau estimates ($w \approx 1.787797$ and $w \approx 1.531063$) currently exist.

(4) Linked classification question — staggered vertex models. Determine which staggered vertex models have a thermodynamic-limit free energy equal to that of the corresponding unstaggered model, extending the proved row-by-row staggered six-vertex equivalence to general staggerings, and in particular decide the case of the row-by-row staggered eight-vertex weight pattern (Eq. (40) of Li et al.).

The problem is completely solved when (1) and (2) are established; (3) and (4) are further source questions within the same transfer-matrix framework, and a submission resolving only (3) or (4) is a partial answer. Numerical methods (Wang–Landau, tensor networks) may support any answer but do not by themselves satisfy the objective, which requires exact results.

References

1. De-Zhang Li, Yu-Jie Cen, Xin Wang, Xiao-Bao Yang, Residual entropy of ice: A study based on transfer matrices, Physical Review B 110, 195414 (2024). DOI: 10.1103/PhysRevB.110.195414, arXiv:2404.13897, https://arxiv.org/abs/2404.13897 (source paper)
2. Xia-Ze Xu, Tong-Yu Lin, Guang-Ming Zhang, Equivalence of residual entropy of hexagonal and cubic ices from tensor network methods, Physical Review B 113, 214416 (2026). DOI: 10.1103/3myh-6pdf, https://doi.org/10.1103/3myh-6pdf, arXiv:2511.22477, https://arxiv.org/abs/2511.22477
3. Zhengyue Chen, Sheng Ran, Hexagonal Stacking Maximizes Proton Configurational Entropy among Ice-I Polytypes, arXiv:2608.07613 (2026), https://arxiv.org/abs/2608.07613
4. Linus Pauling, The structure and entropy of ice and of other crystals with some randomness of atomic arrangement, Journal of the American Chemical Society 57, 2680 (1935). DOI: 10.1021/ja01315a102, https://doi.org/10.1021/ja01315a102
5. J. F. Nagle, Lattice statistics of hydrogen bonded crystals. I. The residual entropy of ice, Journal of Mathematical Physics 7, 1484 (1966). DOI: 10.1063/1.1705058, https://doi.org/10.1063/1.1705058
6. Elliott H. Lieb, Residual entropy of square ice, Physical Review 162, 162 (1967). DOI: 10.1103/PhysRev.162.162, https://doi.org/10.1103/PhysRev.162.162
7. De-Zhang Li, Wei-Jie Huang, Yao Yao, Xiao-Bao Yang, Exact results for the residual entropy of ice hexagonal monolayer, Physical Review E 107, 054121 (2023). DOI: 10.1103/PhysRevE.107.054121, https://arxiv.org/abs/2211.14679
8. Laurens Vanderstraeten, Bram Vanhecke, Frank Verstraete, Residual entropies for three-dimensional frustrated spin systems with tensor networks, Physical Review E 98, 042145 (2018). DOI: 10.1103/PhysRevE.98.042145, https://doi.org/10.1103/PhysRevE.98.042145