Crystalline phases and devil's staircase in qubit spin ice

Motivated by the recent realization of an artificial quantum spin ice in an array of superconducting qubits with tunable parameters [King et al., Science 373, 576 (2021)], we scrutinize a quantum six-vertex model on the square lattice that distinguishes type-I and type-II vertices. We map the zero-temperature phase diagram using numerical (exact diagonalization) and analytical (perturbation expansion, Gershgorin theorem) methods. Following a symmetry classification, we identify three crystalline phases alongside a subextensive manifold of isolated configurations. Monte Carlo simulations at the multicritical Rokhsar-Kivelson point provide evidence for a quantum phase exhibiting a cascade of transitions with increasing flux. By comparing structure factors, we find evidence for the emergence of the fully flippable and plaquette phases in the artificial quantum spin ice.

Figure 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 $i_x+i_y$ of a site's coordinates $(i_x,i_y)$ 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.

Figure 2

(a) The off-diagonal term $\mid ⥀\rangle \langle ⥁\mid +\mid ⥁\rangle \langle ⥀\mid$ 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.

Figure 3

The phase diagram of the quantum six-vertex model, defined by Hamiltonian (3), consists of three phases. For $V\lesssim -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\lesssim V⩽t$. The approximate wave function is a direct product of the $|◯\rangle =\frac{1}{\sqrt{2}}(\mid ⥁\rangle +\mid ⥀\rangle )$ on the resonating plaquettes [ in the loop representation]. The Rokshar- Kivelson point at $V=t$ is quantum critical, and for $0\le t<V$, the isolated manifold appears with configurations having no flippable plaquettes (denoted by crosses).

Figure 4

(a) The $N=16$ square-shaped cluster with 16 sites and periodic boundary conditions is defined by the lattice vectors $g_1=(4,0)$ and $g_2=(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 $\mathbf{m}=(2,0)$ flux sector. (c) The $N=32$-site cluster defined by $g_1=(4,4)$ and $g_2=(-4,4)$ from the $N=2L^2$ cluster family. The presented configuration is in the lowest nonzero $\mathbf{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.

Figure 5

The degeneracy of states (below the solid circles) and maximal value of $N_V$ (above) in the different flux sectors $(m_x,m_y)$ of the 16-site cluster with periodic boundary conditions.

Figure 6

The map shows the total number of flippable plaquettes ($N_V$) as well as the number of type-I and type-II vertices ($N_{\text{I}}$ and $N_{\text{II}}$, where $N_{\text{I}}+N_{\text{II}}=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 below, the four square configurations of Fig. 9 below, and the isolated manifold, shown Figs. 11 and 11 below, are located at the corners.

Figure 7

The classical phase diagram in the parameter space of the $V$ and $\mu$, the chemical potential of the type-II vertices. The isolated and square phases consist only of type-II vertices; they become favorable when $\mu$ 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 below) and four square states (see Fig. 9 below). 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.

Figure 8

The two fully flippable configurations in (a) arrow and (b) fully packed loop representation. The open circles serve as anchor points.

Figure 9

The four square configurations made from type-II vertices in (a) the arrow and (b) the fully packed loop representation. Compared to the fully flippable configurations in Fig. 8, only half of the plaquettes are flippable (grey crosses denote the nonflippable plaquettes). These four states break the translational invariance and partition the lattice into four sublattices denoted by letters A to D. We assign the letters A, B, C, and D to the plaquette in which arrows rotate counterclockwise. The site with the open circle is the same in all clusters and serves as anchor points.

Figure 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 $\mathsf{G}/\tilde{\mathsf{G}}\cong {\mathsf{D}}_{\mathsf{4}}$ group, their action is summarized in Table 1. The coordinate of the A plaquette's lower left corner (magenta square) is (0,0). (c) The symmetry groups of the different phases and the subgroup relations.

Figure 11

The isolated configuration in the $\mathbf{m}=(L,L)$ flux sector in (a) the arrow and (b) the fully packed loop representation. (c) An isolated configuration from the $\mathbf{m}=(0,L)$ flux sector. Flippable plaquettes are absent, and all the vertices are type-II ($N_V=0$ and $N_{\text{II}}=N$). The order of the left- and right-pointing horizontal lines is arbitrary, so the number of states in this flux sector is $\left(\genfrac{}{}{0pt}{}{4}{2}\right)=6$. (d) A mobile leapfrog excitation (“frogton”) with two flippable plaquettes ($N_V=2, N_{\text{I}}=2$) in a cluster with suitable cluster geometry. (e) A configuration in the $\mathbf{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 manifold. The $\mu =V-t$ boundary (solid line) denotes the instability line against frogtons. The isolated states, however, can extend beyond the dashed boundary at $V=t$.

Figure 12

The finite-size scaling of the gap $\mathrm{\Delta }=E_{\text{GS}}^{(L-2,L-2)}-E_{\text{Iso}}$ between the ground-state energy in the $\mathbf{m}=(L-2,L-2)$ flux sector and the $E_{\text{Iso}}=N\mu$ of the isolated state on the $V=\mu +t$ line in the phase diagram. We show three selected cases ($\mu /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 below], 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 fits to the $a{e}^{-bL}$ function. The $b$ values depend weakly on the value of $\mu$.

Figure 13

The ground-state flux sectors in the 32-site cluster. Besides the $\mathbf{m}=(0,0)$ (magenta) and isolated sectors (light gray area), other sectors emerge in a tiny region fanning out from the RK point (the gray area). We will discuss them in Sec. 8.

Figure 14

Exact diagonalization results for the 32-site cluster at $V=0$ as a function of $\mu /t$. (a) The density of the type-I ($n_{\text{I}}=N_{\text{I}}/N$) and II ($n_{\text{II}}=N_{\text{II}}/N$) vertices in the ground state. (b) The expectation values of the squared order parameters and (c) their numerical derivatives with respect to $\mu /t$. We identify inflection points of the $\langle {\tilde{O}}^2\rangle$ (i.e., the extrema of the ${\langle {\tilde{O}}^2\rangle }^{\prime }$) with phase boundaries. (d) The low-energy excitations are defined by their momenta and conjugation parities. The size of the symbols indicates the overlap between the initial states given by Eqs. (30) and (31) and the energy eigenstates. The ground state is fully symmetric. The level crossing of the lowest-lying ${(\pi ,\pi )}_e$ and ${(\pi ,\pi )}_o$ symmetry levels at $\mu /t\approx 0.29$ serves as an alternative indicator of the phase boundary between the plaquette ($-0.78\lesssim \mu /t\lesssim 0.29$) and the fully flippable phases ($\mu /t\gtrsim 0.29$). Avoided level crossings around $\mu /t\approx -0.78$ in the first excitations of the ${(0,0)}_e$ and ${(\pi ,\pi )}_e$ symmetry sectors indicate the boundary between the square ($\mu /t\lesssim -0.78$) and the plaquette phase.

Figure 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 $\mathcal{C}$ parity.

Figure 16

The ground-state phase diagram based on ED for the 32-site 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.

Figure 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=32V/t=0.32$ and $\mu /t=-0.83$, while for $N=36$ are $V/t=0.29$ and $\mu /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\approx -1.75$ and $\mu /t\approx -1.01$. The different Padé approximants do not deviate significantly in the relevant $V/t\lesssim -1.75$ range. The grey dashed line shows the $V=2\mu$ classical phase boundary.

Figure 18

Finite-size scaling of the densities of (a) type-II vertices $n_{\text{II}}$ and (b) flippable plaquettes $n_V$ in the RK wave function for (0,0) flux sector. We present 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).

Figure 19

The expectation value of the densities of (a) type-II vertices $n_{\mathrm{II}}$ and (b) flippable plaquettes $n_V$ in the RK-wave function grows (decreases) linearly with the square of the total flux, $m^2={\left|\mathbf{m}\right|}^2=m_x^2+m_y^2$. 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_x^2+m_y^2$ 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.

Figure 20

The finite-size dependence of the gaps ${\mathrm{\Delta }}_V\left(\mathbf{m}\right)={\langle N_V\rangle }_{\mathbf{m}}-{\langle N_V\rangle }_{(0,0)}$ and ${\mathrm{\Delta }}_{\text{II}}\left(\mathbf{m}\right)={\langle N_{\text{II}}\rangle }_{\mathbf{m}}-{\langle N_{\text{II}}\rangle }_{(0,0)}$ in a few selected flux sectors close to $\mathbf{m}=(0,0)$. The gaps collapse to a common value in the thermodynamic limit when divided by $m^2/m_0^2$. $m_0^2$ is the square of the unit flux, for the $L^2$ cluster $m_0^2=4$ and for the $2L^2$ cluster $m_0^2=8$.

Figure 21

The dependence of the gap of flippable plaquettes and type-II vertices between the smallest nonzero- and the zero-flux sector, $\langle \mathrm{\Delta }{N}_V\rangle ={\langle \mathrm{\Delta }{N}_V\rangle }_{(0,2)}-{\langle \mathrm{\Delta }{N}_V\rangle }_{(0,0)}$ and $\langle \mathrm{\Delta }{N}_{\text{II}}\rangle ={\langle \mathrm{\Delta }{N}_{\text{II}}\rangle }_{(0,2)}-{\langle \mathrm{\Delta }{N}_{\text{II}}\rangle }_{(0,0)}$ on the aspect ratio of the cluster. We calculated the gap for $N=256$ site rectangular clusters of different shapes: $L_x\times L_y=4\times 64,8\times 32,16\times 16,32\times 8$, and $64\times 4$. We run 25 Monte Carlo simulations for each cluster with $5\times {10}^8$ steps each. We show $-\langle \mathrm{\Delta }{N}_V\rangle$ for visual convenience. The large variation of the gap on the shape of the clusters is almost perfectly accounted for by the aspect ratio $L_y/L_x$ while keeping $L_x{L}_y$ constant: $\langle \mathrm{\Delta }{N}_V\rangle \approx -1.52L_y/L_x$ and $\langle \mathrm{\Delta }{N}_{\text{II}}\rangle \approx 1.16L_y/L_x$.

Figure 22

(a) The phase diagram around the $V=t$ and $\mu =0$ RK point. The parameter $\theta$ is defined by Eq. (50). A first-order phase transition occurs between the isolated states and the plaquette phase at $\theta =1.264\pi$. For $\pi /4<\theta <0.291\pi$, when $\mu$ 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 $\mathbf{m}=(0,0)$. The boundary $\theta =\pi /4$ (line $V=\mu +t$ for $\mu >0$) of the isolated states is exact and also holds away from the RK point. The $\mathbf{m}=(0,0)$ flux sector is the ground state for $0.291\pi <\theta <1.264\pi$. (b) We plot $m^2=m_x^2+m_y^2$ of the flux sectors with minimal energy as a function of the parameter $\theta$ in the vicinity of the RK point [$\delta \to 0$ in Eq. (50)] for $\pi /4<\theta <0.291\pi$, where we expect the devil's staircase. The flux monotonically decreases with increasing $\vartheta$. The plot suggests a finite-width plateau at $m^2/2N=1/4$, corresponding to the flux sector $\mathbf{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 $\mathbf{m}=(L/2,L/2)$ plateau indicates a tiny but finite width $\mathrm{\Delta }{\theta }_{1/2}=0.002\pi$ in the thermodynamic limit.

Figure 23

(a) The bond correlation function $C^h(x,y)$ [Eq. (56)] in real space for the 36-site cluster (the horizontal bond at the center is at (1/2,0) in lattice coordinates). It is equal to $C^{hh}(x,y)$ on the horizontal bonds, defined by Eq. (57a), and to $C^{vh}(x+\frac{1}{2},y-\frac{1}{2})$ on the vertical bonds [see Eq. (57b)]. From left to right, we present correlations in the square, plaquette, fully flippable phases, at the RK point, and in the classical disordered manifold. The area of the disks is proportional to the value of $S(x,y)$; blue indicates positive and red negative values. (b) The density plot of the structure factor $S\left(\mathbf{q}\right)$. The green dashed square encloses the Brillouin zone of the 36-site square cluster, while the green dotted line is the boundary of the extended Brillouin zone containing $2N=72\mathbf{q}$ points. (c) The $S\left(\mathbf{q}\right)$ along the path $\mathbf{q}=(\pi ,\pi )\to (0,0)\to (2\pi ,0)\to (\pi ,\pi )\to (\pi ,0)$ drawn in white in the leftmost panel in (b). The structure factor diverges with system size at $\mathbf{Q}=(\pi ,0)$ in the square phase (first column) and at $\mathbf{Q}=(2\pi ,0)$ in the fully flippable phase (third column). The structure factor is diffuse in the plaquette phase, with a peak centered at $\mathbf{Q}=(2\pi ,0)$. At the RK point (fourth column), the value of the structure factor strongly depends on the direction we approach the $(\pi ,\pi ), S(\pi +\delta ,\pi +\delta ):S(\pi ,\pi +\delta ):S(\pi +\delta ,\pi -\delta )=0:1:2$ as $\delta \to 0$, demonstrating the nonanalytic behavior of the pinch point. In the disordered manifold (fifth column), we see subdivergent lines along $\mathbf{Q}=(\pm 2\pi ,q)$ and $\mathbf{Q}=(q,\pm 2\pi )$.

Figure 24

The magnetic structure factor $\tilde{S}\left(\mathbf{q}\right)$ 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 $\mu /t=-1$, (b) the plaquette phase for $\mu =0$, and (c) the fully flippable for $\mu /t=0.5$, keeping $V=0$ in all three cases. We find divergent Bragg peaks in the square and the fully flippable phases, while $\tilde{S}\left(\mathbf{q}\right)$ remains diffuse in the plaquette phase. (d) $\tilde{S}\left(\mathbf{q}\right)$ at the RK point ($V=t$ and $\mu =0$) displays the pinch points at $\mathbf{Q}=(2\pi z_x,2\pi z_y)$, where $z_x,z_y\in \mathbb{Z}$. (e) The magnetic structure factor of the $V=t=0$ classical disordered manifold. The amplitudes of the subdivergent horizontal and vertical lines are given out in Appendix pp6-s4. The green dotted square denotes the boundary of the extended Brillouin zone, the same as in Fig. 23.

Figure 25

The flux sectors, their degeneracy (below the circles), and max value of $N_V$ (above the circles) for the (a) 32-site and (b) $N=36$-site clusters with periodic boundary conditions. For visualization purposes, only one irreducible octant is shown. The full diagram follows from symmetry.

Figure 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),(c)] A nonflippable 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 nonflippable vertices.

Figure 27

The excited states $|X\rangle$ and $|Y\rangle$ appearing in the fourth-order perturbation expansion Eq. (D3) of the ground-state energy in the fully flippable phase. Flipped plaquettes are in magenta and gray crosses indicate nonflippable plaquettes. The numbers on the long arrows show the number of transitions, where $N$ is the system size.

Figure 28

Same as Fig. 27, but for the square phase.

Figure 29

Comparison of the structure factor $S\left(\mathbf{q}\right)$ in the plaquette phase calculated by ED (open circles for N=16, 32, and 36 sites) with the $S\left(\mathbf{q}\right)$ of the factorized plaquette wave function, Eq. (F8) (solid line). We present two cases, (a) one close to the square phase with $\mu /t=-0.5$, and (b) one close to the flippable phase, where $\mu /t=0$ ($V=0$ for both).