Polarization Plateaus in Hexagonal Water Ice $I_h$

TL;DR

This paper investigates how external electric fields induce polarization plateaus and phase transitions in hexagonal water ice I$_h$, revealing complex degeneracy reduction and algebraic correlations through analytical and numerical methods.

Contribution

It introduces a minimal model to study polarization plateaus and order-disorder transitions in water ice I$_h$ under external fields, connecting ice rules to dimer models and Kasteleyn theory.

Findings

Identification of polarization plateaus at specific field directions.

Observation of an order-disorder transition of Kasteleyn type.

Protons exhibit algebraically decaying correlations and characteristic structure factors.

Abstract

The protons in water ice are subject to so called \emph{ice rules} resulting in an extensive ground state degeneracy. We study how an external electric field reduces this ground state degeneracy in hexagonal water ice I$_h$ within a minimal model. We observe polarization plateaus when the field is aligned along the $[001]$ and $[010]$ directions. In each case, one plateau occurs at intermediate polarization with reduced but still extensive degeneracy. The remaining ground states can be mapped to dimer models on the honeycomb and the square lattice, respectively. Upon tilting the external field, we observe an order-disorder transition of Kasteleyn type into a plateau at saturated polarization and vanishing entropy. This transition is investigated analytically using the Kasteleyn matrix and numerically using a modified directed-loop Monte Carlo simulation. The protons in both cases exhibit algebraically decaying correlations. Moreover, the features of the static structure factor are discussed.

Tables

Table 1. Table 1: Location of scattering features in reciprocal space between different conventions.

scattering feature	$𝒌$ in [Tang et al., 2011]	${𝒌}_{Ih,or}$
$S_{xx}$ log. peak	$(\pi ,0)$	$(\pi ,0,\pi )$
	$(\pi ,2\pi )$	$(3\pi ,0,-\pi )$
$S_{xx}$ pinch point	$(\pi ,\pi )$	$(0,0,2\pi )$

Equations

$${\mathbf{d}}_i=s_{i}d{\hat{\mathbf{e}}}_{\mathbf{i}}.$$

$$H=J\sum _{\alpha }{L}_{\alpha }^2-\sum _i{s}_{i}d\mathbf{E}\cdot {\hat{\mathbf{e}}}_{\mathbf{i}},$$

$$p(i\to j)=\frac{z(i\to j)}{{\sum }_{k}z(i\to k)},$$

$${ϵ}_{\alpha }=-\frac{1}{2}\mathbf{E}\sum _{i\in {⊠}_{\alpha }}{s}_i{\hat{\mathbf{e}}}_{\mathbf{i}},$$

$${ϵ}_{\alpha }=-s_{\alpha ,\text{int}}\mathbf{E}\cdot {\hat{\mathbf{e}}}_{\mathbf{\alpha },\text{int}}-\frac{1}{2}\sum _{\text{ext}i\in D_{\alpha }}{s}_i\mathbf{E}\cdot {\hat{\mathbf{e}}}_{\mathbf{i}},$$

$$\mathbf{E}=(1,0,100)$$

$$\mathbf{E}=\left(\frac{1}{\sqrt{2}},100,\frac{1}{\sqrt{2}}\right)\approx (0.71,100,0.71).$$

$$a_c=1+\frac{1}{a_c}⟹a_c=\frac{1}{2}\left(1+\sqrt{5}\right)\approx 1.618.$$

$$T_K=\frac{-{\sum }_i{s}_i{\hat{\mathbf{e}}}_{\mathbf{i}}\cdot {\mathbf{E}}_{\perp }}{\log a_c}=\frac{{ϵ}_a}{\log a_c}=3.396,$$

$$\sinh \left(-\frac{{ϵ}_a}{k_B{T}_K}\right)=\cosh \left(-\frac{{ϵ}_b}{k_B{T}_K}\right).$$

$$\frac{T_K}{E_{\perp }}\approx 3.385$$

$$\mathcal{Z}=\sum _{\{n_1,n_2,\dots ,n_N\}}g(n_1,n_2,\cdots ,n_N)z_1^{n_1}{z}_2^{n_2}\cdots z_N^{n_N},$$

$$\mathcal{Z}=|\mathrm{Pf}K|$$

$$K=\left(\begin{array}{cc}0& \tilde{K}\\ -{\tilde{K}}^T& 0\end{array}\right).$$

$$\mathcal{Z}=|\det \tilde{K}|.$$

$$\tilde{K}=b+\frac{zw}{b}-aw+\frac{z}{a}=\det \tilde{K},$$

$$\begin{gathered} \frac{1}{N}\log \mathcal{Z}= \\ \frac{1}{(2\pi i{)}^2}\underset{|w|=1}{\oint }\frac{dw}{w}\underset{|z|=1}{\oint }\frac{dz}{z}\log \left(b+\frac{zw}{b}-aw+\frac{z}{a}\right). \end{gathered}$$

$${\rho }_a=\frac{N_a-N_{a^{-1}}}{N}=a\frac{\partial }{\partial a}\log \mathcal{Z}.$$

$${\rho }_a=\frac{1}{(2\pi i{)}^2}\underset{|w|=1}{\oint }\frac{dw}{w}\underset{|z|=1}{\oint }\frac{dz}{z}\frac{-aw-\frac{z}{a}}{b+\frac{zw}{b}-aw+\frac{z}{a}}.$$

$$z_0=0z_1=ab\frac{aw-b}{aw+b},$$

$${\rho }_a=\{\begin{array}{ll}1& \text{if}(a>b+\frac{1}{b}+\frac{1}{a})\\ 0& \text{if}(b>\frac{1}{b}+a+\frac{1}{a})\text{or}(\frac{1}{b}>b+a+\frac{1}{a})\\ -1& \text{if}(\frac{1}{a}>b+\frac{1}{b}+a).\end{array}$$

$$\begin{gathered} {\rho }_a=\frac{1}{2\pi i}\underset{|w|=1}{\oint }\frac{dw}{w}\frac{aw}{b-aw} \\ +\frac{1}{2\pi i}\underset{\begin{array}{c}{w}_1\\ |w|=1\end{array}}{\overset{w_2}{\int }}\frac{dw}{w}\frac{-aw}{aw-b}-\frac{b}{aw+b}. \end{gathered}$$

$${\rho }_a=\Theta \left(1-\frac{b}{a}\right)-\frac{1}{\pi }{\varphi }_0+\frac{1}{2\pi i}\log \frac{b^2-a^2+2iab\sin {\varphi }_0}{b^2-a^2-2iab\sin {\varphi }_0}$$

$${\varphi }_0=\arccos \left(\frac{\left(a^2+b^2\right)\left(a^2{b}^2-1\right)}{2ab\left(a^2{b}^2+1\right)}\right),$$

$$S_H(\mathbf{k})\propto {\left|\sum _{j\in \text{H sites}}{n}_j{e}^{i\mathbf{k}\cdot {\mathbf{R}}_j}\right|}^2,$$

$$S_H(\mathbf{k})\propto {\left|\sum _j{s}_j\sin \left({\mathbf{a}}_{\gamma (j)}\cdot \mathbf{k}\right)e^{i\mathbf{k}\cdot {\mathbf{r}}_j}\right|}^2,$$

$$S_I(\mathbf{k})\propto {\left|\sum _j{s}_j{e}^{i\mathbf{k}\cdot {\mathbf{r}}_j}\right|}^2.$$

Full text

I Introduction

Water exists in a huge variety of different phases including liquid, gas, as well as numerous crystalline and amorphous phases. The most common crystalline phase of water on earth is the hexagonal phase ice Ih, cf. Fig.1(a). Bernal and Fowler proposed that ice Ih is a molecular solid formed by H2O moleculesBernal and Fowler (1933) which obey the ice rules: a single proton is placed on each of the four oxygen-oxygen link per oxygen site. There are two protons bonded covalently to each oxygen ion forming H2O molecules and the two other protons form hydrogen bonds. As such the ice rules impose local constraints that can be expressed as a zero divergence condition on a lattice and give rise to emergent gauge fields.Huse et al. (2003); Hermele et al. (2004) These local constraints prevent the H2O molecules from ordering and leave a massively degenerate proton subsystem with a finite residual entropy.Pauling (1935); Giauque and Stout (1936) Consequently, ice Ih is considered to be a prototype of geometrical frustration.

Because of the third law of thermodynamics, the equilibrium residual entropy is expected to vanish as temperature approaches zero, however, the nature of the true ground state of ice Ih is not yet settled.Minagawa (1981, 1990); Abe and Tokunaga (1993); Kirov (1996); Jackson and Whitworth (1997); Iedema et al. (1998); Lekner (1998); Rick (2005); Parkkinen et al. (2014) On one hand, a transition with a significant reduction in entropy is observed in KOH doped ice Ih at $T=72K$. Kawada (1972); Tajima et al. (1982); Matsuo et al. (1986); Suga (1997) On the other hand relaxation times grow beyond reasonable observation times for pure ice below $T\approx 110K$,Wooldridge and Devlin (1988); Suga (1997) obscuring any thermodynamic transition at lower $T$. Yet, experimental signatures of correlated proton hopping processes were reportedBove et al. (2009); Yen and Gao (2015) enabling the protons to possibly form a quantum proton liquid,Castro Neto et al. (2006); Benton et al. (2016) which would be a strong analogy to quantum spin liquids.Anderson (1973)

A magnetic equivalent of the ice rules exists in spin ice materials.Harris et al. (1997); Ramirez et al. (1999); Bramwell and Gingras (2001) Due to strong crystal fields, local magnetic moments are aligned parallel or antiparallel to links in and out of the tetrahedrons of a pyrochlore lattice. The local moments have to satisfy a modified ice rule: two spins point into a tetrahedron and two point out. Analog to ice, the constraint is equivalent to the zero divergence condition of classical electrodynamics, Gauss’ law, on a lattice.

In both cases, water ice or spin ice, the local constraints manifest themselves in certain features of scattering experiments. The ice rules cause an angle dependence of the scattering intensity around lattice Bragg points resulting in a bow-tie-like structure.Youngblood et al. (1980); Li et al. (1994) Such scattering maps contain information on proton-proton correlations, which have been computed for ice Ih using a large-$N$ theoryIsakov et al. (2015) or a projection operator approach.Benton et al. (2016)

In this article, we study the reduction of the residual entropy by applying an external electric field on the extensively degenerate ground state of the proton configurations. Generically, the degeneracy is expected to decrease upon applying an external field, but details depend on the direction of the field and which lattice symmetries the field breaks. The degeneracy remains extensive in certain field directions, e.g., when the field confines the remaining degrees of freedom to a stack of decoupled planar layers such as for pyrochlore spin ice in a $[111]$ magnetic field.Matsuhira et al. (2002); Moessner and Sondhi (2003) A small tilt of the field results in a transition into a long range ordered state.Moessner and Sondhi (2003).

We consider an external electric field that couples to the proton subsystem via a Stark term. The proton configurations are described by a minimal model of ice Ih satisfying the ice-rules: a proton can only take one of two positions on a link between two neighboring oxygen atoms and only a single proton per link exists. This allows to describe the proton configurations as Ising-like spins aligned parallel to the links.

We focus on fields along two easy axes causing a partial ordering of the protons. The remaining degrees of freedom reside on decoupled layers mappable to hardcore dimer models either on an honeycomb or a square lattice, respectively. The former is known from pyrochlore spin-ice,Moessner and Sondhi (2003) whereas the latter has no analog in spin-ice and, thus, is the main focus of this work. The field effectively reduces the dimensionality of the system to uncorrelated two-dimensional layers. While the ground state degeneracy is reduced, it remains extensive. A slight tilt of the field lifts the remaining degeneracy within each layer and leads to an order-disorder transition first described by Kasteleyn.Kasteleyn (1963) Using the Kasteleyn-matrix approachKasteleyn (1963); Cohn et al. (2001); Kenyon et al. (2006); Kenyon (2009) for the emergent planar dimer models, we investigate the critical behaviour in the thermodynamic limit and compare to numerical simulations using a modified directed-worm Monte-Carlo simulation on the three-dimensional model. Moreover, the static structure factor is obtained from the Monte-Carlo simulation. It shows typical bow tie or pinch point features,Lacroix et al. (2011) that get supplemented by logarithmically diverging peaks within the emergent dimer regime.

The article is structured as follows: In section II, the structure of ice Ih and the model are introduced. An explanation of the modified directed-worm Monte-Carlo method follows in Sec. III. Results based on thermodynamic observables are presented in Sec. IV following a discussion of the emergent dimers on a hexagonal lattice in Sec. IV.1 and the square lattice in Sec. IV.2. The latter is extended by utilizing the Kasteleyn-matrix to derive the dimer densities, cf. Sec. V. Section VI contains correlations and diffuse scattering maps.

II Model

Ice Ih has a bipartite hexagonal lattice with four oxygen atoms and eight protons per primitive cell (Fig. 1). The oxygen atoms form a four-fold coordinated tetrahedral structure of class P63/mmc. The following symmetries exist: a sixfold rotation axis, a twofold rotation axis perpendicular to the z-axis and a horizontal reflection plane. The structure can be considered as a stacking of planes in an ABAB pattern, where B is the reflection of A with respect to a horizontal plane.

Within the unit-cell, eight distinct oxygen-oxygen links exist. Each proton has two possible positions arranged symmetrically around the midpoint of a link. We neglect the displacement of the protons away from the link. If we assume a single proton per link, which is one of the ice rules,Bernal and Fowler (1933) the proton configuration can be represented as a configuration of electric dipoles $\bm{d}_{i}$

[TABLE]

The dipoles can be represented as effective Ising spins $s_{i}$ with respect to a local z-axis given by the direction $\bm{\hat{e}_{i}}$ of the oxygen-oxygen links.

We consider the following minimal model in terms of electric dipoles or Ising spins, respectively:

[TABLE]

where $L_{\alpha}=\sum_{i\in\boxtimes_{\alpha}}s_{i}$ is the sum over effective Ising spins on a tetrahedron $\boxtimes_{\alpha}$. In choosing the limit $J\rightarrow\infty$, we restrict the dipole configuration to the ice-rule manifold. The second term represents the Stark-coupling of the moment to an external electric field, $\bm{E}$. The direction of the unit vectors $\bm{\hat{e}_{i}}$ along the links is chosen such that a positive $s_{i}$ points from a site on sublattice A to a site on sublattice B. In the remainder, the constant $d$ is set to unity.

We apply the field along a) the c-axis $[001]$ and b) $[010]$ with respect to the orthorhombic unit cell, see Fig. 1 for an illustration. These field directions still leave an extensively degenerate ground state with finite entropy density. In the following, we investigate the case where the field is slightly tilted towards a) [100] and b) [101], such that a single lowest energy state exist. For the remainder of this article, the notation $\bm{E}=\bm{E}_{0}+\bm{E}_{\perp}$ is used for the external field, where $\bm{E}_{0}$ is the component along the easy axis and $\bm{E}_{\perp}$ denotes the tilt component with $\bm{E}_{\perp}\ll\bm{E}_{0}$.

III Modified Directed-Worm Monte Carlo Sampling

The spin- or proton system is simulated numerically using a directed loop Monte-Carlo method.Barkema and Newman (1998) Since a single spin flip would violate the second ice rule, loops of consecutive spins are created and flipped instead. In doing so, the simulation remains in the defect-free manifold. In order to keep the efficiency high, the sampling enters already at the loop creation: two defects are created by choosing a single spin randomly with uniform probability and flipping it. One of the defects executes a weighted random walk with weights depending on the energetics of the original local configuration and the possible new local configurations. As soon as that defect meets the other defect, the loop is closed and all corresponding spins will be flipped. This is also known as the long loop algorithm, which has the property of spanning on average a fraction of the whole system.Jaubert et al. (2011) A loop which winds once or multiple times around the system connects different winding sectors.

In order to satisfy detailed balance, backtracking of a defect must also be taken into account. If the defect tracks all the way back to its initial creation, no update occurs and a new initial spin randomly chosen.

In most cases, only the change of the local configuration on a single tetrahedron is taken into account when computing the probabilities to move the defect on site. However, this choice makes the algorithm highly inefficient for the field along $[010]$. Consider two neighboring vertices connected by a spin perpendicular to the field: the defect steps forward at the first vertex with high probability as the change in configuration leads to a gain in energy. At the next vertex, the defect has only a small probability to move forward as it would cost nearly the same energy and moving back is more likely. Thus, with decreasing temperature the probability of oscillating between the two vertices increases until the probability converges to one. This issue can be circumvented by considering double vertices of two neighboring tetrahedra sharing this particular spin perpendicular to the field. By doing so, the energy gain and cost are taken into account within one step.

The weighted random walk of the defect will be performed based on local heat bath sampling. Moving a defect through the smallest unit, either a single tetrahedron or a double tetrahedron, will change its local configuration from $i$ to $j$. The probability is given by

[TABLE]

where $z(i\rightarrow j)=e^{-\beta(\epsilon_{j}-\epsilon_{i})}$ is the Boltzmann weight of a local configurational change. The sum in the denominator runs over all possible configurations $k$ which can be achieved by forward steps and changing the local configuration ($i\neq k$) or by backtracking and retaining the local configuration ($i=k$). The energy $\epsilon_{i}$ of a local configuration is either given by

[TABLE]

where the sum runs over all links of a single vertex $\boxtimes_{\alpha}$, or, for the modified algorithm, by

[TABLE]

summing the Stark energy contribution from the internal spin of a double vertex $D_{\alpha}$. The (external) spins are weighted by a factor $1/2$, since they are shared by two neighboring (double) tetrahedron.

The protocol for the simulation is as follows: An initial configuration is chosen that obeys the ice rules everywhere, loop creations and spin flips at infinite temperature follow in order to randomize the initial state. After then equilibrating at a temperature of choice, the auto-correlation time is measured and determines how many loop iterations will be done between each measurement at the particular temperature. This procedure is repeated for all temperatures of interest. The field is held constant over the whole procedure. Simulations are done on systems with up to $160\times 160\times 80$ hexagonal unit cells and with periodic boundary conditions.

IV Polarization Plateaus

The numerically obtained polarization and the specific heat are shown in Fig. 2 for (i) $\bm{E}_{0}$ along $[001]$ with a tilt $\bm{E}_{\perp}$ along $[100]$, and in Fig. 3 for (ii) $\bm{E}_{0}$ along $[010]$ with a tilt $\bm{E}_{\perp}$ along $[101]$. The applied fields are

[TABLE]

$|\bm{E}_{\perp}|$ is set to be two orders of magnitude smaller than $|\bm{E}_{0}|$.

As a consequence of the two energy scales set by $\bm{E}_{0}$ and the tilt $\bm{E}_{\perp}$, three temperature regimes can be distinguished. For temperatures $T\gg|\bm{E}_{0}|$ the physics is solely determined by the ice rules and, consequently, the polarization is close to zero. Upon cooling down to $|\bm{E}_{\perp}|\ll T<|\bm{E}_{0}|$, a crossover into an intermediate polarization plateau is observed. In the intermediate plateau, the spin sublattice with the maximal projection along the field gets pinned. The other spins remain free and form decoupled, two-dimensional layers. Thus, a dimensional reduction of the physics from three to two dimensions occurs due to the external field. The polarization is determined by the fixed spins and is illustrated as horizontal dashed lines in Fig. 2 and 3.

For both field directions, the crossover between high-temperature and intermediate plateau is in close relation to spin ice.Isakov et al. (2004) Starting from the intermediate plateau, the polarization can only be reduced by introducing string-defects. Such string-defects are composed of concurrent spins and span the entire system in order to not violate the ice rules. Spins along a string-defect are flipped once a temperature is reached at which the gain in entropy compensates for the increase in internal energy induced by the defect. The two-dimensional layers, which are decoupled at low temperatures, begin to interact at this crossover temperature and three-dimensionality gets restored. Further increase of the temperature results in multiple string-defects that exhibit an interaction between string-defects of entropic origin.

A second polarization plateau with saturated polarization is approached at temperatures below the intermediate regime. Near $T_{K}\approx 3.4$ the specific heat diverges for $T>T_{K}$, indicating a continuous phase transition. In contrast, the specific heat vanishes for $T<T_{K}$, which would indicate a first order phase transition. In fact, any fluctuations are suppressed below the transition as they are only caused by string defects spanning the entire system. Such a so called Kasteleyn transition is well known in the context of dimer models.Kasteleyn (1963)

We continue with a more detailed discussion of the intermediate plateau and the Kasteleyn-type transition into the fully polarized phase for the two field directions. A one-to-one mapping from spin configurations within the decoupled layers to hardcore dimers exists. The discussion of the emergent dimer models is given separately for the two field directions studied. The first field direction (i) leads to decoupled layers of dimers on a hexagonal lattice, and the second field direction (ii) gives rise to dimers on a square lattice.

IV.1 Intermediate plateau for $\bm{E}\parallel[001]$: Dimer model on a hexagonal lattice

Applying an electric field along $[001]$ pins the spins parallel to the field, the green spins in Fig. 4(a), as soon as the ratio $|\bm{E}_{0}|/T$ is sufficiently large. The remaining degrees of freedom, the blue spins, reside on the edges of two dimensional honeycomb lattices with alternating two-in-one-out or one-in-two-out vertices111Such a partial ordered state of ice Ih was proposed early on in a different context.Rundle (1953) We can map the blue spins to dimers on the same honeycomb lattice by assigning a dimer to each minority spin, which are highlighted in orange in Fig. 4(b). The hexagonal layers form an ’ABAB’ stacking with inverted two-in-one-out or one-in-two-out vertices. Apart from the stacking, this intermediate state is identical to the well studied kagome ice state of pyrochlore spin ice in a $[111]$ magnetic field.Moessner and Sondhi (2003); Isakov et al. (2004)

Pinning one out of four spin per site reduces the entropy, however, the number of dimer coverings of the hexagonal lattice still grows exponentially and results in a residual entropy. In case of equally weighted dimers, the value of the entropy is reduced to $S/k_{B}\approx 0.1615$ per siteWu (1968) from the total residual entropy of Ih, $S/k_{B}\approx 0.410$ per site.Nagle (1966); Isakov et al. (2004); Vanderstraeten et al. (2018)

Closely packed configurations of dimers on the edges of a honeycomb lattice exhibit a transition between a long-range ordered phase and a disordered phase,Kasteleyn (1963); Wu (1968) The disordered phase exhibits an algebraic decay of two-point correlations.Moessner and Sondhi (2003)

The transition occurs once the Boltzmann weight $e^{-\beta\epsilon_{\alpha}}$ of placing a dimer on one edge exceeds the sum of the weights on the remaining two edges, e.g., $a\geq b+c$. Here, $\beta$ is the inverse temperature, $\beta=1/k_{B}T$, with the Boltzmann constant $k_{B}$ set to unity. $\epsilon_{\alpha}=-\sum_{i\in\alpha}s_{i}\bm{\hat{e}_{i}}\cdot\bm{E}_{\perp}$ is the energy of the local configurations of spins on a vertex $\alpha$. The critical weight can be calculated by considering the corresponding spin configuration for each dimer and the direction of the tilt-component of the field.

For the tilt field under consideration, $\bm{E}_{\perp}\parallel[100]$, the critical condition is

[TABLE]

In our simulation, we normalize the dipole moment of the spin $\bm{d}_{i}=s_{i}\bm{\hat{e}_{i}}$. In doing so, the critical temperature $T_{K}$ reads as

[TABLE]

where $\epsilon_{a}$ is the energy of the spin configuration that is equivalent to the energy of a dimer at the edge with weight $a$.

IV.2 Intermediate plateau for $\bm{E}\parallel[010]$: Dimer model on a square lattice

We now turn to the field direction that has no analog in pyrochlore spin ice. If the field is applied along $[010]$, the intermediate polarization plateau can be understood as layers of hardcore dimers on a rhombic lattice. Each layer is separated by the spin sublattice that is pinned by the field [green spins in Fig. 5(a)]. Thermal fluctuations will occur only within two-dimensional layers.

Each layer is composed of double vertices with two neighboring nodes of the same kind, either both of three-in-one-out or both of one-in-three-out configuration. See also Fig. 5(c). By combining two such nodes to a double vertex one internal spin and four external spins are obtained. Three of the external spins are pointing out (in) and one in (out). Thence, a one-to-one mapping to dimers exists by assigning a dimer to each minority spin. In doing so, one ends up with hardcore dimers on a rhombic lattice, see Fig. 5(c).

Such a dimer model still exhibits an exponentially growing ground state degeneracy upon increasing system size. The residual entropy is reduced, but remains finite and acquires a value of $S/k_{B}=G/\pi\approx 0.291$ per site, where $G$ is Catalan’s constant.Fisher (1961); Kasteleyn (1961)

A tilted field removes the remaining entropy. Similar to the hexagonal case discussed before, a concrete condition for the ordering transition on the square lattice is known:Cohn et al. (2001) the dimers order once the Boltzmann weight $e^{-\beta\epsilon_{i}}$ of a dimer exceeds the sum of the weights of the other three dimer orientations. Note that due to the field, the dimers with the same orientation, horizontal or vertical, come in pairs with inverse weight. By swapping all spins of a given configuration its energy changes sign resulting in an inverse weight.

Assuming the ordering of weights $a>b\geq\frac{1}{b}>\frac{1}{a}$, the critical condition reads $a=b+\frac{1}{b}+\frac{1}{a}$. By inserting the corresponding energies, cf. Eq. (3), we arrive at a transcendental equation

[TABLE]

As before, $k_{B}$ is set to one. For a tilt field $\bm{E}_{\perp}$ applied along $[101]$, we obtain a critical temperature of

[TABLE]

in units of the tilt field $\bm{E}_{\perp}$. $T_{K}$ is shown in the polarization plot, cf. Fig. 3, as a dashed vertical line and agrees with the simulation.

In the following section, the Kasteleyn theory is applied to obtain exact expressions results for the effective dimer model.

V Kasteleyn-Matrix Approach for $\bm{E}\parallel[010]$: Dimers on a Square Lattice

The general goal is to provide an expression for the partition function of hardcore dimers densely packed on edges of a graph Kasteleyn (1963)

[TABLE]

where the Boltzmann weight $z_{i}=e^{-\beta E_{i}}$ belongs to edge $i$ with an occupation number $n_{i}=\{0,1\}$, and $g(~{})$ is a generating function.

The partition sum $\mathcal{Z}$ for dimers on any planar graph can be written in terms of a Pfaffian Kasteleyn (1963)

[TABLE]

of a skew-symmetric matrix $K$. The Pfaffian of a skew-symmetric matrix is defined via the determinant as $(\mathrm{Pf}~{}A)^{2}=\det A$. The matrix $K$ is the directed adjacency matrix of a graph—here lattices—with weights $z_{i}$ and an additional phase factor $\phi_{i}$. $\phi_{i}$ has to be chosen such that a loop around a face with $0\mod 4$ edges acquires a negative sign and a positive sign around a face with $2\mod 4$ edges,Kasteleyn (1963) This ensures, that all the terms generated in $\mathrm{Pf}~{}K$ have the same sign and Eq. (9) is valid.

A simplification exists for bipartite lattices: If the nodes of the lattice are numbered such that the nodes of sublattice A get indices $\{0,1,\ldots,N/2-1\}$ and the nodes of B indices $\{N/2,\ldots,N-1\}$, then $K$ has the structure

[TABLE]

Here, $\tilde{K}$ denotes the Kasteleyn matrix, which is the (reduced) adjacency matrix. The matrix $\tilde{K}$ incorporates only the connections from even to odd sites, i.e. rows of the adjacency matrix represent nodes of sublattice A and columns represent nodes of B. In contrast to $K$, the reduced $\tilde{K}$ is not defined with respect to a directed lattice.

Since $\det~{}K=(\det~{}-\tilde{K}^{T})(\det~{}\tilde{K})=(\det~{}\tilde{K})^{2}$, Eq. (9) simplifies to

[TABLE]

Figure 6 illustrates the convention of the fundamental unit cell and the edge weights. With this choice, $\tilde{K}$ reduces to a scalar and reads

[TABLE]

with complex variables $z=e^{i\phi}$ and $w=e^{i\psi}$.

Following Refs. [Kenyon et al., 2006; Kenyon, 2009], the partition sum, $\mathcal{Z}$, can be expressed as a double integral in the complex plane

[TABLE]

A direct computation of (13) for arbitrary weights $a$ and $b$ is difficult, if feasible at all. The logarithm, $\log z$, of a complex variable $z$ has a branch cut for $\{z\in\mathbb{R},z\leq 0\}$ and, thus, the residue theorem cannot be employed.

However, given the partition sum, thermodynamic quantities are defined as derivatives of $\log\mathcal{Z}$, resulting in terms suitable for the residue theorem. Let us focus on the density of dimers on a horizontal bond. If a dimer is placed right or left of an even site (filled dot in Fig. 6, its weight is inverted, $a\rightarrow a^{-1}$. Thence, the densities on either of the horizontal bond are correlated and its difference can be derived

[TABLE]

By moving the derivative $\partial/\partial a$ inside—the integration paths do not depend on $a$—we get

[TABLE]

The inner integral with respect to $z$ is performed using the residue theorem. First order poles are situated at

[TABLE]

of which $z_{0}$ is always within the integration path, and $z_{1}(w)$ depends on the second integration variable. Along the integration path, $|w|=1$, of the outer integral, $z_{1}(w)|_{|w|=1}$ is either completely inside, partially inside, or completely outside of the integration path, $|z|=1$, of the inner integration.

See Fig. 7 for an illustration of the three different cases. Which of the three cases applies depends on the choice of the weights $a,b$. The cases of $z_{1}(w)$ being either completely inside or outside are trivial in that the residue theorem can also be applied for the integration with respect to $w$. These cases occur if one of the weights is larger than the three remaining ones, e.g., $a>b+\frac{1}{b}+\frac{1}{a}$, and the dimer configuration is long-range ordered. Consequently, the density $\rho_{a}$ is constant

[TABLE]

In the remaining case the pole $z_{1}(w)$ is partially inside $|z|=1$. Thence, the outer integral does only run along the path of $w$ for which $|z_{1}(w)|<1$ resulting in

[TABLE]

The first term comes from the residue at $z_{0}=0$, and the second term from the residue at $z_{1}(w)$. The first term is easily solved using the residue theorem again: a pole $w_{0}=b/a$ is either inside or outside $|w|=1$ and results in a step function $\Theta(x)$, which is zero if $x\leq 0$ and one if $x>0$. For the second term, an expression of the indefinite integral exists. Consequently, the result of the definite integral is simply the difference of the indefinite integral at the final, $w_{f}=e^{i\phi_{0}}$, and the initial point, $w_{i}=e^{-i\phi_{0}}$, of the path. Thence,

[TABLE]

with

[TABLE]

being the angle at which $|z_{1}(e^{\pm i\phi})|=1$.

By mapping back to the spin model, the polarization is obtained as $\bm{P}=d\sum_{i}\left<s_{i}\right>\bm{\hat{e}_{i}}$, where $\left<s_{i}\right>$ is the thermal average of the Ising variable $s_{i}$, which depends on $\rho_{a}(a,b)$ and $\rho_{b}(a,b)$. Figs. 2 and 3 include both the numerical results and the exact values for the polarisation in the thermodynamic limit. The latter exhibits a sharp cusp at the transition temperature. This cusp rounds off within the numerical simulations with a more pronounced rounding effect for small systems.

We now turn to a discussion of the specific heat, $C$, which is obtained in two ways: firstly from the Kasteleyn-matrix approach by $C=\partial E/\partial T$, where $E$ is the energy due to the Stark-coupling to the external field, and secondly from the fluctuations of the energy within the directed-loop Monte Carlo simulation. For $T>T_{K}$, $C$ follows a critical behaviour as $C\propto(T-T_{K})^{-1/2}$ for a range of $T-T_{K}<0.1$, cf. Fig. 8. The critical exponent, $\alpha=1/2$, is in agreement with the literature.Kasteleyn (1963); Bhattacharjee et al. (1983)

For the finite system studied numerically, the transition shifts slightly. Beyond the transition of the finite system, the numerically obtained specific heat has the same critical bahaviour as extracted from the Kasteleyn-matrix approach. For $T<T_{K}$ the specific heat drops to zero if being in the thermodynamic limit. In finite systems, however, fluctuations exist near $T_{K}$ following an Arrhenius-type activation energy. The energy is given by the total cost of shifting subsequent dimers along a string spanning the entire system and, thus, are suppressed exponentially in system size and inverse temperature.

In conclusion, the specific heat exhibits a strongly asymmetric behaviour typical for a Kasteleyn transition: the specific heat diverges as $(T-T_{K})^{-1/2}$ for $T>T_{K}$, whereas $C=0$ at $T<T_{K}$ since fluctuations are fully suppressed due to the hard-dimer constraint that in our case originates from the ice rules. Finite size systems mitigate this asymmetry slightly as the energy cost to introduce defects decreases for smaller systems.

VI Correlations and Static Structure Factor

In this section, we study the static structure factor (SSF), which allows for comparison with scattering experiments like inelastic neutron scattering. The SSF is the Fourier transform of the equal-time two-point correlations and reveals details about the microscopic degrees of freedom and how they are correlated. Here, we focus on the diffuse part of the SSF capturing the contribution from the proton disorder. The static oxygen lattice would only contribute to the Bragg peaks and is neglected.

Following the notation of Ref. [Li et al., 1994], the general structure factor of the protons is

[TABLE]

where, the sum runs over all proton sites, $\bm{R}_{j}$, of which two per link exist. $n_{j}$ denotes the corresponding occupation number, which can either be zero or one.

By making use of the first ice rule, i.e. one proton per link, the proton configuration can be rewritten in terms of Ising variables $s_{j}$ residing at the middle of a link, $\bm{r}_{j}$, and the proton shift distance combined with the direction of the link, $\bm{a}_{\gamma(j)}$. Here, $\gamma(j)$ maps the general label $j$ of a link to one of the distinct links with label $\beta$ within the unit-cell. Hence, $\bm{a}_{\gamma(j)}$ is a purely geometrical factor, whereas the Ising variables $s_{j}$ contains all the information about the degrees of freedom and their correlations. Then, Eq. (20) becomes

[TABLE]

where the sum runs over all links $j$. The sine factor leads to a suppression of any diffuse scattering within the first and partly the second Brillouin zone.

Equation (21) simplifies when considering only effective Ising spins $s_{j}$ and neglecting the geometrical factor $\bm{a}_{\gamma(j)}$

[TABLE]

We fix the convention for the $s_{j}$ such that a positive $s_{j}$ refers to a spin pointing from a site on the sublattice A to a site on B. Even though the two SSF look very different, they include the same information on the correlations and contain similar scattering features.

In the following, we discuss the SSF in the different regimes. The discussion is separated into the high-temperature regime and the intermediate regimes depending on the field direction.

VI.1 SSF at high temperatures

The SSF at high temperatures or with small external field, respectively, is shown in Fig. 9, where thermal fluctuations exceed the Stark energy and the SSF is comparable to the zero-field case. Good agreement is found with experimental data based on neutron scattering Li et al. (1994) and the theoretical analysis done by Wehinger et al. Wehinger et al. (2014), Isakov et al. Isakov et al. (2015), and Benton et al. Benton et al. (2016).

Both SSF—that of the protons, $S_{H}(\bm{k})$, as well as that of the spins, $S_{I}(\bm{k})$—exhibit pinch points with a bow tie structure of the intensity around lattice Bragg points, e.g. at $\bm{k}=(2,~{}2,~{}0)$ or $(0,~{}4,~{}2)$ for $S_{H}(\bm{k})$ and at $\bm{k}=(0,~{}2,~{}1)$ or $(0,~{}0,~{}2)$ for $S_{I}(\bm{k})$. The bow tie structure is typical for local constraints like the ice rules, that are equivalent to a divergence free condition on a lattice. The ice rules translate to Gauss law, $\nabla\cdot\bm{E}=0$, that eliminates the longitudinal component of the field. Regarding the microscopic degrees, the ice rules project out any longitudinal correlation of the otherwise unconstrained proton displacements. Broad features occur, that are most noticeable in $S_{I}(\bm{k})$ within the $[hk0]$-plane, e.g. near $\bm{k}=(1,~{}0,~{}1)$. These broad features are the diffuse remnants of pinch points in the $[hk1]$-plane and become pinch points in the effective dimer model of the intermediate regime, due to the dimensional reduction to two dimensions.

VI.2 SSF Within Intermediate Plateau for $\bm{E}\parallel[001]$

In the intermediate plateau upon applying a field along the $[001]$ axis, the SSF, exhibits the sixfold symmetry of the honeycomb lattice the remaining degrees of freedom are confined to, cf. Fig. 10. Pinch points with bow ties remain in the presence of the field. The pinch points occur due to the local constraints of dimers, which itself is a consequence of the initially imposed ice rules.

More strikingly, the SSF exhibits peaks at $\bm{k}=(4/3,~{}0,~{}0)$ (in $S_{I}(\bm{k})$) and equivalent points generated by the $C_{6}$ symmetry. These peaks are caused by algebraically decaying correlations, $C(r)\propto\frac{1}{r^{\alpha}}$, with $\alpha=2$ for classical dimers on the honeycomb lattice.Moessner and Sondhi (2003) As such the intensity of these peaks grows logarithmically with system size. The correlations given in [Moessner and Sondhi, 2003] lead to a similar SSF as in case of the effective spin model,222Yet, the structure factor plotted in [Moessner and Sondhi, 2003] differs. The intensity at and near the peaks ($k=(4/3,~{}0,~{}0)$ and symmetry-equivalent points) is suppressed in their plot. We obtain the same SSF, if we change a sign in the phase factor $e^{i\bm{k}\bm{r}}$ of the Fourier transform. cf. Fig. 10 (a) and (c).

When considering the SSF of the protons, cf. Fig. 10(b), the form factor due to the proton positions shift the pinch points and peaks, yet, the qualitative features remain similar.

Upon lowering the temperature, the effect of the tilt field becomes apparent, cf. Fig. 10(d). The peaks split and shift towards a nearby pinch point. Two peaks superimpose in the SSF stemming from ABAB stacking of honeycomb layers. The tilt field prefers different edges in the two layers. Consequently, one peak moves towards the $\bm{k}=(1,~{}1,~{}0)$ pinch point and a second peak moves towards the $\bm{k}=(2,~{}0,~{}0)$ pinch point. This is illustrated in Fig. 10(d), where the SSF in the left half is obtained by summing only over the spin variables in the lower honeycomb layer, whereas the right half includes only the upper honeycomb layer. Upon transitioning into the ordered phase the peak merges at the pinch point with another peak. The two peaks form a point-symmetric pair with respect to the same pinch point.

VI.3 SSF Within Intermediate Plateau for $\bm{E}\parallel[010]$

The SSF within the intermediate polarization plateau upon applying a field along the $[010]$ axis is shown in Fig. 11. The SSF of the effective spin model, $S_{I}(\bm{k})$ as in Eq. (22), is shown in (a), and the SSF of the full proton system, $S_{H}(\bm{k})$ as in Eq. (21)), is shown in (b).

Regardless of the choice, pinch points with a bow-tie like structure occur as a signature of the ice rules. Some of the initially broad features at high temperature sharpen and become pinch points, e.g., at $\bm{k}=(3,~{}0,~{}1)$ for $S_{H}(\bm{k})$ and at $\bm{k}=(1,~{}0,~{}1)$ for $S_{I}(\bm{k})$. Similar to the first field direction, the SSF exhibits peaks in between the pinch points, which are caused by the algebraic decay of correlations.

In order to make the similarity to a dimer model more striking, let us consider a SSF, where only the external spins of a double-vertex are considered. In doing so, the reduced SSF $S_{I,red}(\bm{k})$ becomes equivalent to the SSF of dimers on a square lattice.Fisher and Stephenson (1963); Youngblood et al. (1980); Fradkin et al. (2004); Tang et al. (2011) The asymptotic correlations of dimers on a square lattice have been calculated for a smaller unit cell with only two Boltzmann weights, one for horizontal and one for vertical dimers.Youngblood et al. (1980) Nevertheless, we can use the result in the limit of equal weights which corresponds to not applying a tilt field. Fig. 11 allows for a comparison of the structure factor based on our simulation (c) with the one based on the asymptotic correlations (d). Most notably, both show the peaks at $\bm{k}=\left(\frac{1}{2}+n,~{}0,~{}\frac{1}{2}+m\right)$ with $n,m\in\mathbb{Z}$. They are caused by algebraic correlation of dimers of the same orientation Tang et al. (2011) decaying algebraically with $1/r^{2}$. The broad square-like feature of increased intensity has its origin in the correlation of dimers of different orientation. Thus, we clearly find a dimensional reduction of the three-dimensional ice to a tow-dimensional model of dimers on the edges of a square-like lattice.

As a remark, due to combining two oxygen nodes of different neighboring sublattices of the original ice Ih lattice, the initially defined orientation of two of the external spins is opposite to the other two spins. This leads to a shift of $2\pi$ along the direction of alternation. In our case this shift occurs along $[001]$ and has already been considered in Fig. 11(c).

A second difference becomes apparent, when considering the enlarged unit cell occurring here. In an earlier work, only two different Boltzmann weights were considered: one weight for horizontal and one weight for vertical dimers.Youngblood et al. (1980) In ice Ih the bipartite property of the lattice is reflected in the dimer weights and we have to double the unit cell consisting instead of two nodes and four inequivalent edges. A rotation of $\pi/4$ is needed in order to convert two of the small unit cells into the bigger one. The reciprocal lattice vectors change accordingly in length and direction. The position of the features in reciprocal space obtained in Ref. [Tang et al., 2011] with the positions in the reciprocal space of ice Ih, cf. Table 1. The diffuse scattering is not affected apart from the shift and rotation.

Upon applying the tilt field, the logarithmically divergent peaks split into a pair of peaks due to the doubled unit cell and the two species of horizontal and vertical dimers. Furthermore, the weight of the peaks becomes asymmetric since one of the two dimers of each orientation is energetically preferred. Its occupation number and the intensity of the corresponding peak increases. Upon further decreasing the temperature and approaching the Kasteleyn transition, two peaks merge at a pinch point. This behaviour is not captured by the asymptotic correlations derived in Ref. [Youngblood et al., 1980].

By further reducing the temperature below the Kasteleyn temperature the system exhibits a transition into a single ordered configuration, which does not exhibit any diffuse scattering.

VII Conclusion

We studied a minimal model of the proton subsystem of hexagonal water ice Ih coupled to an electric field. The field is applied along two different easy axes, resulting in qualitative different microscopic descriptions in terms of emergent dimers. Upon slightly tilting the field, the proton configuration enters a long-range ordered phase with zero entropy. The transition into the ordered state is of Kasteleyn-type.

In case of applying the field along the $[001]$ axis, a dimer model on the honeycomb lattice emerges. Whereas, upon applying the field along the $[010]$ axis, the remaining degrees of freedom are described by a dimer model on the square lattice. The first case is in close analogy to the intermediate magnetization plateau in spin ice with an magnetic field along $[111]$. Its transition to the ordered state is determined in the original work of Kasteleyn.Kasteleyn (1963) The second case is a result of the structure of ice Ih and does not have an analogy to pyrochlore spin ice. The transition into an ordered state upon tilting the field is also described by the Kasteleyn matrix. An analytic expression for the dimer densities is given, which is directly related to polarisation, the internal energy, and the specific heat. The polarization has a cusp at the Kasteleyn transition and the specific heat exhibits a one-sided divergence with a critical exponent of $\alpha=1/2$ for $T>T_{K}$.

In both cases, the application of a field leads to a dimensional reduction of the thermally fluctuating degrees of freedom with a simultaneously reduced power in the algebraic correlations. In the effective two-dimensional models of dimers the correlations decay with $1/r^{2}$ and result in peaks in the static structure factor. Due to the algebraic decay, these peaks diverge logarithmically with system size.

For a sharp Kasteleyn transition to occur, strings of concurrent spins have to span the entire system along which the protons can be swapped without breaking the ice rules. Those dynamics are unlikely to occur in a realistic ice sample, as the energy barrier is to high and thermal activation is suppressed exponentially. However, by introducing defects, i.e. doping with KOH, the necessary dynamics can be restored, but the defects will also round the characteristic cusp of the polarization at the Kasteleyn transition.Jaubert et al. (2008)

Correlated quantum tunneling processes were discussed Drechsel-Grau and Marx (2014); Benton et al. (2016) and may have been observed experimentally Bove et al. (2009); Yen and Gao (2015) in ice Ih at zero field. Within the intermediate polarization plateau, the partial ordering restricts the quantum tunneling of protons into two-dimensional layers. The emergent dimer models change to quantum dimer models with local ring exchange terms. For a broad class of quantum dimer models on bipartite lattices, an ordered phase is to be expected,Moessner and Sondhi (2001) due to the quantum order-by-disorder mechanism. Depending on whether rotation and/or translation symmetry of the ground state is broken, e.g. staggered, columnar, or plaquette dimer configuration, the corresponding proton configurations of each layer can be qualitatively different, e.g. ferroelectric, antiferroelectric, or paraelectric. In contrast, the saturated plateau is completely frozen and quantum tunneling will be suppressed completely.

VIII Acknowledgements

We are grateful to Ludovic Jaubert, and Nic Shannon for insightful discussion. This work was supported in part by Deutsche Forschungsgemeinschaft (DFG) via SFB 1143, the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter – ct.qmat (EXC 2147, project-id 39085490), and Germany’s Excellence Strategy via EXC 2111 (project-id 390814868), FP acknowledges the support of the DFG Research Unit FOR 1807 through grants no. PO 1370/2- 1, TRR80, the Nanosystems Initiative Munich (NIM) by the German Excellence Initiative, and the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement no. 771537).