Crystalline phases and devil’s staircase in qubit spin ice

,


I. INTRODUCTION
Water ice is a beautiful example of how an everyday material can inspire and advance many areas of physics, from classical statistical physics to qubit systems.Oxygen ions in ice form a four-fold coordinated lattice.The protons are located on the links connecting the nearest neighbor oxygen ions: two protons bond to the oxygen covalently, and two with hydrogen bonds.It is the ice rule identified by Bernal and Fowler in 1933 [1].Because the lengths of the covalent and hydrogen bonds are different, arrows can be used to show the position of the protons about the center of the link.It leads to the "twoin, two-out" formulation of the ice rule, which can be considered a local divergence-free condition.The ice rule allows for six different proton configurations around an oxygen ion.Because these six configurations correspond to six different vertices in the arrow representation [see Fig. 1(a)], the fundamental model of ice is known as the six-vertex model (6VM).Remarkably, the number of states satisfying the ice rule increases exponentially with the size of the system, forming a In the fully packed loop representation, the thick bonds are occupied (i.e., part of a loop), whereas the thin bonds are empty.The figure shows the correspondence between the arrows and the occupied or vacant bonds for "even" sites, while the occupancies are reversed for "odd" sites.The "even/odd" refers to the parity of the sum ix + iy of a site's coordinates (ix, iy) in the lattice.(c) In the alternative Baxter representation, the bonds are colored blue if the arrows on the bond are reversed compared to the reference vertex in the first column.
manifold.Pauling estimated the degeneracy of the manifold as W ice = (3/2) N = 1.5 N , where N is the number of oxygen ions, resulting in a finite residual entropy [2].Lieb solved the two-dimensional six-vertex model on the square lattice exactly and got W 2D = (4/3) 3N/2 ≈ 1.5396 N [3], which is very close to Pauling's estimate.Furthermore, Baxter noted that the correlations decay algebraically [4].
The six-vertex model successfully describes physical systems in which divergence-free conditions arise.For instance, Anderson proposed that the frustration in magnetite leads to charge disproportionation, where configurations following the Bernal-Fowler rule minimize the Coulomb energy [5].Another example is the "square ice" substance KH 2 PO 4 (KDP), a quasi 2D material in which the vertex configurations do not have equal energies [6,7].In particular, the discovery of spin ice materials brought the field to flourish [8].In spin ice, the Ising-like magnetic moments of rare earth ions form the highly frustrated pyrochlore lattice of corner-sharing tetrahedra with spins at the corners.Arrows representing these spins realize the low energy two-in, two-out configurations [9].Since then, spin ice physics was also accomplished in fabricated arrays of nanomagnets [10][11][12][13].
The six-vertex model is a convenient starting point for realizing topological defects due to its correlated ground state manifold.Flipping an arrow creates a "three-in, one-out" and a "one-in, three-out" vertex.These vertices are not part of the six-vertex model; the underlying model (e.g., Ising model) determines their dynamics.They correspond to fractional charges in the model for magnetite [14] and to the emergence of magnetic monopoles in the spin ice systems [15], experimentally confirmed in Ref. [16].
The quantum six-vertex model emerges by allowing tunneling between configurations that obey the ice rule by adding an off-diagonal term to the classical Hamiltonian.This term reverses arrows arranged tail-to-nose around an elementary square plaquette, as shown in Fig. 2a, though longer loops are also allowed.Chakravarty used such an expression to describe d-wave superconductors in Ref. 17.
The quantum six-vertex model appears in the perturbative expansion of the S = 1/2 Heisenberg model on lattices of corner-shared tetrahedra in the limit of large easyaxis exchange anisotropy.For example, on the 2-dimensional checkerboard lattice, the quantum term gives rise to a gapped phase, where arrows on alternating square plaquettes resonate [18,19].On the 3-dimensional pyrochlore lattice, Hermele et al. argued that the effective theory is a Maxwellian U (1) action with gapless "photon"-like excitations [20], confirmed numerically in Refs.21-23.This situation may arise in certain spin-ice materials: Tb 2 Ti 2 O 7 [24], Pr 2 Sn 2 O 7 and Pr 2 Zr 2 O 7 [25], and Yb 2 Ti 2 O 7 [26] are all suitable candidates (see Ref. 27 for a review about the quantum spin ice).The stability of the U (1) spin liquid against ordered phases in quantum spin ice was considered in Refs.28 and 29, together with experimental signatures.Other examples include the isotropic S = 1/2 Heisenberg model with four-site ring exchange on checkerboard and pyrochlore lattices, with the constraint of exactly one singlet bond on each tetrahedron [30].By extending the fundamental model, Ref. [31] investigated the quantum effects in the "square ice" KH 2 PO 4 .Ref. [32] studied finite temperature effects in quantum square ice.More recently, the dynamical properties of the model came under scrutiny: it exhibits dynamical quantum phase transitions [33] and quantum manybody scars [34].
The connection between gauge theories and the quantum six-vertex model, initially discussed in Ref. [20], was further explored in Ref. [35], where the two-dimensional quantum six-vertex model was found to be a confining lattice gauge model.The model is also known as the (2+1)-dimensional U (1) quantum link model [36].Refs.37 and 38 considered the model from a gauge field theory point of view.This relationship inspired the concept of engineering arrays of Rydberg atoms as simulators of U (1) lattice gauge theories in various geometries [39][40][41].
The motivation for our research comes from the recent implementation of the quantum six-vertex model in a quantum annealing system by King et al. [42].Their setup consisted of superconducting flux qubits arranged in an array that physically realized the transverse-field Ising model on the checkerboard lattice.Four ferromagnetically coupled qubits formed a single logical spin, representing an Ising spin.Antiferromagnetic two-body couplers between qubits belonging to adjacent logical spins provided a tunable antiferromagnetic interaction between the Ising spins.They implemented two inequivalent couplers that enabled tuning the parameters of the Ising model into the range described by the six-vertex model and lifting the degeneracy between type-I and type-II vertices (see Fig. 1a).Quantum fluctuations induced by the transverse field led to tunneling between the six-vertex configurations.Thus, the V/t isolated states resonating plaquettes fully flippable 1 −0.37 FIG. 3. The phase diagram of the quantum six-vertex model, defined by Hamiltonian (3), consists of three phases.For V ≲ −0.37t, the twofold degenerate fully flippable phase arises, where the number of flippable plaquettes is maximal.Every second plaquette (denoted by magenta circles) resonates when −0.37t ≲ V ⩽ t.The approximate wave function is a direct product of the |⃝⟩ = 1 .The Rokshar-Kivelson point at V = t is quantum critical, and for t < V , the isolated manifold appears with configurations having no flippable plaquettes (denoted by crosses).minimal model of their setup involves the tunneling term where the sum is over the elementary square plaquettes of the lattice, and |⟳⟩ and |⟲⟩ denote square plaquettes with the clockwise and anticlockwise orientation shown in Fig. 2(a), and a chemical potential to distinguish the two types of vertices (the NII operator counts the number of the type-II vertices).It is also convenient to introduce the term, where NV counts the number of flippable plaquettes.
The full model we will consider in this paper is then Let us briefly review the known limiting cases of the Hamiltonian.
H II is the Hamiltonian of the Rys F model [43].It is a wellknown problem in statistical physics [4] and exactly solvable by Bethe Ansatz [44,45].It exhibits two phases at zero temperature.For µ > 0, the two-fold degenerate ground state consists of alternating type-I vertices -this is the antiferroelectric phase.If µ < 0, configurations with only type-II vertices span the disordered phase's highly degenerate ground state manifold.
Following the footsteps of Rokshar and Kivelson [46], Shannon et al. went beyond the pure quantum six vertex model H t of Chakravarty [17] and introduced the Hamiltonian [18].In the fully-packed loop representation, the Hamiltonian (3) has precisely the same form as the quantumdimer model of Rokshar and Kivelson, except for the Hilbert space: it acts on dimers in one case and loops in the other.Just like in the quantum-dimer model, the exact ground state of the model is an equal-weight superposition of all connected states when V = t.This is a quantum critical point with algebraically decaying correlations.It separates the subextensively degenerate ground state manifold of isolated (also called disconnected) states from the resonating plaquette phase (see Fig. 3 for a sketch of the phase diagram).Configurations in the isolated manifold consist of type-II vertices only and have no flippable plaquettes; thus, H t annihilates them.The resonating plaquette phase is similar to the one in the quantum dimer model [47], but every alternating square plaquette resonates, so it is two-fold degenerate only.Exact diagonalization studies estimated the lower boundary of the plaquette phase as V /t ≈ −0.3727 [18] and V /t = −0.359(5)[37], a gauge-invariant matrix product states calculation located the transition point at V /t = −0.37(3)[38], and quantum Monte Carlo at V /t = −0.35(3)[41].Below this boundary, flippable plaquettes minimize the energy but without resonance.One can think of this "fully flippable phase" as the antiferroelectric phase of the Rys F model dressed with quantum fluctuations.This phase corresponds to the Néel phase in the XXZ model on the checkerboard lattice [18].
The present study aims to extend the phase diagram introduced above and shown in Fig. 3 by including the term with the chemical potential for the type-II vertices, Eq. (1b), noting that for V = 0 and µ > 0, the plaquette phase is known to persist up to µ/t = 0.288 [19].We will derive the phase diagram of the quantum spin ice Hamiltonian (2) in the complete V -µ plane and the structure factor at zero temperature.It will allow us to get an insight into the results of the qubit-engineered quantum spin ice of King et al. [42].
The paper is organized as follows.We describe the classical six-vertex configurations in various representations and the flux sectors in finite-size clusters with periodic boundary conditions in Sec.II.In Sec.III, we construct the t = 0 classical phase diagram of the model.In Sec.IV, we systematically classify the symmetries of the model, construct order parameters, and write the Landau free energy for phases with zero topological flux.Sec.V discusses the properties of isolated states.In Sec.VI, we present numerical results (exact diagonalization) to reveal the ground state phase diagram of the quantum model.As an independent check, in Sec.VII, we use perturbation theory to calculate corrections to the ground state energies of the classical phases and deduce some of the phase boundaries.In Sec.VIII, we sample the wave function at the Rokhsar-Kivelson point with a Monte Carlo method and explore the phases emanating from this quantum critical point.We also characterize the emergent quantum electrodynamics.Structure factors in different phases are evaluated in Sec.IX and compared with the ones observed in qubit quantum spin ice.We conclude with a summary of results in Sec.X.Finally, appendices A-F contain some details of our calculations.

A. Finite clusters
We study the six-vertex model on finite clusters with periodic boundary conditions on the square lattice.Their size N and geometrical symmetries characterize these clusters.We focus on two families having the full D 4 point group symmetry of the lattice.We refer to the ones generated by the g 1 = (L, 0) and g 2 = (0, L) lattice vectors as the N = L 2 family, see Fig. 4(a).The g 1 = (L, L) and g 2 = (−L, L) lattice vectors define the N = 2L 2 family, shown in Fig. 4(c).The periodicity of the ground states, as we will see later on, requires even values for L. We consider lattice sites translated by an integer multiple of the g 1 and g 2 identical, r + n 1 g 1 + n 2 g 2 ∼ = r, where n 1 , n 2 ∈ Z.

B. Representations of the 6-vertex configurations
The classical 6-vertex model has a long history and applies to many systems, each prompting a convenient representation.Below we review some of them.
a. Arrow representation: This is the original representation of water ice.A configuration is represented as a directed graph, with arrows showing the direction of the edges (bonds).To satisfy the ice rule, every vertex has two inward and two outward pointing arrows, demonstrated in Fig. 1

(a).
A plaquette is flippable if the arrows around the elementary square point clockwise or counterclockwise.The flip itself corresponds to changing the directions of the arrows around a plaquette; see Fig. 2(a).This representation is meaningful for calculating neutron scattering cross section detailed in Sec.IX. b.Fully packed loop representation: In Ref. [5], Anderson described magnetite as a charge-frustrated material using the Ising model.The charge frustrated Fe +2.5 build a pyrochlore lattice consisting of corner-shared tetrahedra.The minimal Coulomb energy corresponds to two 2.5 + δq and two 2.5 − δq charged ions on each tetrahedron.These are represented as occupied and empty bonds, shown in Fig. 1(b) for the two-dimensional model on the checkerboard lattice.Identically charged bonds form closed loops in a finite system.A plaquette is flippable if, as we go around it, oppositely charged bonds meet at each vertex; see Fig. 2(b).The bonds exchange their charges by a plaquette flip, just like in the quantum dimer model [48] describing short-range resonating valence bonds.
We measure the occupancy of a bond by n r = ±1, where r is the coordinate of the center of the bond.Assuming that the horizontal bonds are along the x and the vertical along the y direction, the following relations hold between the arrow and fully packed loop representations: The integer-valued (i x , i y ) are the coordinates of the vertices, and the bond lengths are set to 1. c. Baxter-and alternative Baxter-representation: In his textbook [4], Baxter chose an isolated configuration formed by identically oriented horizontal and vertical bonds as a reference configuration, shown in Fig. 11(a).Then he highlighted all the bonds in a configuration that pointed in the opposite direction compared to the reference.Here we use the same principle but choose one of the fully-flippable configurations as the reference (therefore, we call it the alternative Baxter representation); see Fig. 1(c).For instance, a vertex is type-II if two highlighted bonds meet there, and a plaquette is nonflippable if it has both highlighted and non-highlighted bonds [Fig.2(c)].This representation helps us to identify the mathematical structure of the configuration space.
d. Faraday loop representation: Type-II vertices can be associated with local dipole moments.Drawing these dipole moments as arrows, they form closed loops in ice rule obeying configurations [49].They help study the thermodynamic properties of ice systems and provide a way to approach magnetic monopoles.
e. Height representation: The local divergence-free constraint at the vertices enables us to transcribe an arrow configuration to integers on the plaquettes.Since we are not using it in our work, we only refer to Refs.[50] for details.
All of the representations above constitute a basis where both the NV and NII operators are diagonal, and the quantum flipping term H t is strictly off-diagonal.We refer the reader to Ref. [51] for a comprehensive account of the various sixvertex model representations.

C. Flux sectors
In a cluster with periodic boundary conditions, for each 6-vertex configuration, we can count the net flux of arrows through any given vertical (m x ) or horizontal (m y ) cut.Since the local flips do not change the net flux, the vector m = (m x , m y ) defines a set of winding numbers that the Hamiltonian conserves.States having the same index pair (m x , m y ) form a flux sector.We can generate configurations in different flux sectors by flipping arrows on a directed loop crossing the cluster's boundaries, as illustrated in Fig. 4(b).In the N = 2L 2 clusters, the horizontal and vertical cuts are the diagonals of the rotated square, and the minimal nonzero flux sector is the m = (2, 2) shown in Fig. 4(c).Let us note that the m defining the flux sectors is proportional to the total magnetization in the Faraday loop description [49], and the quantum term mixes all the configurations having the same total magnetization (except the isolated states).
Besides the geometric symmetries described by the point and translation group of the lattice or cluster, there is an internal symmetry, the charge conjugation C [37].It reverses the occupation of the bonds in the fully packed loop representation and commutes with the Hamiltonian, [C, H] = 0.In the arrow representation, it reverses the direction of all the arrows.As a consequence, the flux sector of a configuration changes signs under charge conjugation, Cm = −m.In the alternative Baxter representation, C changes highlighted edges into unlighted ones and vice versa.
We mostly use numerical means to calculate the ground state properties.To this end, we shall generate all the possible ice-rule configurations in a given cluster.Based on empirical findings on small clusters, we assume that the flux sectors are ergodic if their classical states have at least one flippable plaquette (ergodicity has been proven for the m = (0, 0) flux sector in Ref. 20).Therefore, it is enough to find a configuration from a flux sector, since applying local flips will generate all the configurations within the sector.In the case of the  8, the four square (Sq) configurations of Fig. 9, and the isolated manifold, shown Fig. 11(a) and (c), are located at the corners.N = L 2 clusters, one can find a systematic way to construct initial configurations using only type-II vertices.In these configurations, the directions of the arrows along a horizontal or vertical line are all the same (but the directions may differ from line to line).Turning all the arrows on a line changes the flux sector by one unit, allowing access to the desired flux sector.
In Fig. 5, we present the number of configurations in each flux sector for the N = 16 site cluster.The dimension of the Hilbert space in the m = (0, 0) flux sector is a modest 990.Data for larger site clusters are presented in Fig. 25 in Appendix A. We just note that the dimensions of the (0, 0) flux sectors are 962 734 for N = 32 and 5 482 716 for N = 36, these are easy to diagonalize by the Lánczos method.
Plotting the possible N V and N II values of the configurations, we find these values are not independent.Fig. 6 shows a map for the N = 32 and N = 36 site clusters.The ice rule and the periodic boundary conditions constrain the number of allowed type-II vertices and flippable plaquettes to a triangle The classical phase diagram in the parameter space of the V and µ, the chemical potential of the type-II vertices.The isolated and square phases consist only of type-II vertices; they become favorable when µ is negative and V is positive.The fully flippable phase consists of type-I vertices, and all plaquettes are flippable, gaining energy when V is negative.There are two fully flippable states (shown in Fig. 8) and four square states (see Fig. 9).The degeneracy of the isolated phase is subextensive and increases exponentially with the linear size of the cluster.The boundary between the square phase and isolated manifold (thick red line) hosts the disordered manifold of the Rys F model.
in the N II -N V plane.In Appendix B we derive the inequalities which determine the triangle boundaries for a cluster with N vertices (i.e., N sites).

III. CLASSICAL PHASE DIAGRAM
Below we derive the phase diagram in the classical limit of the Q6VM where t vanishes.The Hamiltonian is diagonal in the basis of both the six-vertex and fully packed loop configurations shown in Fig. 1.The energy of a configuration depends only on the number of type-II vertices N II and flippable plaquettes N V as Since the energy in Eq. ( 6) is linear in both N II and N V , three phases emerge in the minimization procedure corresponding to the three corners of the triangle.Fig. 7 shows the classical phase diagram.The first phase consists of isolated configurations having no flippable plaquettes and only type-II vertices so that (N V , N II ) = (0, N ).The energy is then We devote section V to the properties of the isolated manifold.The second one is the fully flippable phase (FF).It maximizes the number of the flippable plaquettes and contains type-I vertices only (Fig. 8), so that (N V , N II ) = (N, 0) and the energy is Equating the two energies above, we get the V = µ phase boundary between the fully flippable phase and the isolated manifold.It results in an extensively degenerate boundary carrying configurations from the side of the triangle by N V +N II = N .Its degeneracy may allow the quantum term to induce further phases.The fully flippable and isolated phases appeared in the isotropic (i.e., µ = 0) limit of the Q6VM studied in Ref. [18] as the doubly-degenerate Néel and the sub-extensively-degenerate quasi-collinear phase.
In addition to these known phases, we identified a third classical phase called the square phase (Fig. 9), where only half of the plaquettes are flippable, and all the vertices are type-II, (N V , N II ) = (N/2, N ).Its energy is This phase is 4-fold degenerate and breaks the translational symmetry.
Comparing the E cl Sq to E cl FF , we get the V = 2µ phase boundary between the fully flippable and square phase.Similarly, V = 0 is the classical boundary between the isolated and square phases.
The Rys F model corresponds to V = 0.It has two phases, the "anti-ferroelectric" for µ > 0 and the "disordered" for µ < 0 [4,44,45].The former matches the fully flippable phase and the latter the phase boundary between the isolated and the square phases (thick red line in Fig. 7), with a subextensive degeneracy W disordered = 4 L in N = L 2 clusters.The disordered manifold consists of all the configurations with only type-II vertices, including the square and the isolated ones.All arrows on a horizontal or vertical line point in the same direction in these configurations, the directions on different lines do not correlate.
The classical phases, particularly square one, motivate a quadripartite division of the lattice.Selecting the position of the flippable plaquettes with a counterclockwise direction of arrows, the four square states define the A, B, C, and D kind of plaquettes, see Fig. 9.This partition of the plaquettes allows us to write the configurations in the fully flippable and square phases as They may serve as variational wave functions and initial states for the Lánczos algorithm when we study the system with exact diagonalization in Sec.VI.
So far, we discussed systems with periodic boundary conditions.Extending the results for open boundary conditions or infinite systems size requires further discussion.For example, Ref. [52] considers a system with domain wall boundary conditions that fix a flux sector.Quantum dynamics then splits the Hilbert space into smaller fragments within the selected flux sector (Krylov spaces).They give the number of these and show that an RK-like exact eigenstate exists in each fragmented space, among others.

IV. SYMMETRY GROUPS AND ORDER PARAMETERS
We now turn to the t ̸ = 0 case, where the off-diagonal terms of quantum origin appear in the Hamiltonian.Previous studies determined the phase diagram as a function of V /t for µ = 0 [18,37,38,41].We apply numerical and analytical approaches to extend the phase diagram with the µ/t axis.But to identify the different phases in the (0, 0) flux sector, we need to know the symmetries they break and the respective order parameters.In this section, we systematically construct the order parameters from symmetry considerations using the mathematical tools of group theory.We conclude this section by formulating the Landau free energy and discussing the order of phase transitions.

A. Symmetry groups of the different phases
For convenience, we work in the packed loop representation below.All of the three phases -the fully flippable, the square, and the plaquette-are invariant to the translations by the t 2x = (0, 2) and t 2y = (2, 0) lattice vectors, and the mirror symmetries σx and σy with vertical and horizontal axes that split the squares into half, see Fig. 10(a).The order of the group T 2 formed by translations t 2x and t 2y is N/4 in a cluster with periodic boundary conditions respecting the division into four sublattices.The two orthogonal reflections σx and σy generate a point group isomorphic to D 2 with four elements.So the symmetry group G = D 2 × T 2 that preserves any of these three phases has | G| = N elements in the N = L 2 and N = 2L 2 type clusters.
On the other hand, the Hamiltonian commutes with all the elements of the wallpaper group G of the square lattice, which is p4m in the IUCr notation.The phases mentioned above break the symmetries of the quotient group The σ ′ x+y and σ ′ x−y reflections generate the symmetry group of the plaquette states, Since there is no subgroup relation between the D 2 and D ′ 2 symmetry groups, the phase transition between the fully flippable and the plaquette phase is first order according to Landau's criterium.However, the symmetry groups of the square phases are both subgroups of the D ′ 2 group of the plaquette phase (see Fig. 10(c)), so the transition between the plaquette and the square phases can be continuous.

B. Order parameters
To characterize the different phases, we construct order parameters below using the transformation properties of the vertices and characters, tabulated in Tab.I.It gives the transformation properties of the various ordered phases, from which we can calculate the characters and their irreducible representations.The fully flippable phase belongs to A 1 ⊕ B 1 representation, the square phase to A 1 ⊕ B 2 ⊕ E, and the plaquette as A 1 ⊕ B 2 .To distinguish them, we need to construct order parameters that transform according to the B 1 , B 2 , and E irreducible representations.Using the vertex operators, we find the following irreducible representations at a site i The ôFF i transforms as the B 1 , the ôPl i as the B 2 , and ôSq i as the two dimensional E irreducible representation of the D 4 .Here, the n i ( ) = | ⟩⟨ | is 1 if the site i contains vertex and zero otherwise, and similarly for other vertices.
Since the phases are invariant under G, the order parameters Ô transform as the trivial irreducible representation of G and are given as the sum over the elements of G acting on the local ôi operators, Ô = g∈ G gô i .Performing the sum, we TABLE I.The action of the D4 group elements, represented by permutations, on the fully flippable, plaquette, and square states, and the six possible vertices.The bottom row shows the action of the charge conjugation operator C. The definition of the fully flippable and the square states in the t → 0 limit is provided by Figs. 8 and 9.In the case of the plaquette phase, AD denotes the quantum state in which the A and D plaquettes resonate, and BC, where the B and C resonate.The last five rows are the characters of the irreducible representations.We grouped the rows belonging to the same conjugacy class.
construct the order parameters of the various phases as The values the order parameters take in different phases are summarized in Tab.II.Let us mention that Furthermore, the number operators, are invariant under D 4 (they transform as the A 1 irreducible representations).Eqs. ( 13) and ( 14) define the order parameters using vertices, unlike the quantum dimer model, where the order parameters depend on the occupation on bonds [53].We may ask ourselves why cannot we follow the same construction.To this end, let us denote by n l the occupation of the bonds indexed by l in Fig. 10(b); it is one if occupied by a loop segment and -1 if not.The operators n 1 ( ), n 2 ( ) ,n 3 ( ), and n 4 ( ) belong to the A 1 ⊕B 1 ⊕E representation.We may write the local order parameters for the fully flippable and square phases as However, the plaquette order parameter cannot be expressed as a linear operator in bond densities since we cannot combine n l s to transform according to the B 2 irreducible representation.The vertex operators in Eq. ( 13) are bilinear in bond occupation; they span a larger operator space, the n i ( ) and n i ( ) transforms as A 1 ⊕ B 1 and the n i ( ),n i ( ), n i ( ), and n i ( ) as A 1 ⊕ B 2 ⊕ E. The construction of the plaquette order parameter requires B 2 .Let us mention that the staggered flippability, which is bilinear in bond occupations, is also an obvious choice for a plaquette order parameter [19,38].Height representation is yet another tool to construct order parameters [37,41].
To complete the analysis, the charge conjugation acts on the order parameters as Only the plaquette states are invariant to charge conjugation, with a symmetry group enlarged to having eight elements.
The fully flippable phase breaks C but is invariant under the fourfold rotation combined with the charge conjugation, CC 4 .The symmetry group of this phase is then The square phase also breaks charge conjugation, but they are invariant to the combination of the C and a reflection.For example, the symmetry group C s of the SqA and SqD states (see Fig. 10(c)) is extended to The quotient of the groups defined in Eqs. ( 19) and ( 21) is isomorphic to D 2 /C s ∼ = C 2 .Even though we added the charge conjugation, a single generator remains broken at the phase transition from the plaquette phase to the square phase, preserving the possibility of a continuous transition.We will construct and analyze the Landau free energy in the next subsection to see how this happens.

C. Landau free energy
Once we identified the order parameters and their transformation properties, we can write down the free energy invariant under the D 4 quotient group.Including up to quartic terms, its form is The coefficients c are some functions of the couplings.The non-geometric charge conjugation symmetry C further restricts the allowed terms.For instance, since the O FF is odd under C [see Eq. (18a)], the inclusion of the C removes the Sq,1 term from the Landau free energy F. We are then left with where we introduced the parametrization of the square order parameter.The appearance of the angle ϕ in the cubic term, together with the plaquette order parameter, is the consequence of the hierarchy of symmetry breaking presented in Fig. 10.For example, we can develop the SqA and SqD states from the plaquette phase resonating on the AD sublattice by breaking one of the reflections.In the Landau functional language, sin 2ϕ will be fixed to 1 throughout the phase transition between the AD plaquette and the SqA and SqD phases when the square order parameter becomes nonzero.Similarly, breaking the BC plaquette phase into SqB and SqC sets sin 2ϕ = −1.This agrees with ϕ = π/4 + nπ/2 (n ∈ Z) corresponding to the four classical square states.Moreover, it implies that c 0,1,2 < 0.
Because of the many undefined coefficients, it is not easy to describe the phase diagram and the order of the phase transitions.We only mention that the Landau free energy allows both first and second-order transition for the phase boundary between the square and the plaquette phase, depending on the sign of the quartic term.In Appendix C, we present a simple variational wave function to describe the plaquette-square transition.It displays both a first and a continuous boundary separated by a tricritical point.
Let us note that the charge conjugation symmetry is particular for the 6-vertex model, as half of the bonds are occupied in the fully packed loop representation.Hence, the charge conjugation symmetry is absent in the quantum dimer model.

V. THE MANIFOLD OF THE ISOLATED STATES
Configurations without flippable plaquettes form the isolated manifold.They are disconnected from other configurations by local flips and consist only of type-II vertices if periodic boundary conditions are imposed.The quantum term H t [Eq.(1a)] annihilates any of them, so their energy remains the classical one, E Iso = µN .
Ref. [18] provides a recipe to construct all the configurations in the isolated manifold.If all the arrows along the horizontal or vertical lines in a configuration point in the same direction, with the proviso that either the horizontal or the vertical lines must be oriented alike, none of the plaquettes is flippable.See Fig. 11(a) and (c) for examples.It implies that at least one of the flux indices of these isolated configurations has to be extremal.
This recipe helps to determine the degeneracy of this manifold.Let us consider the N = L 2 cluster as an example.We then have two choices of fixed orientation, both in the case of horizontal and vertical lines, resulting in a factor of four (Fig. 11(a) illustrates one of these four states).In the nonfixed direction, each of the L lines can point in two directions, giving 2 L possibilities altogether.Supposing that l lines point in one and L − l in the other direction defines configurations in the (±L, L − 2l) or (±L, L − 2l) flux sectors, with degen- eracies L l .Considering the double counting of (±L, ±L) flux sectors, we end up with 4 × 2 L − 4 isolated states in the manifold.The degeneracy exponentially grows with the linear size of the system, so it is subextensive.We can apply similar considerations to clusters with other geometries.
Knowing their energies allows us to determine an exact region in the phase diagram where they are the ground states.For this purpose, we will use Gerschgorin's theorem, which claims that |H ii − ε i | ≤ j̸ =i |H ij | for a finite cluster Hamiltonian H with eigenvalues {ε i }.We use a basis where both NV and NII operators are diagonal, and the quantum term H t is strictly off-diagonal.Then for the ith configuration having N V flippable plaquettes and N II type-II vertices the diagonal term is H ii = V N V + µN II .The sum over the off-diagonal terms j̸ =i |H ij | = N V t, as the configuration connects to exactly N V other ones, each with −t amplitude (note that we chose t > 0).So we can write Resolving the absolute value, the ε i s become bounded as Let us denote by ∆ i = ε i − µN the energy gap between the eigenvalue ε i and the energy of the isolated manifold.Then Isolated states are ground states while the gap ∆ i ≥ 0, which is satisfied when (V −t)N V −µ(N −N II ) ≥ 0. We need to find the region in the parameter space of V , µ, and t, where this inequality holds, provided that Eqs. ( 5) constrain the possible N V and N II values into the triangle shown in Fig. 6.It is a simple linear optimization problem analogous to finding the classical phase diagram with V replaced by V − t.The extrema occur at the corners of the triangle, where (N V , N II ) = (0, N ) (isolated phase), (N/2, N ) (square phase), or (N, 0) (fully flippable phase).Eventually, we conclude that the isolated configurations form the ground state manifold when Fig. 11(f) shows this region.The approach above does not tell us whether the isolated states remain the ground states outside this region.Below we will consider excitations that become soft at the positive µ = V − t boundary of the region above, promoting it to a phase boundary.Following Ref. [18], we can calculate the energy of a leapfrog excitation (we call it 'frogton'), depicted in Fig. 11(d).It exists in a cluster of a suitable skew shape, defined by the g 1 = (L, 2) and g 2 = (0, L) lattice vectors, for example.The boundary condition introduces a correlated kink-antikink pair on the line along which the arrows are reversed [the magenta line in Fig. 11(d)].A frogton comprises two flippable plaquettes and two type-I vertices, so the diagonal energy is 2(V − µ).Flipping one of the two plaquettes, the excitation hops with an amplitude −t and acquires a bandwidth of 4t centered at the diagonal energy.The minimal energy of this variational state gives an upper bound for the gap, ∆ = 2(V − µ − t).Combined with Eq. ( 28), the gap closing along the line for µ ≥ 0 in the phase diagram gives the exact boundary of the isolated phase.
In a square-shaped cluster with N = L 2 geometry, the flux sector that closes the gap originates from the "staircase-like" excitation shown in Fig. 11(e) in the m = (L − 2, L − 2) flux sector.Unlike frogtons, where the kink-antikink pair is confined within a lattice spacing, in this cluster, the number of the kinks and antikinks is not conserved, making the calculations more complicated.Exact diagonalizations up to L = 10 in N = L 2 clusters reveal that the finite size gap between the ground state energy in the (L − 2, L − 2) flux sector and the energy of the isolated manifold exponentially decreases with the system size along the V = t + µ line, as demonstrated in Fig. 12. Therefore we may conclude that the gap closes along the whole V = t + µ line in the thermodynamic limit.We will further scrutinize this question in Sec.VIII D.

VI. PHASE DIAGRAM FROM EXACT DIAGONALIZATION
The number of the ice rule configurations grows exponentially with the system size N , restricting the size of the clusters that can be studied numerically.Using the conservation of the fluxes (m x , m y ), we diagonalized 16, 32, and 36 site clusters having the full D 4 point group symmetry of the square lattice.We used the standard Lánczos algorithm to calculate the energy of the ground state and low-lying excitations.
To begin, we scanned the V /t and µ/t parameter space to reveal which flux sectors give the ground state in the N = 16 and 32 site clusters.Fig. 13  (d) (d) 15.The same as Fig. 14 for N = 36.The positions of the phase transitions are unchanged compared to the 32-site cluster.We also show all the states in the selected energy window, indicating their C parity.
from the m = (0, 0) and the isolated manifold for positive µ values.Computations on 36 site clusters also confirmed the existence of this phase.In this section, we discuss the phases in the (0,0) and isolated sectors, and we will present results about the grey region in Sec.VIII.Identification of the isolated phase is numerically straightforward.Since we know the exact energy of the isolated states, E Iso = µN , it is enough to calculate the lowest energy levels of the non-isolated flux sectors with ED and compare them to E Iso .Distinguishing the different phases in the (0,0) flux sector is more challenging.We applied two methods, (i) one based on the order parameters and (ii) one based on energy level spectroscopy [18,37].
Re (i), we computed the expectation values of the squares of the order parameters given in Eqs.(14), ⟨GS| Ô2 FF |GS⟩, ⟨GS| Ô2 Plaq |GS⟩, and ⟨GS| ÔSq • ÔSq |GS⟩, by calculating the ground state wave function |GS⟩ numerically.We show the result along the V = 0 line in Figs.14(b) and 15(b) for the 32 and 36 site clusters.Since the squares of the order parameters transform as the A 1 irreducible representation in Tab.I and contribute to the energy [c.f. the Landau functional in Eq. ( 23)], we estimate the phase boundaries by inflection points of their expectation values.For this purpose, we calculate the extrema of the ∂⟨ Ô2 ⟩/∂(µ/t) in Figs.14(c) and 15(c).According to this criterium, the square phase is realized for µ/t ≲ −0.78, the plaquette phase between −0.78 ≲ µ/t ≲ 0.3 − 0.4, and the fully flippable phase for 0.3 − 0.4 ≲ µ/t (the reason behind the uncertainty for the upper boundary is that the inflection points of the ⟨GS| Ô2 FF |GS⟩ and ⟨GS| Ô2 Plaq |GS⟩ are not at the same µ/t values).We note that the plaquette order parameter is also finite in the square phase, as expected from the discussion of the hierarchy of the symmetry breaking in Sec.IV B.
Re (ii), energy level spectroscopy provides a less direct but more powerful tool to estimate the phase boundaries in a finite-size cluster [54].Since the Hamiltonian preserves all the symmetries of the model, its wave functions transform as the irreducible representations of the underlying group.We can identify a set of low-lying states in the spectra, belonging to specific irreducible representations of the D 4 quotient group.They collapse into a degenerate ground state manifold in the thermodynamic limit.The symmetry-breaking states are linear combinations of the wave functions in this manifold.Characterization of these low-lying states gives the basis for level spectroscopy, a tool for detecting the various phases.
As an example, let us consider the fully flippable phase.The classical |FF1⟩ and |FF2⟩ configurations manifestly break the σ ′ x+y reflection symmetry, but we may linearly combine them into the For a finite cluster, the energy of the (0, 0) e level is different from the (π, π) o level, but we expect the gap between them to vanish in the thermodynamic limit when the fully flippable phase is realized.We can associate the linear combinations above with the irreducible representations listed in Tab.I.More generally, the (0, 0) e belongs to the trivial A 1,e irreducible representation, the (π, π) e transforms as B 2,e , the (π, π) o transforms as the B 1,o , and the (π, 0) o and (0, π) o span the two-dimensional E o .Since the charge conjugation C commutes with all the point group elements in D 4 , the irreducible representation of the {1, C} × D 4 group are simply the irreducible representation of the D 4 appended with the even (e) or odd (o) parity with respect to C. We note that one shall be careful with the interpretation of the momentum labels, as they can be different for the arrow and fully packed loop representation (for example, the fully flippable state is translationally invariant in the fully packed loop representation but not in the arrow representation).Above, we used the arrow representation.
We prepare states with appropriate momentum and parity and use them as input to the Lánczos code since the iterations in the method preserve their symmetry.Figs.14(d) and 15(d) show the momentum and parity resolved spectra for the most important (0, 0) e , (π, π) e , (π, π) o , and (π, 0) o symmetry sectors [(π, 0) o and (0, π) o are degenerate].The ground state is always in the (0, 0) e sector.For µ/t ≳ 0.29, the first excited state is in the (π, π) o sector, just like in the decomposition (30b), and it becomes degenerate with the ground state as µ/t increases and approaches the classical limit.For µ/t ≲ 0.29, the first exited state is in the (π, π) e sector, and a level crossing occurs at µ/t ≈ 0.29.Since the phases break discrete symmetries, their spectrum is gapped.The low-lying excitations belonging to different irreducible representations on the two sides of the level crossing lead to different symmetry breaking and, thus, phases in the thermodynamic limit.We identify the phase for µ/t ≳ 0.29 as fully flippable.Ideally, a finitesize scaling should be performed to accurately determine the phase boundary, as in Ref. [18] for the case µ = 0.However, the positions of the level crossings in 1/N are not monotonic for general values of the parameters.Nevertheless, our value agrees well with the result of quantum Monte Carlo simulation [19], t/µ ≈ 3.47 (i.e., µ/t ≈ 0.288).
Next, consider the limit for negative µ, where the square phase appears.The combinations of the four classical square configurations give the following momentum and parity eigenstates: In Figs.14(d) and 15(d), four states with precisely these momenta and parities are quasi-degenerate for µ/t ≲ −0.78, supporting the realization of the square phase in this region.As µ increases from µ/t ≲ −0.78, the energy levels of the μ/t V/t approximate wave functions.The |PlAD⟩ − |PlBC⟩ is in the (π, π) e symmetry sector.Instead of a level crossing, an avoided level crossing characterizes the square-plaquette phase transition.
To further elucidate the nature of the excitation spectrum, we calculated the dynamical correlation functions of the order parameters, using the Lánczos method, where X are excited states and α = FF, Sq, Pl denotes the order parameter.We first calculated the ground state |GS⟩ for a given value of parameters µ and V , applied the operator Ôα to the |GS⟩, and then used Ôα |GS⟩ as the initial state for the second run of the Lánczos procedure.The algorithm then computes the matrix elements in the definition of S α (ω).Not surprisingly, the largest matrix elements are for the lowest-lying excitations of the momenta and parities corresponding to the symmetry of the phase.In Figs.14(d) and 15(d), the size of the open symbols is proportional to the values of the matrix elements.
Using the above criteria, we established the phase diagram in the parameter space of V /t and µ/t, Fig. 16.We determined the first-order boundaries between the fully flippable and the plaquette and between the fully flippable and the square phases by following the positions of level crossings and the boundary between the square and the plaquette phase following the positions of the inflection points.We also checked that the parameter values of the avoided level crossings coincide with those of the inflection points.We obtained the boundary to the isolated state manifold by comparing their energy to the ground state energies in the (0, 0) flux sectors.
The main consequence of the quantum fluctuations is the appearance of the plaquette phase with resonating alternating plaquettes that fills up the central region of the phase diagram.The plaquette phase extends along the V ≈ µ line to larger positive values of V , following the unknown phase with finite flux sector ground states that separate it from the isolated states.Otherwise, for large |V | and |µ| values, the phase diagram is consistent with the classical one, shown in Fig. 7.In the next section, we will confirm the validity of some of these phase boundaries using perturbation theory.

VII. PHASE BOUNDARIES FROM THE PERTURBATION THEORY
We use the size-consistent Rayleigh-Schrödinger perturbation theory below to estimate the ground state energies in the fully flippable and square phases.We calculate the secondand fourth-order corrections in H t to the ground state energy of the configurations drawn in Figs. 8 and 9 in the t → 0 limit.We get an estimate for the first-order phase boundary between these phases by comparing their energies.Furthermore, comparing the energy of the square phase to E Iso = µN will provide the corresponding phase boundary.
The perturbation series calculation is straightforward for the square and fully flippable phases in any representation in which the H cl is diagonal.We give details in Appendix D. The ground state energies up to the fourth order in the flip-ping amplitude t are Solving the E FF = E Sq , we get the following Padé approximants for the phase boundary between the fully flippable and the square phases, All three Padé approximants equally satisfy the energy equation up to the fourth order.The orders [m/n] of the Padé approximants above denote the power of t in the numerator (m) and the denominator (n).
Similarly, the E Sq = E Iso equation provides the phase boundary between the square and the isolated phases.The Padé approximants are Padé[4/0], Fig. 17 shows these approximants together with the numerical results of the ED calculation for both phase boundaries.The comparison of different orders of Padé approximants allows us to estimate the convergence of the perturbation series: the different lines are essentially superimposed on each other in the relevant domains, indicating a rapid convergence of the series.The perturbation expansion also agrees well -typically within two decimal places -with the phase bounds extracted from ED calculations on finite clusters.

VIII. THE ROKHSAR-KIVELSON POINT AND THE LIQUID PHASE
The exact diagonalization of the 32-site cluster shows that the flux sector of the ground state gradually increases from the m = (0, 0) in the plaquette phase as we approach the manifold of isolated states (the dark grey area in Fig. 13).To gain a deeper insight into the properties of this phase, wedged between the plaquette phase and the isolated states and emanating from the quantum critical RK point, we perform a Monte Carlo evaluation of the RK wave function [20,55].This enables the extension of cluster sizes to up to 576 sites close to the RK point.

A. First order perturbation around the Rokshar-Kivelson point
The RK point is a particular point in the phase diagram since the exact ground state wave function |RK(m)⟩ is known: it is the equal amplitude superposition of the configurations in an N V diagonal basis within the flux sector m [46].The ground state energy is the same in all the flux sectors.But the expectation values of the flippable plaquettes NV and NII operators depend on m.We use the Hellman-Feynman theorem at the RK point to estimate the splitting of the ground states.In the first order, we approximate the lowest energy of the perturbed RK Hamiltonian in a given flux sector with the formula The denote the expectation values of the number of flippable plaquettes and type-II vertices.Comparing these energies, we can figure out the flux sector of the ground state (a similar argument appeared in Ref. [56] for the quantum dimer model on the triangular lattice and in Ref. [57] for a quantum dimer model on the honeycomb lattice).

B. The Monte Carlo method
The Monte Carlo method uses the RK wave function to evaluate the expectation values by random sampling [20].We started the simulation from a configuration formed by arrows directed along horizontal and vertical lines since this allowed the selection of the flux sector and generated new configurations by randomly flipping plaquettes.We discarded the first 5 million configurations to reach thermalization.Following thermalization, we measured N II and N V after every N step and updated their averages, where N is the system size.The number of elementary steps in a Monte Carlo run was typically between 5 × 10 7 and 10 9 flips, depending on the system size and the statistical error.After we exported the averages and continued measuring another four times.Repeating the procedure above five times, we collected N MC = 25 average value pairs for each flux sector.Denoting by m i (i = 1, . . ., N MC ) the averages (means) from the Monte Carlo runs, we estimate the statistical error by the standard error of the mean, given by the formula Here the is the mean value of the m i averages.

C. Dependence of expectation values on flux sectors
For Monte Carlo calculations, we used clusters with N = L 2 geometry up to N = 576 sites and N = 2L 2 geometry up to N = 512.The finite size dependence of the n V = ⟨N V ⟩/N and n II = ⟨N II ⟩/N densities is shown in Fig. 18, together with ED data for the N = 32 and N = 36 sites to check the consistency of the data.The finite size fitting functions are as follows.The density of the type-II vertices and flippable plaquettes for the N = L 2 class of clusters (red curves in Fig. 18) The same for the clusters with the N = 2L 2 geometry (green curves in Fig. 18) The expectation values for both cluster geometries extrapolate to the same values in the thermodynamic limit well within the error bars.
It is instructive to evaluate the ratio ⟨n II ⟩/⟨n I ⟩ in the RK wave function for different flux sectors.The density of the type-I vertices is ⟨n I ⟩ = 1 − ⟨n II ⟩, so for m = (0, 0) we get ⟨n II ⟩ (0,0) ⟨n I ⟩ (0,0) = 0.6199 ± 0.0001 0.3801 ± 0.0001 = 1.6308 ± 0.0006 .( Were the vertices uncorrelated, we would expect ⟨n II ⟩/⟨n I ⟩ = 2 instead of 1.63.This ratio improves as the flux increases, the limiting case being the isolated states with type-II vertices only. Inspired by the flux dependence of the energy in the U (1) liquid in 3D [22,23], we plotted how the expectation values depend on the flux sectors for different cluster sizes in Fig. 19.The plot reveals the linear dependence of ⟨n II ⟩ m on m 2 /N for not too large values of the flux m and the slope appears to be independent of the size of the cluster.An additional factor of 2 compensates for the geometry of the clusters.
To further elucidate the linear dependence on m 2 , we plotted the finite size dependence of the gaps divided by m 2 /m 2 0 in Fig. 20.The reason to divide m 2 by m 2 0 is the smallest nonzero flux (i.e., unit of flux).It is m 0 = (0, 2) with m 2 0 = 4 in N = L 2 cluster with even L.In the N = 2L 2 cluster, the unit of flux is m 0 = (2, 2) and, therefore, we divide m 2 by 8.The gaps in both cluster geometries are the same when threaded by the unit flux m 0 , independently of the system size.We collected the finite size scaled values of the ∆(m)m 2 0 /m 2 in Tab.III.The numbers extracted from the lowest flux sectors are identical within the error bars, and we may conclude that Putting together with the m = (0, 0) values in Eqs.(41), we get the behavior of the expectation values in the thermodynamic limit These are the solid black lines in Fig. 19.

D. The quantum electrodynamics of the RK wave function
How do we understand the scaling of the expectation values at the RK point?The local 2-in/2-out constraint at the vertices represents a divergence-free field.Associating the arrows with an electric field, one can build a kind of emergent quantum electrodynamics (QED) in the spin ice systems via the Gauss law, leading to a gapless U (1) spin liquid phase.It has been widely studied in 3D models within the context of quantum spin ice [20][21][22][23], and found that the U (1) liquid extends beyond the RK point.In 2D, the gapless liquid phase is usually at the RK point only.All flux sectors have the same energy, and as we leave the RK point, a gap opens in the ground state flux sector, with possible exceptions [57][58][59][60], as we will see later.The expectation values of the n V and n II do not follow the behavior of the energy.While the energy values are degenerate, the ⟨n V ⟩ and ⟨n II ⟩ values vary with the flux.The energy of the electric field is proportional to where the integral is over area A. We neglect the "magnetic" part of the emergent QED.
In an effective theory, the average electric field on the lattice E is proportional to the m.Let us reverse a single arrow to see how the electric field emerges.It creates two vertices that are neither type-I nor type-II: a vertex with a 3-in/1out (charge) and a vertex with a 1-in/3-out (anti-charge) arrows.These vertices can be considered as fractional charges, spinons, or monopoles, according to the actual physical problem we apply the Q6VM model.Moving one of these defects (charges) across the periodic boundary by reversing other ar-rows and eventually annihilating them makes a loop we considered e.g. in Fig. 4. It changes the flux sector and introduces a finite electric field E when the arrows are coarse-grained.The strength of the average field is proportional to the density of the flux lines, E = qm/L for the N = L 2 shape cluster and E = qm/2L for the N = 2L 2 shape cluster, where q is the charge of the monopole.Taking the area A = N , squaring the E and replacing it into Eq.( 46), we get for the energy At the RK point, the degeneracy implies ε = 0.
The expectation values of the ⟨n II ⟩ m and ⟨n V ⟩ m in Eq. ( 45) are also quadratic function of the electric field.To further corroborate this statement, in Fig. 21, we show the dependence of the expectation values of n II and n V for a fixed number of sites but changing the aspect ratio of the rectangular cluster.The gap varies according to our expectations, the flux density.
Since, in the vicinity of the RK point, the energy follows Eq. ( 37), comparing the flux-dependent part with Eq. ( 47), we get For energies above the small gap in the m = (0, 0) sector, we expect the system to follow the energy of the emergent QED with a tunable ε permittivity.This is in the spirit of Ref. [23], which considers the emergent QED in the 3D quantum spinice model.

E. The phase boundaries emanating from the RK point
To describe the small perturbations around the RK point, we parametrize the V − t and µ with the angle θ where δ is some small energy scale.Next, for a value of θ and system size, we calculate the energy in Eq. ( 37) for each flux sector m, and find for which it is minimal.The result of this energy minimization is presented in Fig. 22(a) for the full circle around the RK point.We recovered the phase boundaries anticipated from the ED calculations: the first-order phase transition between the isolated states and the plaquette phase for negative values of µ and the liquid phase for µ > 0.
To determine the boundary between the isolated states and the plaquette phase more precisely, we compare the energy density µ of the isolated states with the energy of the m = (0, 0) sector using Eq.(37), that leads to the following equation: For the angle θ 1st we then get in the thermodynamic limit, using the extrapolations given in Eqs.(41), so θ 1st = (1.26427±0.00006)π,taking into account that both V − t and µ are negative at this boundary.Let us note that in the denominator the 1 − ⟨n II ⟩ (0,0) = ⟨n I ⟩ (0,0) , the density of the type-I vertices appears.In Fig. 22(b), we zoom in on the tiny region where the liquid phase appears.The isolated manifold gives the ground state up the θ = π/4, in full agreement with Eq. (28) in Sec.V.The flux sector first appears next to isolated manifold is the m = (L − 2, L − 2) in the clusters having N = L 2 sites.
To get the boundary between the plaquette and the liquid phases, we compare the energies of the m = (0, 0) and small m flux sectors.This involve the gaps ∆ V (m) and ∆ II (m) and provide the condition so that in the thermodynamic limit for values of m tabulated in Tab.III.This translates to θ 2nd = (0.2912 ± 0.0007)π, indicated by the tiny arrow in Fig. 22(b).The window for the liquid state is thus tiny, about θ 2nd − π/4 ≈ 0.041π.We note that the value of θ 2nd is where the permittivity ε in Eq. ( 49) changes sign.
It is difficult to resolve the precise character of the liquid phase.It is unclear whether the topological sectors increase continuously or whether we are faced with an infinite sequence of incommensurate states, exemplifying a "devil's staircase" (also called "Cantor deconfinement") [57][58][59][60].Possible evidence for the latter scenario is the plateau at half maximum flux, m = (L/2, L/2).Analysis of finite-size scaling suggests a finite width of the plateau, which is about 5% of the size of the liquid phase, see Fig. 22(c).It is adjacent to the m = (L/2 − 2, L/2 + 2) and m = (L/2, L/2 − 2) flux sectors; however, this does not follow assuming a perfect m 2 dependence of the expectation values on the flux.

IX. STRUCTURE FACTORS
In this section, we determine the zero-temperature correlation functions and the structure factors using exact diagonalization.We will first discuss the structure factor in the fully packed loop representation and then the magnetic structure factor in the arrow representation.The latter will allow us to compare our results to the ones observed in the artificial spin ice with superconducting flux qubits by King et al. in Ref. [42].

A. Correlations in fully packed loop representation
We define the correlation function with respect to a horizontal bond as where |GS⟩ is the translationally invariant ground state.n r measures whether the bond centered at r is occupied (n r = 1) or not (n r = −1).Since the |GS⟩ transforms according to the trivial irreducible representation, the ⟨GS|n r |GS⟩ = 0 and we can use the above definition of C h (x, y).The x and y values are either both integers or half-odd integers.We calculate the ground state wave function |GS⟩ using the Lánczos algorithm in N = 16, 32, and 36 site clusters with periodic boundary conditions for a few selected parameters, representing the different phases.The exact diagonalization provides a fully symmetric |GS⟩ in a finite cluster with periodic boundary conditions (the |GS⟩ is in the (0, 0) e symmetry sector, see Figs. 14 and 15).For practical purposes, we introduce the where the coordinates (i x , i y ) are integers.The C hh (i x , i y ) is a correlation function between horizontal bonds, C vv (i x , i y ) between vertical bonds, and C vh (i x , i y ) and C hv (i x , i y ) between orthogonal bonds.They provide sufficient information to obtain both the density and the magnetic correlation function in real and reciprocal space.How they behave under the action of the point group symmetries is described in Appendix E. Fig. 23(a) displays the bond-bond correlation function in the ordered square, plaquette, and fully flippable phase, as well as for the quantum-disordered RK point from ED calculations on the 36-site cluster, and for the disordered phase boundary in the classical phase diagram.While the bond-bond correlations decay rapidly in the plaquette phase and at the RK point, the long-range pattern of the ordered loops manifests itself in the classical square and fully flippable phase.In the square phase, when the central horizontal bond at (i x , i y ) = (0, 0) in Fig. 23(a) is occupied (blue disk), all the horizontal bonds in the same columns are also occupied and the next column of horizontal bonds is empty (red disks), in full accordance with Fig. 9(b).For the fully flippable phase, the occupied horizontal bond invokes the occupation of the horizontal bonds (all the horizontal bonds are blue, and all the vertical bonds are red), c.f. Fig. 8(b).In the disordered manifold of the classical Rys-F model (the boundary between the isolated and square phase in Fig. 7), the bond-bond correlations are finite only along a line (otherwise, the average over the disorder nulls the correlations).
The Fourier transform of the real space correlation function is the structure factor where both r and r ′ in the sum run over the N horizontal and N vertical bonds.Separating the vertical and horizontal bonds, we arrive to the expression.Using the symmetry properties described by Eqs.(E1)-(E4), one can show that the S(q) is real and satisfies the full D 4 point group symmetry in the q space.Since the centers of the horizontal and vertical bonds form a square lattice rotated by 45 • and lattice constant 1/ √ 2, the S(q) is periodic in the reciprocal space for momenta shifts by (2π, 2π) and (2π, −2π): S(q) = S(q + (2π, 2π)) and Figure 24 displays the spin structure factor S(q) for the typical examples shown in Fig. 23(b).As for S(q), Bragg peaks dominate the spin structure factor in the square and the fully flippable phase, and diffuse scattering is seen for the plaquette phase and the RK-point.In the case of the RK point, the pinch points move to Q values that are now multiples of 2π.The sub-divergent lines of the disordered manifold get a momentum-dependent intensity from the dipolar factor in Eq. ( 63), in addition to the momentum shift.We can now compare our zero-temperature results for the structure factors with those observed in the artificial quantum spin-ice experiment.Our magnetic structure factor S(q) should correspond to that in the "strong coupling" regime in Fig. 2 in Ref. 42 (J = J MAX in their notation), where the frequencies (densities) of type-III and IV vertices are negligible.We shall bear in mind that the experiment dealt with a system with open boundary conditions at finite temperature while we work with periodic boundary conditions at zero temperature.
Let us start with the topmost row in Fig. 2 in Ref. 42, the "degenerate ice" with J ⊥ /J ∥ = 1.At this point, the observed N I : N II = 1 : 2 ratio between the number of type-I and II vertices corresponds to their statistical weights, as the local Hilbert space on a site consists of two type-I and four type-II vertices, so neither of the vertices is favored in this case.The same ratio also occurs at the RK-point, where the ground state wave function is an equal-amplitude superposition of the allowed vertex configurations.We find excellent agreement with our Fig.24(d) regarding the overall distribution of the weight.However, there is one notable difference: the pinch points at the Q = (2n x π, 2n y π) show no weights in our Fig.24(d), while in Ref. 42 they are finite.A possible resolution of this discrepancy is that the weight at the center of the pinch point comes from finite flux sectors.The ⟨n II ⟩ (0,0) /⟨n I ⟩ (0,0) ≈ 1.631 in Eq. ( 42) also supports this idea, to recover the measured ⟨n II ⟩/⟨n I ⟩ ≈ 2 we shall take into account other flux sectors as well.However, it is also possible (c) V=0, μ/t=0.5 Fully flippable Type-II disordered q x /2π 0 S(q) FIG.24.The magnetic structure factor S(q) in the arrow representation, given by Eq. ( 63) and calculated by exact diagonalization on the 36-site cluster for (a)-(d) and analytically for (e).The ordered phases are (a) the square phase for µ/t = −1, (b) the plaquette phase for µ = 0, and (c) the fully flippable for µ/t = 0.5, keeping V = 0 in all three cases.We find divergent Bragg peaks in the square and the fully flippable phases, while S(q) remains diffuse in the plaquette phase.(d) S(q) at the RK-point (V = t and µ = 0) displays the pinch points at Q = (2πzx, 2πzy), where zx, zy ∈ Z. (e) The magnetic structure factor of the V = t = 0 classical disordered manifold.The amplitudes of the sub-divergent horizontal and vertical lines are woked out in Appendix F 4. The green dotted square denotes the boundary of the extended Brillouin zone, the same as in Fig. 23(b).
that the temperature was so high that it concealed the quantum fluctuations and the complete disorder is a consequence of the almost identical Boltzmann-weights of ice rule configurations.Ref. 62 studied this case with classical Monte Carlo and noted the same structure factor.Let us proceed to the second row in Fig. 2 in Ref. 42, denoted by "Type-I bias".The choice of parameters resulted in slightly more type-I vertices than type-II vertices, with a ratio N I : N II ≈ 0.54 : 0.46.According to our ED calculations shown in Figs.14(a) and 15(a), the experimental N I : N II ratio corresponds to the µ/t ≈ 0.12 in the plaquette phase when V = 0, close to the fully flippable boundary at µ/t ≈ 0.29 (we find ⟨N I ⟩ : ⟨N II ⟩ = 0.6 : 0.4 in our calculations at the phase boundary).In the experimental structure factor, we can identify both the Bragg peaks of the fully flippable phase with a diffuse scattering around, just like in our Fig.24(c), and the finite rhombi-like structures at (q x , q y ) = (0, ±3π) and (±3π, 0) typical for the plaquette phase shown in Fig. 24(b).Therefore, it is hard to identify the ground state unambiguously.While the fully flippable phase is classical, the emergence of the diffusive structures compatible with the plaquette phase may indicate quantum effects at play.Lastly, we turn to the third row in Fig. 2 in Ref. 42, which shows the case where N I : N II ≈ 0.19 : 0.81 in the strong-coupling limit ("Type-II bias" with J ⊥ /J ∥ = 0.98).We can reproduce the reported ratio of occupancies setting µ/t ≈ −0.9 in our V = 0 exact diagonalization calculation, which is in the square phase but close to the plaquette phase.The magnetic structure factors are quite different: our calculation reveals apparent Bragg peaks with some diffuse scattering, while the experimental plot shows strong intensities along straight lines.But such structures appear in the finite temperature plots in Fig. 3(b) in Ref. 62 and our calculation of the magnetic structure factor in the V = t = 0 disordered manifold.It implies that the thermal fluctuations destroy the quantum-mechanical order, and the observed state is a mixture of different flux sectors.The fact that quantum features are revealed in some cases and not in others may be related to the excitation energies in different phases, and how they compare to the temperature.However, this is only a hypothesis at this stage and would require further study.
Above, we compared the magnetic structure factors S(q).The equivalent of the S(q), defined by Eq. ( 58), is shown in Fig. S4 in the Appendix of Ref. 42 (we note that the coordinate axes in Fig. S4 are rotated by 45 • in contrast to the magnetic structure factors S(q) and the indicated Brillouin zone corresponds to our extended Brillouin zone.).

X. SUMMARY
We studied the ground state properties of a quantum six vertex model on the square lattice that distinguishes the type-I and type-II vertices.We established the zero-temperature phase diagram and the static correlation functions in real and momentum space using analytical and numerical methods.
Regarding the classical (t = 0) model, we found three extended phases in the µ-V parameter space.As discussed in sections II-III, the twofold degenerate fully flippable phase contains only type-I vertices.The fourfold degenerate square phase and the sub-extensive manifold of isolated configurations consist of only type-II vertices.Figure 7 summarizes the classical phase diagram.The fully flippable phase is the analog of the "antiferroelectric" phase in the Rys F model [43], while the subextensive boundary between the square phase and the isolated manifold is equivalent to the "disordered" phase identified in Refs.[44,45].
Classification of the configurations by the number of flippable plaquettes and type-II vertices revealed that the two numbers are correlated and form a triangle presented in Fig. 6.Configurations at the corners of the triangle define the three classical phases mentioned above.In Sec.IV, we identified the symmetries broken in these phases, constructed diagonal order parameters as irreducible representations of the D 4 point group using vertices, and wrote down the Landau free energy for the model taking into account the charge conjugation symmetry.
To derive the phase diagram for the quantum model, we applied the Lánczos method to diagonalize the model on finitesize clusters (N ≤ 36) with periodic boundary conditions.We calculated the expectation values of the order parameters and followed the level crossings in the low-energy spectra in a wide range of parameters (see Figs. 14 and 15 for V = 0 case).Figure 16 summarizes our phase diagram.Interestingly, the boundary of the plaquette phase does not confine to small V /t and µ/t but extends along the line V = µ + t for large values, together with a gapless liquid phase emanating from the quantum critical Rokhsar-Kivelson point.The extension of the quantum behavior results from the highly degenerate V = µ classical boundary so that t competes with V − µ.A rapidly converging perturbation expansion in t up to the fourth order confirmed the accuracy of the phase boundaries.Beyond the numerical methods, we applied Gerschgorin's theorem to obtain the exact boundary between the isolated and the liquid phases in Sec.V. A simple variational treatment indicated the existence of a tricritical point between the square and plaquette phases, allowed by the Landau-free energy expansion.
We applied the Hellman-Feynman theorem to reveal the nature of the liquid phase and the splitting of the ground state degeneracy of the multicritical Rokhsar-Kivelson point in Sec.VIII.Using Monte Carlo simulations for clusters up to 576 sites, we calculated the expectation values of the densities of flippable plaquettes and the type-II vertices scanning through the flux sectors.We obtained the phase boundaries emanating from the quantum critical point.The liquid phase unveils itself as a possible manifestation of the "devil's staircase" [58,59], with evidence for a finite-width plateau at half the maximum flux sector.The results also allowed us to study emergent quantum electrodynamics and to determine the electrical permittivity near the RK point.
Finally, in Sec.IX, we presented the zero-temperature structure factors for the various ordered phases and the RK point.We compared our results with the experiment on the artificial spin-ice system formed by superconducting qubits of Ref. [42].While many features of our calculation and the measurement agreed, there were also some that we could not interpret using our zero-temperature calculation.TABLE IV.The degeneracy (deg.) of the (0, 0) flux sector for different system sizes N .The third column is the calculation based on Lieb's formula, and the ratio with the actual degeneracy.The fifth and sixth columns are the estimate from Pauling's formula and the ratio.While the ratio is stable for Lieb's estimate, the Pauling formula underestimates the number of configurations.case of equality, so the estimate is sharp.2N V + N II ≤ 2N : Out of four plaquettes around a type-II vertex, at least two are non-flippable.Let us fix N II and find a configuration where the number of non-flippable plaquettes N − N V is minimal.The best strategy is to densely pack type-II vertices so that four type-II vertices block the same plaquette, just like in the square configurations.Detaching a type-II vertex would generate two new flippable plaquettes.Consequently, we get a lower bound for the number of nonflippable plaquettes: 2N II /4 ≤ N − N V .From this, Eq. (5b) follows.
N ≤ N V + N II : Any non-flippable plaquette has at least two type-II vertices; in a specific case, it has four [solid circles in Fig. 26(b)-(d)].The vertex type of the open circles is not determined; they can be type-I or II.We can construct loops of non-flippable plaquettes by connecting them via the type-II vertices denoted by the solid circles in the case of (b) and (c), or the loops intersect (d).These loops must close in a finite system and contain the same number of type-II vertices and non-flippable plaquettes.So if there are non-intersecting loops only, N − N V ≤ N II because there might be unaccounted type-II vertices on empty circles.It is equivalent to Eq. (5c).The intersection of the loops does not violate the inequality.
We recall that boundary conditions fundamentally influence the above inequalities.If we have open boundary conditions or infinite system size, it may not be true that the number of the non-flippable plaquettes is N − N V , even the number of the plaquettes and the vertices might be different.The condition for a second-order phase transition is e 2 = 0 and e 4 > 0. This happens along the line until the e 4 changes sign and becomes negative at µ = −3t/4 and V = −t (we note that e 6 is positive for V < −t/2).It signifies a tricritical point where the phase transition changes from a continuous (for −t < V ) to a first-order one (for V < −t).Unfortunately, our ED calculation is unsuitable for confirming the tricritical point's existence beyond the variational approach.Perturbation theory is a convenient tool to determine energy corrections in the quasi-classical regime caused by quantum fluctuations.The ground state energy is expanded in powers of t as where ε (0) is the classical energy and the second and fourth order corrections, following e.g.Ref. [64], read Above, we get |X⟩ by flipping a single plaquette in the classical ground state configuration |0⟩.The |Y ⟩ ̸ = |0⟩ are states constructed by flipping additional plaquettes in |X⟩.Since the energy is an even function of t [65], the correction terms proportional to an odd power in t must vanish.
Let us start with fully the flippable states and choose |0⟩ = |FF1⟩ in Eq. (10) and the 4th order is Let us turn to the square states.Choosing |0⟩ = |SqA⟩ in Eq. (10) and repeating the same steps for the square phase, we obtain the |X⟩ and |Y ⟩ intermediate states shown in Fig. 28, with energies where L x and L y are the horizontal and vertical lengths of the periodic cluster, and z x , z y ∈ Z (we set g Dis = 2 since the disordered manifold is translational invariant as a set of ice configurations).
All the structure factors we computed above reflect all the symmetries of D 4 point group, just like for the S(q) and S(q) calculated from the fully symmetric ground state of a finite cluster with periodic boundary conditions.Furthermore, all of the S(q) above satisfy the sum rule given by Eq. ( 60).

FIG. 1 .
FIG. 1. (a)The two-in, two-out configurations in the six-vertex model.While all of them are equivalent in three dimensions, they become distinguishable for the two-dimensional ice and split into type-I and type-II vertices.(b) In the fully packed loop representation, the thick bonds are occupied (i.e., part of a loop), whereas the thin bonds are empty.The figure shows the correspondence between the arrows and the occupied or vacant bonds for "even" sites, while the occupancies are reversed for "odd" sites.The "even/odd" refers to the parity of the sum ix + iy of a site's coordinates (ix, iy) in the lattice.(c) In the alternative Baxter representation, the bonds are colored blue if the arrows on the bond are reversed compared to the reference vertex in the first column.

FIG. 2 .
FIG. 2. (a)The off-diagonal term |⟲⟩⟨⟳| + |⟳⟩⟨⟲| reverses the direction of arrows along a directed loop around a square plaquette.(b) In the fully packed loop representation, the off-diagonal term acts on configurations where a pair of opposite bonds is empty and the other is occupied.(c) In the alternative Baxter representation, the four gray colors become blue and vice versa.

1 FIG. 4 .
FIG. 4. (a) The N = 16 square-shaped cluster with 16 sites and periodic boundary conditions is defined by the lattice vectors g1 = (4, 0) and g2 = (0, 4).The shown ice-rule obeying configuration, with 12 type-I and four type-II vertices, is in the (0, 0) flux sector.All the plaquettes are flippable except for the four denoted by crosses.(b) Reversing the arrows along the magenta loop crossing the boundaries, we get a configuration in the m = (2, 0) flux sector.(c) The N = 32 site cluster defined by g1 = (4, 4) and g2 = (−4, 4) from the N = 2L 2 family of cluster.The presented configuration is in the lowest non-zero m = (2, 2) flux sector.It originates from a periodic configuration of arrows where all the plaquettes are flippable (we call it a fully flippable configuration later on) by reversing the arrows along the magenta path.

FIG. 5 .
FIG. 5.The degeneracy of states (below) and maximal value of NV (above) in the different flux sectors (mx, my) of the 16-site cluster with periodic boundary conditions.

FIG. 6 .
FIG.6.The map shows the total number of flippable plaquettes (NV ) as well as the number of type-I and type-II vertices (NI and NII, where NI + NII = N ) in the classical basis of the six-vertex configurations.The convex hull is a triangle, and the two fully flippable (FF) configurations shown in Fig.8, the four square (Sq) configurations of Fig.9, and the isolated manifold, shown Fig.11(a) and (c), are located at the corners.

FIG. 10 .
FIG. 10.The elements of the symmetry group on the square lattice.(a) The symmetries that leave all three phases in the (0, 0) flux sector invariant.(b) The elements of the G/ G ∼ = D4 group, their action is summarized in TableI.The coordinate of the A plaquette's lower left corner (magenta point) is (0, 0).(c) The symmetry groups of the different phases and the subgroup relations.

0FIG. 11 .
FIG. 11.The isolated configuration in the m = (L, L) flux sector in (a) the arrow and (b) the fully packed loop representation.(c) An isolated configuration from the m = (0, L) flux sector.Flippable plaquettes are absent, and all the vertices are type-II (NV = 0 and NII = N ).The order of the left-and right-pointing horizontal lines is arbitrary, so the number of states in this flux sector is 4 2 = 6.(d) A mobile leapfrog excitation ('frogton') with two flippable plaquettes (NV = 2, NI = 2) in a cluster with suitable cluster geometry.(e) A state in the m = (L − 2, L − 2) flux sectors.It results from reversing the arrows in the isolated state shown in (a) along the closed magenta loop crossing the boundary.(f) The minimal extent of the isolated state manifold.The µ = V − t boundary (solid line) denotes the instability line against frogtons.The isolated states, however, can extend beyond the dashed boundary at V = t.

FIG. 12 .
FIG. 12.The finite size scaling of the gap ∆ = E (L−2,L−2) GS − EIso between the ground state energy in the m = (L − 2, L − 2) flux sector and the EIso = N µ of the isolated state on the V = µ + t line in the phase diagram.We show three selected cases (µ/t = 1, 2, and 4) in N = L 2 class of clusters, calculated by the Lánczos method.The size of the Hilbert space is 90 for L = 4 [see Fig. 5], 2 772 for L = 6 [see Fig. 25(b)], 51 480 for L = 8, and 923 780 for L = 10 in the (L − 2, L − 2) flux sectors.The gap exponentially vanishes with the linear size of the cluster L (note the logarithmic vertical axis), and the solid lines show the fit to the ae −bL function.The b values depend weakly on the value of µ.

FIG. 13 .
FIG.13.The ground state flux sectors in the 32 site cluster.Besides the m = (0, 0) (magenta) and isolated sectors (white), other sectors emerge in a tiny region fanning out from the RK point (the gray area).We will discuss them in Sec.VIII.

FIG. 14 .
FIG. 14. Exact diagonalization results for the 32-site cluster at V = 0 as a function of µ/t.(a) The density of the type-I (nI = NI/N ) and II (nII = NII/N ) vertices in the ground state.(b)The expectation values of the squared order parameters and (c) their numerical derivatives with respect to µ/t.We identify inflection points of the ⟨ Õ2 ⟩ (i.e., the extrema of the ⟨ Õ2 ⟩ ′ ) with phase boundaries.(d) The low-energy excitations are defined by their momentum and charge conjugation parity.The size of the symbols indicates the overlap between the initial states given by Eqs.(31) and the energy eigenstates.The ground state is fully symmetric.The level crossing of the lowest lying (π, π)e and (π, π)o symmetry levels at µ/t ≈ 0.29 serves as an alternative indicator of the phase boundary between the plaquette (−0.78 ≲ µ/t ≲ 0.29) and the fully flippable phases (µ/t ≳ 0.29).Avoided level crossings around µ/t ≈ −0.78 in the first excitations of the (0, 0)e and (π, π)e symmetry sectors indicate the boundary between the square (µ/t ≲ −0.78) and the plaquette phase.

FIG. 16 .
FIG.16.The ground state phase diagram based on ED for the 32site cluster.The phase boundaries agree with the N = 36 result up to two decimal places.We present the diagram for the smaller systems because of the better resolution.The points denote the positions of the (avoided) level crossings and inflection points of the squared order parameters.The boundary between the phase with finite flux sectors and the isolated states (solid line) is the exact one of Eq. (29).The red point denotes the Rokshar-Kivelson (RK) point.

FIG. 17 .
FIG. 17.Comparison of the phase boundaries calculated by ED for a 32-site cluster and by the fourth-order perturbation expansion in t.We show the Padé approximants of the perturbation series to estimate its convergence.(a) ED data and Padé approximants of the boundary between the square and the isolated phases near the triple point.The coordinates for the triple point are cluster dependent, for N = 32 V /t = 0.32 and µ/t = −0.83,while for N = 36 are V /t = 0.29 and µ/t = −0.87.(b) The perturbational curve gives the boundary between the fully flippable and the square phases.It meets the lines (FF-Pl and Sq-Pl) from the ED at the triple point V /t ≈ −1.75 and µ/t ≈ −1.01.The different Padé approximants do not deviate significantly in the relevant V /t ≲ −1.75 range.The grey dashed line shows the V = 2µ classical phase boundary.

2 FIG. 18 .
FIG.18.Finite-size scaling of the densities of (a) type-II vertices nII and (b) flippable plaquettes nV in the RK-wave function for zero flux sectors.We show ED data for smaller size clusters (N = 36 and 32) and Monte Carlo data for larger (up to N = 576 and 512).The solid lines show the fit to the finite size corrections, Eq. (41).

FIG. 19 .+ m 2 y
FIG. 19.The expectation value of the densities of (a) type-II vertices nII and (b) flippable plaquettes nV in the RK-wave function grows (decreases) linearly with the square of the total flux, m 2 = |m| 2 = m 2 x + m 2 y .The evaluation is exact numerically for N = 32 and 36, and we sampled the RK wave function by Monte Carlo for larger system sizes of up to 576 sites.For consistency, we divide m 2 x + m 2 y by N for the N = L 2 clusters (N = 36, 64, 256, 576, denoted by red colors in the plots) and by 2N for the N = 2L 2 size clusters (N = 32, 128, 512, green in the plots).The black straight lines show the extrapolation to the thermodynamic limit.

4 FIG. 22 .
FIG. 22. (a)The phase diagram around the V = t and µ = 0 RK point.The parameter θ is defined by Eq.(50).A first-order phase transition occurs between the isolated states and the plaquette phase at θ = 1.264π.For π/4 < θ < 0.291π, when µ is positive, the two phases are separated by a region in which the flux sectors interpolate from the isolated manifold with maximal flux to the m = (0, 0).The boundary θ = π/4 (line V = µ + t for µ > 0) of the isolated states is exact and also holds away from the RK-point.The m = (0, 0) flux sector is the ground state for 0.291π < θ < 1.264π.(b) We plot m 2 = m 2x + m 2 y of the flux sectors with minimal energy as a function of the parameter θ in the vicinity of the RK-point [δ → 0 in Eq. (50)] for π/4 < θ < 0.291π, where we expect the devil's staircase.The flux monotonically decreases with increasing ϑ.The plot suggests a finite-width plateau at m 2 /2N = 1/4, corresponding to the flux sector m = (L/2, L/2).Calculations were done for up to 576 sites in the N = L 2 type clusters.(c) The finite size scaling of the width of the m = (L/2, L/2) plateau indicates a tiny but finite width ∆θ 1/2 = 0.002π in the thermodynamic limit.

C
. Comparison with the structure factors in the artificial quantum 6-vertex model

FIG. 26 .
FIG. 26.The four symmetrically inequivalent vertex configurations at the corners of a plaquette.(a) If the plaquette is flippable, the vertices are undetermined (open circles).(b) and (c) A non-flippable plaquette has at least two or (d) four type-II vertices (closed circles).The remaining vertices are undetermined (open circles).The red lines are the Faraday loops passing through the non-flippable vertices.

FIG. 27 .
FIG.27.The excited states |X⟩ and |Y ⟩ appearing in the fourthorder perturbation expansion Eq. (D3) of the ground state energy in the fully flippable phase.Flipped plaquettes are in magenta and gray crosses indicate non-flippable plaquetes.The numbers on the long arrows show the number of transitions, where N is the system size.

TABLE II .
The value of the order parameters in different phases.Let us note that the OPl is finite in both the plaquette and the square phase.