Hyperhomogeneity of hydrogen-disordered ice and its origin: the residual entropy compatible with the disparity in hydrogen bond energy

The residual entropy is one of the most crucial properties for the existence of a large number of ice polymorphs. The residual entropy has been estimated by Pauling assuming that there is no large difference between the hydrogen bond energies in ice. This simple model accurately predicts the entropy change of the phase transition between a hydrogen-disordered ice phase and its hydrogen-ordered counterpart. This fact is, however, incompatible with another fact that the difference in the pair interaction energies involved in hydrogen bonds in an ice phase can be larger than the thermal energy of a few kJ/mol. Here we propose a mechanism that reconciles them by considering the equality of the binding energy in each molecule rather than the pair interaction energy of the proximate pair. The topological feature of ice, called the ice rules, allows us to replace the interactions of a water molecule with the other individual molecules by that with the collections of the dipoles represented by directed cycles consisting of O-H vectors. This resummation reveals that molecular environments in ice are extremely homogeneous thereby providing a solid basis for Pauling's model. difference entropy between the phases, called residual entropy, is successfully calculated by Pauling’s the mutual orientation of all molecules irrespective the difference in the pair interactions. indicating that the binding energy of a water molecule in ice is mostly determined by the interactions with the dipole vectors that share the same cycles with the central water molecule. As shown by the blue curves, the SD of the interaction energy decreases with increasing r when the molecule-molecule interactions are accumulated. In contrast, as shown by the orange curves, the SD increases slightly with r beyond 0.7 nm when the molecule-cycle interactions are accumulated. This suggests that there are no correlations even at r = 0.7 nm between molecule-cycle interactions.


Hyperhomogeneity of the binding energy
Computer modelling allows us to estimate the intermolecular interaction between water molecules, which is capable of reproducing large differences in hydrogen bond energy in ice form. [4][5][6] We show the distribution functions of the interaction energy of a water molecule with nearest neighbors in the leftmost panels of Fig. 1a for various ice polymorphs at 0 K. The standard deviation (SD) of the interaction energy is 7 to 14 kJ mol -1 for ice Ih, III, V, and VI, and ~ 20 kJ mol -1 for ice VII. These SD values seem too large to consider that all water molecules are energetically equivalent compared with the magnitude of the hydrogen-bonding energy, ~ 25 kJ mol -1 . The large energy discrepancies in ice manifest the experimental observation and theoretical calculation that the molecular positions are irregularly deviated from the lattice point due to the influence of the orientation of surrounding molecules, which is known as site-disorder of ice. [6][7][8][9] The cumulative sum of the interaction energies of a water molecule with surrounding molecules is defined as , (1) where the pair interaction energy between water molecules i and j and the summation runs over all the molecules within a distance, r, from the central molecule, i. (∞) is called the binding energy. The cumulative interactions, I i (r), of various ice types are plotted in the main panels of Figure 1a. In general, the standard deviation of the sum of random and statistically independent values must increase with increasing the number of the accumulated samples (see Figure S1 of SI). The width of the distribution decreases with increasing r in Figure 1a. This indicates that there are some correlations between the interactions with nearest neighbors and the interactions with distant molecules as if the molecules in the middle to long range are oriented in such a way that they counteract the deviation by the interaction energy with the proximate environment. The SD of I i (r) at distance r, σ(r), converges to a very small value at a long distance of ~ 5 nm for all ice phases as shown in Fig 1b. The correlation is quite long-ranged for such a molecular system. The behavior of ice shown in Figure 1 can be called hyperhomogeneity.

Representation by a directed graph
We explore the mechanism of the hyperhomogeneity in the hydrogen-disordered ice only in terms of Coulombic interactions between point charges located on atoms for simplicity. In ice, every water molecule accepts two hydrogen bonds and donates two hydrogen bonds. In graph theory, such a network can be represented by a simple regular directed graph with in-and out-degrees of two. A directed edge of the graph corresponds to a hydrogen bond and a vertex represents the position of the negative charge in a water molecule. 5 A schematic diagram of the twodimensional model ice is shown in Figure 2a. There are numerous directed cycles in the graph of ice. A directed cycle is a cycle whose edges are traversed in the same direction. (In two dimensions, they are limited to clockwise or counterclockwise directions.) In hydrogen-disordered ice, the variety of the molecular arrangements in ice can be characterized by the spatial distribution of the directed cycles, just like the case of ice polymorphism where the structures of ice are characterized by the spatial arrangements of the undirected cycles (rings) in the hydrogen-bond network. 10 We focus on the directed cycles to examine the effect of hydrogen-disordering. All 4-membered and 6membered directed cycles in the model ice are shown in Figures 2b and 2c. One can find many other larger cycles in the graph. Here we deal with irreducible cycles only; we get rid of any cycles that pass through the same vertex 4 more than once because they are reducible to smaller cycles. (Figure 2b) We also eliminate the cycle that spans the entire cell, that is essential for an ice to be polarized. 11

Interaction energy with cycles
Each water molecule can also be viewed as being made up of two oxygen-to-hydrogen dipoles. In this paper, we will take this view to treat the two dipoles in a water molecule separately. The moments of all dipoles are identical. In general, the sum of a series of vectors that make up a ring is zero. Because the hydrogen bonds in each ice are all roughly equal in length and the O-H dipoles are parallel to the O-O vectors, the net dipole moment of a directed cycle of hydrogen bonds is approximately zero. The interaction between a water molecule and a directed cycle in ice can be reduced to the dipole-quadrupole interaction, which is much short-ranged compared with the dipole-dipole interaction. 5 We attempt to remove all dipole-dipole interactions by introducing the directed cycles. It can be achieved by an appropriate resummation of the interactions by the simplest and most tractable way. We enumerate all the cycles not larger than n-membered and make a set of them. The minimum size of the cycle depends on the type of ice, and the maximum size, n, is determined by the following procedure.
Let the set of all directed cycles up to n-membered ones in a given digraph be C = {c 1 , c 2 , ···}. A subset of C consists of cycles passing through the edge e i is denoted by C i . If one sets n appropriately to find a sufficiently large number of cycles, the weights w j for the cycle c j are determined such that , (3) exactly for all edges. (Figure 2d) To this end, all edges must belong to one or more cycles. In the case of Figure 2, all 4-membered and 6-membered cycles are sufficient to incorporate all the edges in the graph and to determine the weights. We increase n one by one until appropriate w is obtained. For ice Ih, n=12 is sufficient in most cases (it also depends on the molecular arrangements in the ice structure). A directed graph of ice is expressed by a linear combination of directed cycles contained in C. (The procedure to determine the weights is described in Methods.) That is, the interaction of a dipole with other dipoles is replaced by the weighted sum of the interactions of the dipole with cycles. It should be noted that a set of weights exactly satisfying eq. 3 can be obtained only for ice structures without polarization nor defects.
We calculate the cumulative interaction between water molecule i and cycles located within a distance r from the water molecule: , where r i (c j ) is the distance between the vertex i and directed cycle j, which is defined as the shortest one among the distances between the vertex i and any vertex in cycle j. The distance is defined to be zero when the vertex i is a member of cycle j. ψ ij is the sum of the Coulombic interactions between molecule i and all the dipole vectors composing directed cycle j. As seen in Figure 3, the cumulative interaction at r = 0 is close to that at long distances, indicating that the binding energy of a water molecule in ice is mostly determined by the interactions with the dipole vectors that share the same cycles with the central water molecule. As shown by the blue curves, the SD of the interaction energy decreases with increasing r when the molecule-molecule interactions are accumulated. In   Figure S4 of SI, the strong interaction of this type of cycle is equivalent regardless of the number of edges in a cycle. This fact guarantees that the interaction with these cycles is equivalent irrelevant to the arrangement of the surrounding hydrogen bonds. Figure S3 demonstrates that the sum of the interaction energies of a water molecule with cycles (a) and cycles (b) is almost the same as the whole potential energy of the water molecule and the contribution from cycles (c) is negligible.

Effects of Defects
Water molecules at an interface or defect cannot be a member of a directed cycle. The interaction of such a molecule with another molecule is a dipole-dipole interaction, which is long-ranged compared with the interaction between a water molecule and a directed cycle. As a result, the distribution of the accumulated interaction energy becomes wide. Figure 4 shows the cumulative interaction of water molecules for the ice slab in vacuum. The SD of the cumulative interaction remains large for water molecules near the ice/vapor interface. The convergent tendency of the cumulative interactions recovers as we move from the surface to the interior of ice.
7 arrangements. The SD of each of ten types of water molecules is close to that of ice Ih ( Figure S5 of SI). The SD is also larger for ice VI and VII with an interpenetrating double network. This is because there are no hydrogen bonds between the two networks to constrain the orientation of the molecules, and therefore strong attractive and repulsive forces are exerted depending on the orientation of the proximate molecules. This is consistent with the occurrence of site disorder in these ices. 7 8

Discussion
We find a hyperhomogeneity in hydrogen-disordered ice that reconciles the local heterogeneity called the site disorder with Pauling's model to evaluate the residual entropy by the small disparity in the binding energy. The distribution function of the interaction energy between a water molecule and nearby water molecules is fairly wide.
The corresponding distribution becomes substantially narrow when every water molecule is divided into two dipoles and the interactions of a molecule with other molecules are replaced by the interactions with directed cycles consisting of the dipoles. Such a resummation is allowed due to the ice rules: the hyperhomogeneity is an intrinsic property of ice.
There are efficient ways of accumulation to obtain the Madelung constant for an ionic crystal. If the Coulombic interactions are accumulated in an inappropriate order, the values will vary very widely depending on the order of summation, but the final value should be identical. 12 Similarly, we find that there is an efficient way of accumulation to obtain a well converged potential energy of a water molecule in hydrogen-disordered ice. We demonstrate that the molecule-cycle interaction is strong when the cycle includes the water molecule. The contribution from the other cycles is quite small and only acts as a noise. On the other hand, the SD of the cumulative interactions becomes very large in the medium distance and indicates a slow convergence if the interactions are accumulated molecule by molecule in order of proximity.
It is expected that there are large energetic frustrations associated with the energy inequality due to the absence of directed cycles at interfaces and near impurities. Such frustrations may trigger a solid-solid phase transition such as hydrogen-ordering in ice Ih that occurs only when specific electrolytes are doped. [13][14][15][16] The orientation of water molecules at the ice surface has been observed by advanced spectroscopic methods such as SFG (sum-frequency generation). 17 When interpreting experimental measurements, one may assume implicitly that the water molecules at the surface are energetically equivalent. However, unlike inside the bulk ice, the distribution of the binding energies of water molecules is wide at the surface. Due to this heterogeneity, the Pauling model can break down for the ice surface and there may be significant biases in molecular orientation. 18,19 The orientation of spins in magnetic substances called "spin ice" also obeys a rule that is topologically equivalent to the ice rules. 20,21 There may be a similar effect on the spin-spin interactions in spin ices. The residual entropy is observed for crystalline carbon monoxide but there is no high order regularity that can correspond to the ice rules.
( Figure S1 of SI) The hyperhomogeneity cannot exist for such a system.

Modeling of Ice
We prepare 30 hydrogen-disordered structures for each ice phase using the GenIce tool. 11 All the structures obey the ice rule and have zero net dipole moment. The number of molecules in a cell is ~1000, which is large enough for the present analyses. The TIP4P/Ice water model is employed. 5 The densities are set to 0.945, 1.1655, 1.2616, 1.389, and 1.688 for ices Ih, III, V, VI, and VII, respectively. Each structure is optimized at a fixed volume using the steepest descent energy minimization method. 22 We generate an ice Ic slab structure in vacuum to examine the effect of the interface. The dimensions of the ice crystal consisting of 1024 water molecules are 2.55x2.55x5.1 nm and the dimensions of the simulation cell are 2.55x2.55x10 nm. The structure is optimized using the steepest descent energy minimization method.
For long-ranged behaviors of SD appearing in Fig. 1b, 1,000 structures with N = 200,000 to 400,000 for each crystalline form are generated and are subject to the steepest descent minimization of the total energy.

Calculation of the weights
The weight w i for each cycle is calculated as follows. First, all 2N hydrogen bonds are labelled for an ice structure made of N molecules. We search for all directed cycles of size n or smaller and label them. The ij element of a nonsquare matrix A is set to 1 if directed cycle j contains hydrogen bond i, otherwise set to 0. If the number of directed cycles is sufficiently greater than 2N, a columnar vector w satisfying A·w = 1 can be determined by the pseudo inverse matrix method, where 1 is a 2N-row columnar vector in which all elements are 1. A sample python code to calculate the weights of the directed cycles in a given structure of ice is provided in the GitHub repository: https://github.com/vitroid/CoverByCycles.
If the size n of the cycle is infinitely large and the ice structure is free from defects and not polarized , it is easy to prove that there exists w that satisfies A·w = 1 . This is because it is always possible to draw a cyclic single stroke path that passes through all edges only once in the graph of such an ice structure, and the weight is unity. Another example is the case that the directed graph is tessellated into directed cycles where all the weights are unity (see Section S6 of SI).

Data Availability
In the interest of verifiable research, the data and code used to generate all Figures in both the main text and the supplementary information are available at https://github.com/vitroid/hyperhomogeneity. (Full disclosure will be made after acceptance of the manuscript.)