Quantum valence bond ice theory for proton-driven quantum spin-dipole liquids

We present a theory of a hybrid quantum liquid state, $\textit{quantum spin-dipole liquid}$ (QSDL), in a hydrogen-bonded electron system, by combining a quantum proton ice and Anderson's resonating valence bond spin liquid, motivated by the recent experimental discovery of a proton-driven QSDL in $\kappa$-H$_3$(Cat-EDT-TTF)$_2$ (a.k.a. H-Cat). In our theory, an electron spin liquid and a proton dipole liquid are realized simultaneously in the ground state called $\textit{quantum valence bond ice}$, while neither of them can be established independently of the other. Analytical and numerical calculations reveal that this state has a volume-law entanglement entropy between spins and dipoles, which is far beyond the (crude) Born-Oppenheimer approximation. We also examine the stability of QSDL with respect to perturbations and discuss implications for experiments in H-Cat and its deuterated analog D-Cat.

Introduction. -Hydrogen bonds occupy a unique position in condensed matter physics, where we can expect a huge quantum paraelectricity due to its lightness. Most hydrogen-bonded systems so far were treated within the (crude) Born-Oppenheimer approximation, where we ignore quantum entanglement between atomic motion and electronic states. However, recently it was recognized that in many hydrogen-bonded materials quantum nature of protons plays an essential role and even affects basic properties of the system [1][2][3][4][5][6]. Indeed, protons can tunnel quantum mechanically between two stable positions of the O-H· · · O bonding and furthermore the hydrogen bonds are subject to the ice rule due to the frustration in a wide range of systems, including ferroelectric KH 2 PO 4 [4,7], hexagonal water ice [5,8], and molecular crystals [6,9,10].
Recently, a new member of hydrogen-bonded material, purely organic κ-H 3 (Cat-EDT-TTF) 2 (called H-Cat), was synthesized [11][12][13][14][15][16][17][18][19]. H-Cat has a quasi-twodimensional (2D) triangular lattice structure of dimers with moderate interlayer coupling, and each layer is connected by hydrogen bonds which also form a triangular lattice. Experimentally, protons are unfrozen down to the lowest temperature and the observed quantum paraelectricity can be attributed to proton tunneling [16]. In sharp contrast to the above mentioned hydrogen-bonded materials, H-Cat is shown to be a dimer Mott insulator, where we can expect a strong electron correlation. Moreover, in H-Cat, spins delocalized over the π-conjugated orbital are also fluctuating without any magnetic order below a temperature scale close to that of the quantum proton motion. This implies that a quantum liquid state of unknown type is realized in H-Cat, where both protons and electrons are quantum mechanically disordered simultaneously and are strongly entangled with each other, which makes this material far beyond the Born-Oppenheimer approximation. On the other hand, the deuterated counterpart D-Cat exhibits a charge den-sity wave (CDW) state below 185 K, where protons and electrons together show a coupled charge order, indicating strong Coulomb interaction between them [11,13,16]. Such a large isotope effect strongly suggests that the proton tunneling is responsible for the hybrid quantum liquid state in H-Cat. However, it is unclear why the proton and spin fluctuations can coexist in the quantum liquid state.
The triangular lattice of H-Cat is a prototypical example of frustrated lattices and such a frustrated lattice can impose constraints on proton and spin configurations, leading to macroscopic degeneracy at a classical level. Quantum fluctuations could realize a quantum liquid state out of the degenerate classical states as originally proposed in Anderson's theory of the resonating valence bond (RVB) state [20][21][22]. The basic mechanism of realizing a quantum liquid was later employed in the theory of quantum spin/water ice where the RVB state is exactly realized by tuning a system into the exactly solvable Rohksar-Kivelson point [21][22][23][24][25][26]. This observation would imply that both the electron spin liquid and the proton dipole liquid can be treated in a unified manner to describe a spin-dipole coupled quantum liquid.
In this Letter, motivated by the experimental discovery of H-Cat, we present a theory of a hybrid quantum liquid state, quantum spin-dipole liquid (QSDL), in hydrogenbonded electron systems by combining the concepts of Anderson's RVB state of electron spins and a quantum ice of protons. Although a hybrid quantum spin liquid itself is discussed in the literature [27][28][29][30][31][32][33], including spincharge coupled systems without hydrogen bonds [34][35][36][37][38], its entangled nature has never been discussed in detail. In our theory, the RVB state is extended to describe the entanglement between spins and dipoles. We give a simple concrete description of the QSDL based on a minimal model which incorporates three key features of the QSDL, i.e. macroscopic classical degeneracy, quantum proton tunneling, and electron correlations. The result-ing QSDL is named a quantum valence bond ice, and the property of this state will be discussed in detail from analytical and numerical perspectives.
Model. -We consider an idealized model for the QSDL defined by the following Hamiltonian for holes and protons, The holes are described by the t-J model where double occupancy is prohibited by the strong onsite interaction [39]: a real t el is a hole hopping parameter inside a tetrahedron denoted by , and J el > 0 is a spin-spin exchange interaction between nearest-neighbor (NN) holes [see Fig. 1(a)]. c js is an annihilation operator of a hole at the jth site with a spin s, and S j are spin-1/2 operators defined for this hole. The summation over s = ↑, ↓ is implied. σ ij are Pauli matrices defined on each hydrogen bond ij connecting the NN tetrahedra, and two eigenstates with eigenvalues ±1 of the σ z ij operator correspond to the stable positions of a proton on the hydrogen bond ij [16,19]. The tetrahedra sit on vertices of a diamond lattice and σ z ij is defined such that σ z ij = +1 for the sites i(j) belonging to the A(B)-sublattice of the diamond lattice corresponds to a proton position away from the A-sublattice site i. A real t pro represents quantum tunneling between the two positions, and J pro > 0 is a Coulomb repulsion between the protons on an NN pair of bonds ij , kl . g > 0 is a Coulomb repulsion between the hole and the proton. We will later extend the model Eq. (1) and discuss effects of an additional proton-proton interaction and an intertetrahedron hole hopping [see Fig. 1 We consider holes to be quarter-filled for each tetrahedron referring to the dimer Mott insulator H-Cat [11,13,16]. Essentially the same argument can also be applied for an electron quarter-filled system. While H el has a form of the t-J model, H pro is a transverse-field Ising model on the frustrated pyrochlore lattice, and related models have been investigated as a model for quantum spin liquids [21][22][23][24][25][26]. Our theory is based on the frustration of the protons in an ice manifold rather than that of the spins. It should be noted that g and J pro are supposed to be very large, g, J pro ∼ O(10 3-4 ) K because they come from a strong Coulomb repulsion. In addition, we suppose |t pro | J pro and |t el | g by referring to the parameters in H-Cat where t pro ∼ 1-10 K and t el ∼ 10-100 K [11,16,18], which allows a perturbation theory with respect to t pro , as well as t el . The magnetic coupling is assumed to be J el ∼ 10-100 K taken from the value for H-Cat, and one can focus only on the singlet sector of a tetrahedron when the temperature is well below J el . Classical limit. -In the case of t pro = 0, the Hamiltonian is exactly solvable and the ground states are degenerate for all proton configurations obeying a (2 in 2 out) ice rule as shown in Fig. 1(a). In the limit of t el /g → 0, the proton configuration completely couples to the valence bonds of hole pairs, and the position of the valence bond is fixed inside each tetrahedron in the ground states. Since the valence bonds also form an ice configuration, this state may be called a classical valence bond ice. This basically holds true even for a nonzero t el where one can diagonalize the electron Hamiltonian of each tetrahedron in the valence bond basis and the ground state of a tetrahedron will be a dressed valence bond state with a finite excitation gap ∼ g. Thus, every degenerate ground state of the total system at t el = 0 is adiabatically connected from that at t el = 0.
Quantum effects.
-We now introduce a quantum tunneling of protons t pro = 0, which will lift the macroscopic degeneracy of the classical ground states in the same way as in the quantum spin ice [21][22][23][24][25][26]. We can treat the effects of t pro perturbatively when |t pro |/J pro 1 by using the degenerate perturbation theory. The condition |t el |/g 1 makes a physical picture clear, though essentially the same argument is true for a moderate t el . Then, if a proton moves from its ground state position breaking the ice rule, a hole in the nearby tetrahedron will be pushed away from the proton because of the repulsive interaction g between them. In this way, the proton motion is always accompanied by the hole hopping inside a tetrahedron and vice versa. We note that such a charge-correlated hopping leads to an effective coupling between spins and protons despite the absence of an explicit spin-orbit interaction because holes form valence bond pairs due to the magnetic interaction J el . A nontrivial contribution appears from the 6th-order perturbation in t pro /J pro . This process results in the following effective Hamiltonian in addition to a constant energy shift of the degenerate ground states, where the flipping operator is given by the product of hole and proton operators along a hexagonal plaquette loop as shown in Fig. 2: ,12 c 12 with σ ± = (σ x ± iσ y )/2 in the lowest order approximation for t el /g. The sum runs over all hexagonal plaquettes of the diamond lattice, and the coefficient K > 0 is estimated to be O(t 6 pro /J 5 pro ) with a prefactor O(t 6 el /g 6 ) when |t el | is small [40].
The effective Hamiltonian Eq. (5) is formally equivalent to the quantum ice model, which has been studied extensively in the context of the quantum spin ice [21][22][23][24][25][26]. We stress, however, that the physical constituents of the low energy degrees of freedom | , | are very different from those in the quantum spin ice; they are spindipole coupled degrees of freedom in the former, while solely spins in the latter. The ground state of the effective Hamiltonian is known to be the quantum liquid state, which is described by an emergent compact U(1) electromagnetic field which is deconfined in 3 + 1 dimensions. Thus, we have theoretically established that the quantum proton motion can indeed drive the system into a QSDL, namely the quantum valence bond ice, and form a stable three-dimensional (3D) phase. The coupled spins and dipoles behave as gapless "photons" of the emergent U(1) gauge theory, while local proton excitations which break the ice rule with an energy gap ∼ J pro are regarded as "monopoles" [21][22][23][24][25][26]. Entanglement entropy of QSDL. -One of the most characteristic features of the QSDL is a large entanglement between protons and holes.
The corresponding entanglement entropy defined by S EE = −Tr el ρ el log ρ el with ρ el = Tr pro |Ψ Ψ| for the ground state |Ψ can provide useful information of its entanglement.
Now one can estimate that S EE ∼ log(the number of flippable ice states) = O(N tot ), where N tot = N pro + N el and N pro/el is the total number of proton bonds/electron sites. This is easily understood for example based on the explicit calculation of S EE for the Rohksar-Kivelson state [40,41]. Therefore, the QSDL has an entanglement entropy with a volume law in contrast to the area law for a real space bipartition [42,43]. Such a large entanglement between protons and holes is an intrinsic nature of the QSDL, which is far beyond the Born-Oppenheimer approximation.
Stability of QSDL. -Up to here, we have been based on the idealized model Eq. (1) where a minimum number of terms are included to describe the QSDL. However, there could be additional terms in general, and considerations on such terms are important to discuss its relevance in real materials. Here, we discuss stability of the QSDL by introducing two types of perturbations into the model Eq. (1); (i) an additional interaction between protons which favors a classical order by lifting the macroscopic degeneracy of protons, and (ii) an intertetrahedron hole hopping which may also lift the macroscopic degeneracy of electrons. The following discussions will be supplemented with numerical calculations later.
First, we examine effects of additional proton interactions such as J pro σ z ij σ z kl shown in Fig. 1(b), and discuss implications for the qualitative difference in X-Cat (X = H, D). The additional interaction lifts the macroscopic degeneracy of the classical ground states, and an infinitesimal J pro will lead to an ordered state when t pro = 0 by the order-by-disorder mechanism. Indeed, classical (t pro = 0) spin ice models with additional interactions have been extensively studied, where sev-eral ordered states are stabilized depending on types of the additional interactions [44,45]. In the present system, we consider J pro interaction which favors an antiferroelectric (AFE) state of protons, but it competes with the quantum tunneling t pro [9,22,46]. Based on these observations, one can draw a schematic phase diagram of the perturbed model with both J pro and t pro as shown in Fig. 3(a), where the QSDL with a volume-law entanglement entropy remains stable for J pro K, while the AFE state with a small entanglement is favored for J pro K. Note that because g |t el | holes simultaneously show a valence bond solid (VBS) state corresponding to the proton AFE state realized when J pro K. Experimentally, an H bond will have a larger t pro and hence a larger K than those of a D bond. Therefore, an H-bonded system can exhibit a QSDL, while a D-bonded system shows a classically ordered AFE/VBS state. This is a generic feature of proton-driven quantum liquids and therefore may be relevant to the experimentally observed isotope effect in X-Cat, where H-Cat is a QSDL while D-Cat exhibits an AFE CDW order [11,13,16].
Next, we consider the intertetrahedron hole hopping along the hydrogen bonds, t el ij c † iσ c jσ as shown in Fig. 1(b), within a perturbative regime |t el | |t el |. The second order perturbation process with respect to t el will generate two kinds of additional interactions between holes on neighboring tetrahedra when beginning with a Hubbard-like model from which the t-J model Eq. (2) is derived [40]. The first term is a Heisenberg term J el [S i · S j − n i n j /4], and the second one is a density-density interaction J el [n i (n j − 1) + (n i − 1)n j ] with 0 < J el , J el J el for the low-energy manifold. Although the J el S i · S j term introduces a band dispersion for magnetic excitations, it is negligible at low temperature T J el because the magnetic band width is at most ∼ J el , which is much smaller than the singlettriplet energy gap ∼ J el . Effects of the density contribution −J el n i n j /4 will be strongly suppressed by the large interaction g J el between holes and protons because a VBS configuration favored by this term has an energy cost ∼ g. Similarly, the other density interaction J el is also negligible [40]. Therefore, the QSDL is stable with respect to the perturbations and can be realized in more realistic models.
Numerical simulations. -The above discussions on the stability can be supplemented by numerical calculations using exact diagonalization for a finite-size system with perturbations, such as the J pro term. To this end, we consider a 2D version of our model Eq. (1) which consists of 8 proton bonds and 16 electron sites (8 holes) with a periodic boundary condition [40]. Although the system size is very small, this model can reproduce much of the main results and hence supports the above discussions. When J pro is moderately strong, the 2D model shows an AFE order of protons, and holes simultaneously form a VBS state for a large g. This state will have a corresponding long range order in the thermodynamic limit. On the other hand, for a small J pro it can exhibit a QSDL-like state which is a superposition of many flippable ice states. Indeed, we find that out of the total N = 18 ice states there are N flip = 16 flippable ice states, which are connected by the lowest order perturbation in t pro . The QSDL-like state is essentially given by a superposition of these 16 reduced ice states. The AFE/VBS state is clearly distinguished from the QSDL-like state even in the present finite-size system, which is well-captured by the entanglement entropy S EE (t pro , J pro ) between protons and holes as shown in Fig. 3(b). In the limit of t pro → 0, the ground state is given by a superposition of simple AFE/VBS states with fourfold degeneracy due to the translational and rotational symmetries. Therefore, S EE (t pro → 0) → log 4 for such a state, while S EE → log N flip = log 16 deeply inside the QSDL region with a moderately large t pro . This qualitatively reproduces the phase diagram in Fig. 3(a) proposed to describe the difference between H-Cat and D-Cat.
In the 3D thermodynamic system, N flip , N ∼ e Ntot and the present QSDL-like state will be replaced by a true quantum liquid state where a macroscopic number of ice states are involved, resulting in a volume-law entanglement entropy, S EE = O(N tot ). Note that similar crossover behaviors in the finite size system can also be seen in correlation functions [40], and those quantities will characterize the corresponding phase transition in the thermodynamic limit. We can also discuss the effects of an intertetrahedron hopping t el , and numerical results suggest that induced intertetrahedron valence bond interactions have only small effects on QSDL if |t el | |t el |, and QSDL is parametrically stable [40]. Detailed studies on the stability of QSDL would require more elaborate investigations, such as cluster-based calculations [47], and are left for future work.
Summary and discussions. -We have developed a theory of a hybrid quantum liquid, QSDL, by combining the quantum proton ice and Anderson's RVB state, motivated by the recent experimental discovery of the QSDL in H-Cat. We proposed an idealized model and demonstrated that proton tunneling drives the system into the quantum valence bond ice, a QSDL described by an emergent U(1) gauge theory. Our theory sheds light on the essential roles of protons in the realization of QSDLs proposing a coherent understanding of spindipole coupled systems, and thus can provide a basis for future developments not only for X-Cat but also for other QSDL candidate systems [31,48]. For example, a hydrogen-assisted Kitaev spin liquid material H 3 LiIr 2 O 6 was recently discovered [31], and proton tunneling is suggested to be important to stabilize the quantum spin liquid [49][50][51][52].
The large entanglement predicted by our theory characterizes a QSDL in general, and it would lead to spin-dipole coupled response such as a magnetoelectric effect, which could provide an experimental signature of entangled spins and protons. Other experimentally observable quantities characterizing a QSDL are dynamical ones, which can be captured by neutron scattering or nuclear magnetic resonance. Further investigation is necessary to reveal a smoking gun signature for the strong entanglement between the two degrees of freedom. Here we derive the effective low-energy Hamiltonian Eq. (5) in the main text from the original Hamiltonian Eq. (1), based on the degenerate perturbation theory with respect to t pro [1]. The 6th-order perturbation in t pro gives where we decompose the original Hamiltonian H(t pro ) given by Eq. (1) in the main text into H(t pro ) = H 0 + V t , and H 0 = H(t pro = 0). P is a projection operator onto the degenerate ice state manifold with an energy E 0 .
Though an intratetrahedron hole hopping t el has been included in H 0 and the perturbation itself is done for a small t pro /J pro , we concentrate on a principal contribution coming from the lowest order in t el /g. To do so, we take a basis set consisting of product states of the hole and proton sectors {|pro ⊗ |el }. Each hole state |el can be explicitly obtained by diagonalizing the hole Hamiltonian of each tetrahedron for a given proton configuration |pro . By inserting an identity operator 1 1 = |pro ⊗ |el pro| ⊗ el| , we find for two flippable ice states |a and |b , where h ab = a|H eff |b , J = 4J pro , and G = 2g. |el a(b) is the hole wavefunction in the ice state a(b). The overlap matrix element can be calculated numerically and is O(t 6 el /g 6 ) for a small |t el | as shown in Fig. S1. The overlap amplitude takes the same value for any pair of electron states connected by the 6th-order perturbation with respect to t pro , and therefore we simply denote K ≡ −h ab to obtain Eq. (5) in the main text as H eff = −K |a b| . The effective Hamiltonian can explicitly be expressed in terms of hole operators c ( †) is , and |el a el b | = plaquette c † is c js (i, j ∈ ) is nothing but a ring-hopping operator for each hexagonal plaquette loop in the leading order of the perturbation in t el /g.

B. Entanglement entropy for the Rohksar-Kivelson state
The entanglement entropy between protons and electrons can be calculated explicitly for an extreme limit of a quantum spin-dipole liquid (QSDL), the Rohksar-Kivelson (RK) state. The RK state is the exact ground state of the ice model Eq. (5) in the main text with an additional potential µ [| | + | |] fine-tuned to µ = K. The RK state is a superposition of all possible ice states, |RK = (1/ √ N ) C |pro C ⊗ |el C , where C represents an ice configuration and N is the dimension of the ice manifold. |pro C and |el C are proton and electron wavefunctions for a fixed ice configuration C. Eventually, the entanglement entropy between protons and holes is given by S EE = log N and it shows a volume law where t el is an intertetrahedron hopping along a hydrogen bond ij . U > 0 is an onsite interaction and V (< U ) is an intersite interaction inside each tetrahedron.
Assuming U and V to be large, the 2nd-order perturbation in H t gives several terms, J el [S i · S j − n i n j /4] with J el = 2t 2 el /(U − V ), J el [S i · S j − n i n j /4] with J el = 2t 2 el /U , and J el [n i (n j − 1) + (n i − 1)n j ] with J el = t 2 el /V . The first term is an exchange interaction inside each tetrahedron, which is included in Eq. (1) in the main text, while the other two terms are new intertetrahedron interactions. Sufficiently large U |t el | and V |t el | also lead to constraints on possible hole filling so that n i ≤ 1 for all sites i, and i∈ n i = 2 for each tetrahedron assuming quarter filling of holes. The resulting state can be considered as a "tetramer Mott insulator." The perturbed Hamiltonian is now reduced to In the following, we only consider a singlet-pair state of holes on each tetrahedron since there is an energy gap ∼ J el J el , J el between singlet and triplet states. We note that the −(J el /4) i,j∈ n i n j term does not matter in the present system where the total fermion number inside a tetrahedron is fixed.
The J el term is nonzero only when both of the two sites connected by a hydrogen bond ij are occupied, and similarly the J el term is nonzero only when one of the two sites connected by the hydrogen bond is occupied. Effects of the J el term will be strongly suppressed by a large interaction g J el between holes and protons included in the original model because of the energy cost g for a configuration favored by J el . On the other hand, the J el term effectively adds to the interaction g because the J el term also favors the same local charge configurations as the g term does. Therefore, these two potential interactions between singlet pairs on the nearest-neighbor (NN) tetrahedra could effectively be absorbed into the interaction g, or equivalently the intratetrahedron hopping t el . Then, it is concluded that thanks to the stability of the QSDL for a wide range of t el /g the QSDL is also stable with respect to the intertetrahedron hopping t el . This behavior is numerically confirmed by exact diagonalization in the next section.

D. Exact diagonalization
The Hamiltonian used in the numerical exact diagonalization is written as where H is a two-dimensional (2D) version of Eq. (1) in the main text, in which tetrahedra are located on vertices of the square lattice as shown in Fig. S2. H pert describes the perturbation: The J pro term represents an additional interaction between protons with a properly chosen η ij kl = ±1 which favors a certain antiferroelectric (AFE) state, and the J el , J el terms are intertetrahedron interactions between holes taken from Eq. (S6). Note that for a large g holes simultaneously exhibit the valence bond solid (VBS) state out of the classical valence bond ice manifold corresponding to the proton AFE state. We consider a finite-size system shown in Fig. S2, where there are 8 proton bonds and 16 electron sites (8 holes) with a periodic boundary condition.
In this system, there are 18 ice states when t pro = J pro = J el = J el = 0, and 16 are connected by the 2ndorder perturbation in t pro . The 2nd-order tunneling is possible because the linear system size for one direction is 2 and the periodic boundary condition has been imposed. Furthermore, the conservation of the total spin for each tetrahedron, [H, S ] = 0 for S = i∈ S i , greatly reduces the computational cost for the unperturbed Hamiltonian H. Taking this advantage, we focus only on the subspace where S = 0 for every tetrahedron also for the perturbed Hamiltonian H tot because there is an energy gap ∼ J el between the S ≡ 0 sector and the other sectors in H. From now on, we omit J el because the triplet sectors of the original Hilbert space are all projected out.
We calculate an entanglement entropy between protons and holes: where {λ n } are the squared singular values of the Schmidt decomposition of the ground state wavefunction |Ψ into the proton and electron subspaces, The results are shown in Fig. 3 in the main text, and the two ground states, AFE/VBS and QSDL states, are wellcharacterized by S EE . The saturated value S EE log 16 for a large t pro reflects the fact that the ground state is essentially given by the superposition of 16 flippable ice states with a nearly equal weight, and the ground state can be considered as a QSDL in this sense. Note that the QSDL is a stable phase extended in the parameter space, although it is parametrically small as expected from other quantum ice systems [2,3]. A similar characterization is also possible by correlation functions corresponding to the long-range order of the AFE/VBS state. We consider the following two order parameters, where η ij = ±1 characterizes the proton AFE configuration and |VBS is the VBS configuration of the holes on a tetrahedron corresponding to the proton AFE shown in Fig. S2(b). The calculation results are shown in Fig. S3. Note that the AFE correlation function can be evaluated explicitly for an extreme AFE state of protons, being 4 degenerate AFE configurations, and it is given by AFE| O 2 AFE |AFE = 1/2. Similarly, for an extreme QSDL state |QSDL = (1/ N flip ) C∈flip |pro C ⊗ |el C , which is the equal-weight superposition of N flip = 16 flippable ice states, the AFE correlation is given by QSDL| O 2 AFE |QSDL = 1/4. These behaviors are reproduced in numerical calculations. Note that QSDL| O 2 AFE |QSDL → 0 in the thermodynamic limit. The same thing happens for O 2 VBS . Finally, we discuss the effect of the intertetrahedron interactions J el , J el arising from the hole hopping |t el | |t el |. Fig. S4 shows entanglement entropies S EE for different values of (g, J el ) at J el = 0 and (g, J el ) at J el = 0. The other parameters used are t pro = 0.2J pro and J pro = 0. S EE (g, J el ) is decreased only slightly by the introduction of J el because it competes with g, but is negligible for a large g. S EE (g, J el ) is enhanced a little by J el , as it favors valence bond configurations including the flippable ice states. These behaviors are consistent with the qualitative discussion in the previous section. Therefore, we conclude that the QSDL is parametrically stable for a small intertetrahedron hole hopping t el .