Constructing Quantum Spin Liquids Using Combinatorial Gauge Symmetry

TL;DR

This paper introduces combinatorial gauge symmetry, a new local transformation involving spin rotations and permutations, enabling the construction of quantum spin liquids with accessible two-body interactions and potential applications in topological qubits.

Contribution

It defines combinatorial gauge symmetry and links it to Hadamard matrices, providing a novel framework for designing quantum spin liquids with simple interactions.

Findings

Identifies combinatorial gauge symmetry in Hamiltonians with Hadamard coupling matrices.

Shows how this symmetry can be used to construct quantum spin liquids.

Connects the symmetry to automorphisms of Hadamard matrices and topological qubit design.

Abstract

We introduce the notion of combinatorial gauge symmetry -- a local transformation that includes single spin rotations plus permutations of spins (or swaps of their quantum states) -- that preserve the commutation and anti-commutation relations among the spins. We show that Hamiltonians with simple two-body interactions contain this symmetry if the coupling matrix is a Hadamard matrix, with the combinatorial gauge symmetry being associated to the automorphism of these matrices with respect to monomial transformations. Armed with this symmetry, we address the physical problem of how to build quantum spin liquids with physically accessible interactions. In addition to its intrinsic physical significance, the problem is also tied to that of how to build topological qubits.

Tables

Table 1. Table 1: Eight ground states in the classical limit J ≫ Γ , Γ ~ much-greater-than 𝐽 Γ ~ Γ J\gg\Gamma,\widetilde{\Gamma}

${\sigma }_1^z$	${\sigma }_2^z$	${\sigma }_3^z$	${\sigma }_4^z$	${\mu }_1^z$	${\mu }_2^z$	${\mu }_3^z$	${\mu }_4^z$
$\uparrow$	$\uparrow$	$\uparrow$	$\uparrow$	$\downarrow$	$\downarrow$	$\downarrow$	$\downarrow$
$\uparrow$	$\uparrow$	$\downarrow$	$\downarrow$	$\uparrow$	$\uparrow$	$\downarrow$	$\downarrow$
$\uparrow$	$\downarrow$	$\uparrow$	$\downarrow$	$\uparrow$	$\downarrow$	$\uparrow$	$\downarrow$
$\uparrow$	$\downarrow$	$\downarrow$	$\uparrow$	$\uparrow$	$\downarrow$	$\downarrow$	$\uparrow$
$\downarrow$	$\downarrow$	$\uparrow$	$\uparrow$	$\downarrow$	$\downarrow$	$\uparrow$	$\uparrow$
$\downarrow$	$\uparrow$	$\downarrow$	$\uparrow$	$\downarrow$	$\uparrow$	$\downarrow$	$\uparrow$
$\downarrow$	$\uparrow$	$\uparrow$	$\downarrow$	$\downarrow$	$\uparrow$	$\uparrow$	$\downarrow$
$\downarrow$	$\downarrow$	$\downarrow$	$\downarrow$	$\uparrow$	$\uparrow$	$\uparrow$	$\uparrow$

Equations

$$\left(\begin{array}{c}{\vec{\sigma }}_1\\ {\vec{\sigma }}_2\\ {\vec{\sigma }}_3\\ {\vec{\sigma }}_4\end{array}\right)\to \left(\begin{array}{cccc}0& 0& g_1& 0\\ 0& 0& 0& g_2\\ g_3& 0& 0& 0\\ 0& g_4& 0& 0\end{array}\right)\left(\begin{array}{c}{\vec{\sigma }}_1\\ {\vec{\sigma }}_2\\ {\vec{\sigma }}_3\\ {\vec{\sigma }}_4\end{array}\right).$$

$$H=-\sum _s\left[J\sum _{\begin{array}{c}a\in s\\ i\in s\end{array}}{W}_{ai}{\sigma }_i^{\mathrm{z}}{\mu }_a^{\mathrm{z}}+\Gamma \sum _{a\in s}{\mu }_a^{\mathrm{x}}\right]-\tilde{\Gamma }\sum _i{\sigma }_i^{\mathrm{x}},$$

$${\sigma }_i^{\mathrm{z}}$$

$${\mu }_a^{\mathrm{z}}$$

$$W=\left(\begin{array}{cccc}-1& +1& +1& +1\\ +1& -1& +1& +1\\ +1& +1& -1& +1\\ +1& +1& +1& -1\end{array}\right).$$

$$L=$$

$$G_p=\prod _{s\in p}{\mathcal{L}}_s^{(\mu )}\prod _{s\in p}{\mathcal{R}}_s^{(\sigma )},$$

$$E_a^{(\pm )}({\sigma }_1^{\mathrm{z}},{\sigma }_2^{\mathrm{z}},{\sigma }_3^{\mathrm{z}},{\sigma }_4^{\mathrm{z}})=\pm {\left[J^2{\left(\sum _{i=1}^4{W}_{ai}{\sigma }_i^{\mathrm{z}}\right)}^2+{\Gamma }^2\right]}^{1/2}.$$

$$E_a^{(\pm )}=\pm C_0\pm C_2\sum _{i\ne j}^4{W}_{ai}{W}_{aj}{\sigma }_i^{\mathrm{z}}{\sigma }_j^{\mathrm{z}}$$

$$\pm C_4{W}_{a1}{W}_{a2}{W}_{a3}{W}_{a4}{\sigma }_1^{\mathrm{z}}{\sigma }_2^{\mathrm{z}}{\sigma }_3^{\mathrm{z}}{\sigma }_4^{\mathrm{z}},$$

$$H_{\mathrm{eff}}^{\mathrm{star}}=\gamma -\lambda {\sigma }_1^{\mathrm{z}}{\sigma }_2^{\mathrm{z}}{\sigma }_3^{\mathrm{z}}{\sigma }_4^{\mathrm{z}},$$

$$H_{\mathrm{eff}}=-\lambda \sum _s\prod _{i\in s}{\sigma }_i^{\mathrm{z}}-\tilde{\Gamma }\sum _i{\sigma }_i^{\mathrm{x}}.$$

$$H=$$

$$-\sum _p\left[J\sum _{\begin{array}{c}b\in p\\ j\in p\end{array}}{W}_{bj}{\sigma }_j^{\mathrm{x}}{\tau }_b^{\mathrm{x}}+\Gamma \sum _{b\in p}{\tau }_b^{\mathrm{z}}\right].$$

$$F_s=\prod _{p\in s}{\mathcal{L}}_p^{(\tau )}\prod _{i\in s}{\sigma }_i^{\mathrm{z}}.$$

$$H=-\sum _s\left[J\sum _{\begin{array}{c}a\in s\\ i\in s\end{array}}{W}_{ai}{\sigma }_i^{\mathrm{z}}{\mu }_a^{\mathrm{z}}\right]+g\sum _p{G}_p.$$

$$H^{\mathrm{star}}=-J\sum _{\begin{array}{c}a\in s\\ i\in s\end{array}}{W}_{ai}{\sigma }_i^{\mathrm{z}}{\mu }_a^{\mathrm{z}}-\Gamma \sum _{a\in s}{\mu }_a^{\mathrm{x}}$$

$$\gamma$$

$$\lambda$$

$$w_{{\mathbf{\sigma }}_{\ell }}:=\mathrm{tr}\left({\rho }_{\mathbf{\sigma }}|{\mathbf{\sigma }}_{\ell }⟩⟨{\mathbf{\sigma }}_{\ell }|\right).$$

$$H_{\mathrm{eff}}^{3D-\mathrm{toric}}=-\lambda \sum _p\prod _{i\in p}{\sigma }_i^{\mathrm{z}}-\tilde{\Gamma }\sum _i{\sigma }_i^{\mathrm{x}}.$$

$$H_{\mathrm{eff}}^{\mathrm{star}}$$

$$=3\gamma -\lambda \left(B_s^{xy}+B_s^{yz}+B_s^{xz}\right),$$

$$H_{\mathrm{eff}}^{\mathrm{X}-\mathrm{cube}}=-\lambda \sum _s\left(B_s^{xy}+B_s^{yz}+B_s^{xz}\right)-\tilde{\Gamma }\sum _i{\sigma }_i^{\mathrm{x}}.$$

Full text

Quantum liquids of spins are systems where no magnetic symmetry-breaking order should be detectable down to zero temperature Balents (2010), and instead topological order exits Wen (1990). On the theoretical side, there are a number of model Hamiltonians where quantum spin liquid states exist Kitaev (2003, 2006). Gauge symmetries are common in these models, whether discrete or continuous, intrinsic or emergent. Many of these gauge models, such as the ${\mathbb{Z}}_{2}$ toric code Kitaev (2003) and fracton models such as the $X$-cube Castelnovo et al. (2010); Vijay et al. (2016), are defined using multi-spin interactions. Here, we show that exact local ${\mathbb{Z}}_{2}$ gauge symmetries in these models can arise from solely two-spin interactions. That one can generate effective multi-spin interactions in some low energy limit of a two-spin Hamiltonian is not unexpected; what is novel is that the symmetries we discuss are exact. We articulate a notion of combinatorial gauge symmetry that underlies why it is possible to construct local two-spin Hamiltonians with an exact ${\mathbb{Z}}_{2}$ gauge symmetry.

Algebra-preserving transformations and monomial matrices – We start with a set of $N$ spin-1/2 degrees of freedom, such as the familiar spin models on a lattice with $N$ sites. The spin operators are Pauli matrices $\sigma^{\alpha}_{i}$, where $\alpha={\rm x,y,z}$ and $i=1,\dots,N$. Spins on different sites commute, while those on the same site satisfy the usual angular momentum algebra. Let us ask a simple question: which transformations of these $3N$ operators can preserve all commutation and anti-commutation relations? For $N$ bosons or fermions, this is a trivial question to answer; the allowed set of single-particle transformations belong to the unitary group $U(N)$ because either the commutation or anti-commutation relations need to be satisfied. But for spins, the question is harder; one cannot simply mix spatial components of different spins and retain both the intra- and inter-site algebra.

The Hilbert space for $N$ spins is $2^{N}$-dimensional and the allowed operators in this space are $2^{N}\times 2^{N}$ unitary matrices, corresponding to the group $SU(2^{N})$. A generic transformation on the spin operators, $\sigma_{i}^{a}\to U\,\sigma_{i}^{a}\,U^{\dagger}$ preserves the algebra, but also acts simultaneously on many spins: it mixes the $3N$ single-spin operators $\sigma^{a}_{i}$ with the other (multi-spin) $2^{2N}-1-3N$ generators of $SU(2^{N})$. Therefore, if one is to remain with only single-spin terms, one must work with a much smaller subgroup of $SU(2^{N})$. The simplest solution is trivial: only rotate spins individually by restricting the allowed transformations to $SU(2)\otimes SU(2)\otimes\cdots\otimes SU(2)$ or $N$ copies of $SU(2)$. A more interesting and non-trivial solution is to also allow permutations of spins. (If one wishes to connect to quantum gates, these transformations correspond to the combination of one-qubit rotations and the use of two-qubit SWAP gates.)

Any $SU(2)$ transformation on spin $i$ can be represented by a matrix in the rotation group $g_{i}\in SO(3)$ that acts on the spatial components of the vector $\vec{\sigma}_{i}=(\sigma^{\rm x}_{i},\sigma^{\rm y}_{i},\sigma^{\rm z}_{i})^{\top}$. This representation makes it convenient to combine permutations and single-spin transformations into monomial matrices. Monomial matrices are generalizations of permutation matrices such that the non-zero elements in each row and column are group elements, not simply equal to 1. Here is an $N=4$ example:

[TABLE]

It is clear from this form that monomial matrices are orthogonal and that the product of any two monomial matrices is another monomial matrix. The example above can be written as a product of the diagonal matrix Diag($g_{1},g_{2},g_{3},g_{4}$) and a $4\times 4$ permutation matrix.

For arbitrary $N$, the group of monomial matrices is a semidirect product of the group generated by the diagonal matrices Diag($g_{1},\cdots,g_{N}$) and the group of permutations (symmetric group) $S_{N}$. In mathematical literature, this particular form of a semidirect product is sometimes referred to as a wreath product.

To summarize the above: we are pointing out that many-body spin states admit a group of non-trivial transformations on the $3N$ spin components that preserve all spin algebras. When formulated in this way, the combination of local and permutation symmetry will allow us to construct exact lattice gauge theories using only two-body interactions.

Combinatorial gauge symmetry – One particular subgroup of monomial transformations, such as in Eq. (1), is for $SO(3)$ rotations by angle $\pi$ around a given axis, which we take to be $\hat{x}$. This is equivalent to flipping the $z$-component of spin. We shall use this special case to construct a microscopic model with local ${\mathbb{Z}}_{2}$ symmetry. We term our methodology combinatorial gauge symmetry for its relation to monomials and permutations.

Consider the lattice depicted in Fig. 1, where 4 “matter” spins $\mu$ are placed on each lattice site, and “gauge” spins $\sigma$ are placed on the links. A single site (star) is isolated in Fig. 1(a), and contains the 4 matter spins and 4 gauge spins sitting on the links. The gauge spins are shared by neighboring stars, as depicted in Fig. 1(b). Each matter spin couples only to its neighboring gauge spins but not to one another (or other lattice sites). Gauge spins do not couple to each other. We encode all two-spin ($ZZ$) couplings between $\mu^{\rm z}_{a}$ and $\sigma^{\rm z}_{i}$ by a $4\times 4$ matrix $W_{ai}$.

The quantum fluctuations will come from two transverse fields $\widetilde{\Gamma}$ and $\Gamma$ acting on the gauge spins and matter spins, respectively. For generality, we allow $\Gamma$ and $\widetilde{\Gamma}$ to have different magnitudes.

Thus the full lattice Hamiltonian is given by

[TABLE]

where the $s$ are stars on the lattice.

We shall select the interaction matrix $W$ so as to satisfy the monomial transformations as in Eq. (1) that act on the $z$-components of the gauge and matter spins as follows:

[TABLE]

These are monomial transformations that preserve the spin commutation and anticommutation relations, as discussed above. The $L$ (“left”) and $R$ (“right”) matrices act like gauge transformations on the $z$-components of the gauge and matter spins. These monomial matrices have elements $\pm 1$. (Henceforth all monomial matrices will be of this kind.)

The requirement that the Hamiltonian Eq. (2) be invariant with respect to transformations Eq. (3) is equivalent to the requirement that the $W$ matrices be invariant under the automorphism transformation $L^{-1}WR=W$, where $L$ and $R$ are $4\times 4$ monomial matrices Kantor (1969). [The transverse fields are also invariant under the transformation Eq. (3), and we shall return to this point below.]

Hadamard matrices Kantor (1969) satisfy these conditions. These matrices have elements $\pm 1$, and all its columns (or rows) are orthogonal vectors, i.e., $W^{\top}\,W\propto\openone$. (They maximize the determinant of the information matrix $W^{\top}\,W$.) We pick an intuitive form of $W$, where the coupling between $\sigma^{\rm z}_{i}$ and $\mu^{\rm z}_{a}$ is anti-ferromagnetic when $i=a$ and ferromagnetic otherwise:

[TABLE]

All other choices of $W$ are equivalent by symmetry and will not affect the spectrum. Specifically, any two Hadamard matrices $W$ and $W^{\prime}$ are equivalent if there exist monomial matrices $S_{1},S_{2}$ such that $W^{\prime}=S_{1}^{-1}WS_{2}$.

Our model further restricts $R$ to be diagonal because any off-diagonal permutation of gauge spins would deform the lattice. For example, with our choice of $W$ in Eq. (4), the following pair satisfies the conditions above:

[TABLE]

Once we choose an $R$, we determine $L$ uniquely by solving the automorphism condition: $L=W\,R\,W^{-1}$. (The number of $-1$’s in the diagonal $R$ matrix must be even so that the corresponding $L$ is a monomial matrices.) Note that flipping gauge spins, even without permuting them, requires a simultaneous permutation of matter spins.

The automorphism pair $(L,R)$ directly leads to the local $\mathbb{Z}_{2}$ gauge symmetry of the full lattice Hamiltonian (2). Consider an elementary plaquette $p$ depicted in Fig. 2(a) and define the local gauge transformation

[TABLE]

where $\mathcal{L}^{(\mu)}_{s}$ denotes the operator that permutes and flips the matter spins at each corner site $s$ of the plaquette as in Eq. (5): $\mathcal{L}^{(\mu)}_{s}\mu^{\rm z}_{a}\;(\mathcal{L}^{(\mu)}_{s})^{-1}=\sum_{b}\mu^{\rm z}_{a}\;(L^{-1})_{ba}$, as in the transformation in Eq. (3), and similarly for $\mathcal{R}_{s}^{(\sigma)}$ on the gauge spins: $\mathcal{R}_{s}^{(\sigma)}\sigma^{\rm z}_{i}\ (\mathcal{R}_{s}^{(\sigma)})^{-1}=\sum_{j}R_{ij}\sigma^{\rm z}_{j}$. $\mathcal{L}^{(\mu)}_{s}$ is uniquely determined by the local operator $\mathcal{R}^{(\sigma)}_{s}$ that flips the two gauge spins on links emanating from the site $s$ – just as $L$ is determined by $R$. Since here we restrict $R$ to be diagonal, corresponding to only flipping the gauge spins without permuting them, the spin flip $\sigma^{\rm z}\rightarrow-\sigma^{\rm z}$ is simply generated by $\sigma^{\rm x}$. Therefore we have $\prod_{s\in p}\mathcal{R}_{s}^{(\sigma)}=\prod_{i\in p}\sigma^{\rm x}_{i}$, where $i$ runs over all gauge spins in a plaquette. Any two $L$ matrices commute and therefore the plaquette operators do as well, $[G_{p},G_{p^{\prime}}]=0$.

The importance of $G_{p}$ is that it is a local symmetry of the full lattice Hamiltonian (2): $[H,G_{p}]=0$, for all $p$. Invariance of the Ising interaction term follows from the automorphism above, while invariance of the transverse field terms $\Gamma$ and $\widetilde{\Gamma}$ follows from two observations. First, all spin flips by the operator pair $(\mathcal{L}^{(\mu)}_{s},\mathcal{R}^{(\sigma)}_{s})$ can be viewed as $180^{\circ}$ rotations around the $x$-axis, which commute with $\sigma^{\rm x}$ and $\mu^{\rm x}$. Second, the transverse fields are uniform and therefore independent of permutations. Therefore, the Hamiltonian (2) is a gauge theory with a local $\mathbb{Z}_{2}$ gauge symmetry that is generated by $G_{p}$. This symmetry relies on the locking of the permutations contained in the operators $\mathcal{L}^{(\mu)}_{s}$ to the $\mathbb{Z}_{2}$ transformations in the $\mathcal{R}^{(\sigma)}_{s}$, which is another reason that we refer to it as combinatorial gauge symmetry.

One can further construct loop or closed string symmetry operators on the lattice, as shown in Fig. 2(b). For systems with boundaries, one can also associate a symmetry operation to open strings, as depicted in Fig. 2(c). The loop (or string) operator along a path is composed of both the gauge spin flips $\prod_{\ell}\sigma^{\rm x}_{\ell}$, where $\ell$ are the links along the path, as well as the corresponding operations on matter spins $\prod_{s}\mathcal{L}^{(\mu)}_{s}$ applied to each star along the path. In the case of closed paths, the loop operator is equivalent to a product of all plaquette operators $G_{p}$ enclosed by the loop.

Special case: ${\mathbb{Z}}_{2}$ gauge theory – The Hamiltonian Eq. (2) obeys a local ${\mathbb{Z}}_{2}$ gauge symmetry for all values of the parameters $J,\ \Gamma$, and $\widetilde{\Gamma}$. Here we shall obtain, as a particular limit, an effective Hamiltonian with a 4-spin interaction on a star, which lands directly onto the more familiar ${\mathbb{Z}}_{2}$ gauge theory on the square lattice Wegner (1971); Kogut (1979), in the following manner.

Isolate a single star with its 4 spins $\mu$ on the site and 4 gauge spins $\sigma$ on the links, as depicted in Fig. 1(a). Let us freeze for the moment a given configuration of the gauge spins $\sigma^{\rm z}_{i},i=1,2,3,4$ in the $z$-basis. The Hamiltonian (2) for each matter spin $\mu_{a}$ on a star can be viewed as that of a single spin in a magnetic field, whose eigenvalues are functions of $\sigma^{\rm z}_{i}$:

[TABLE]

The expression in Eq. (7) can be written, for any value of $\Gamma$ and $J$, as

[TABLE]

where $C_{0}$, $C_{2}$ and $C_{4}$ are constants that depend on $J$ and $\Gamma$. This expression follows from expanding the square root in Eq. (7) in powers of the $\sigma^{\rm z}_{i}$ and using $(\sigma^{\rm z}_{i})^{2}=1$ and $(W_{ai})^{2}=1$; the binary polynomial inside the square root terminates and the only terms that remain are of the form in Eq. (8). While the expansion is useful in proving the identity between Eqs. (7) and (8), we remark that the result is exact (non-perturbative), because both expressions only take values in discrete sets.

The low energy manifold of states corresponds to the sum over the lowest eigenvalues, $H_{\rm eff}^{\rm star}=\sum_{a=1}^{4}E_{a}^{\,(-)}$, which is separated from the next levels by a gap of size at least $2|\Gamma|$. We thus arrive at the following simple effective Hamiltonian for a single star:

[TABLE]

where the coefficients $\gamma$ and $\lambda$ are functions of $\Gamma$ and $J$ are explicitly given in the Supplemental Material. These relations follow from the consistency between Eqs. (7) and (8). The parity $P\equiv\sigma^{\rm z}_{1}\,\sigma^{\rm z}_{2}\,\sigma^{\rm z}_{3}\,\sigma^{\rm z}_{4}$ for the ground state of Eq. (9) is $P=+1$, since $\lambda>0$. By modifying the matrix $W$, we could flip the sign of $\lambda$ and have instead the $P=-1$ parity sector as the ground state (for example, by flipping the sign of any one column of $W$).

Let us now turn to the low energy effective model for the whole lattice. In the limit $|\Gamma|\gg J$, we find the effective Hamiltonian

[TABLE]

This Hamiltonian is exactly that of the $\mathbb{Z}_{2}$ quantum gauge theory, which supports a topological phase for $\widetilde{\Gamma}/\lambda$ below a threshold. To get the toric/surface code limit, one only has to notice that the lowest order term that survives in a perturbation theory in $\widetilde{\Gamma}/\lambda$ is the term that flips all spins around a plaquette Fradkin and Susskind (1978); Fradkin (2013); Sachdev (2018).

Taking $|\Gamma|\to\infty$, while keeping $\lambda$ fixed, opens an infinite gap to the excited sectors, where at least one $E_{a}^{(-)}$ becomes $E_{a}^{(+)}$. The splitting $2|\lambda|$ between the two parity states within the lowest energy sector remains finite. The expansion of $\lambda$ in the regime of $J\ll\Gamma$ yields $\lambda=12J^{4}/\Gamma^{3}+\mathcal{O}(J^{6}/\Gamma^{5})$. (Note that terms of order $\Gamma$ vanish.) To access this regime we would fix $\lambda$ and tune $J$ such that $J=\left|\lambda\,\Gamma^{3}/12\right|^{1/4}$. Physically, in this limit the matter fields $\mu$ can be “integrated out” to obtain the exact four-spin effective Hamiltonian.

We corroborate the above analytical features with numerical studies in the Supplemental Material. All degeneracies are confirmed to machine precision.

Combinatorial gauge symmetry for both electric and magnetic loops – So far we used the combinatorial gauge symmetry to construct a model with $\mathbb{Z}_{2}$ plaquette operators. Here we shall construct a model with both $\mathbb{Z}_{2}$ plaquette and star operators, just as in the toric code, but still using only at most two-body interactions.

We add another four spin-1/2 degrees of freedom to the center of all plaquettes in addition to the ones on the star. We denote these additional matter spins on the dual lattice as $\tau$, shown in Fig. 3(a). Furthermore, the pairwise interaction couples $\tau^{\rm x}$ and $\sigma^{\rm x}$, i.e. an $XX$ interaction. The full Hamiltonian is:

[TABLE]

In other words, on the dual lattice the spin components are transformed by $X\leftrightarrow Z$ relative to the original lattice. There is no need for a transverse field on the gauge $\sigma$ spins in this model; quantum dynamics is already present through the presence of both $XX$ and $ZZ$ interactions.

By analogy with the plaquette operators $G_{p}$ in Eq. (6), there is a set of star operators $F_{s}$, according to combinatorial gauge symmetry, which exist on the dual lattice [see Fig. 3(b)]:

[TABLE]

The dual, “left” operators $\mathcal{L}^{(\tau)}_{p}$ flip $\tau$ spins in the $x$-basis just like the operators $\mathcal{L}^{(\mu)}_{s}$ in Eq. (6) flip $\mu$ spins in the $z$-basis. By construction, these two operators commute: $[\mathcal{L}^{(\mu)}_{s},\mathcal{L}^{(\tau)}_{p}]=0$. Therefore, we have a star and a plaquette operator that also commute: $[G_{p},F_{s}]=0$, exactly as in the toric code. It is easy to check that the Hamiltonian commutes with both stars and plaquettes: $[H,G_{p}]=[H,F_{s}]=0$.

Given the commuting set of star and plaquette operators, the Hamiltonian in Eq. (11) is equivalent to the toric code in the asymptotic limit of large $\Gamma$, except that it contains only two-body interactions and fields. This is a direct result of the combinatorial gauge symmetry.

Extension to other topological states – Fracton topological phases Chamon (2005); Bravyi et al. (2011); Haah (2011); Vijay et al. (2016) (for a review, see Ref.  Nandkishore and Hermele, 2019) are novel phases of matter with a robust sub-extensive ground state degeneracy and with excitations that are strictly immobile, or constrained to move within a subdimensional manifold. Apart from theoretical interest such as classifications of phases of matter and formulations in terms of higher-rank gauge theories Pretko (2017), fracton systems are also believed to hold promise for fault-tolerant quantum computation, as well as robust quantum memory Haah (2011). In spite of the intensive theoretical investigations on fractonic models, experimental realizations directly in terms of spins have barely been discussed You and von Oppen (2018).

The building blocks of our $\mathbb{Z}_{2}$ gauge theory can also be used to construct 3D models, such as one of the simplest fractonic model, the X-cube Castelnovo et al. (2010); Vijay et al. (2016). The construction with matter and gauge spins parallels closely that in 2D, and we provide details for the construction of both the 3D toric code and the X-cube model in the Supplemental Material.

Summary and outlook – We have argued that many-body spin states admit a combinatorial gauge symmetry and we have used it to construct quantum spin liquids out of only two-body and single-body terms. The symmetry holds exactly for all ranges of parameters in the Hamiltonians that we have constructed. This presents an alternative path to explore quantum spin liquids in systems without four-body (or higher) interaction terms. Our approach may prove useful in the quest for topological qubits (via surface codes), which can potentially be hosted by certain quantum spin liquids.

Acknowledgments – The work by C. C. and Z.-C. Y. is supported by the DOE. The part of the work centered on topological phases of matter is supported by DOE Grant No. DE-FG02-06ER46316; the part of the work centered on quantum information science is supported by Grant No. DE-SC0019275.

Supplemental Material for ”Constructing Quantum Spin Liquids Using Combinatorial Gauge Symmetry”

In this Supplemental Material, we discuss: (1) large $J$ limit of Hamiltonian Eq. (2) in the main text; (2) numerical studies of the gauge-matter Hamiltonian; (3) construction of 3D toric code; and (4) construction of the X-cube model.

.1 Large $J$ limit of Hamiltonian Eq. (2)

The classical limit of Hamiltonian (2) has finite $J$ and no transverse fields $\Gamma=\widetilde{\Gamma}=0$. It is straightforward to verify that the ground state of each star on the lattice is eightfold degenerate with energy $-8J$ and parity $P=+1$. The eight ground state configurations of the gauge and matter spins on each star are shown in Table 1.

The energy gap is $4J$. Note that the matter spins are “slaved” to the gauge spins.

Now apply small transverse fields $\Gamma\ll J$ and $\widetilde{\Gamma}\ll J$ to both the matter and gauge spins. In degenerate perturbation theory the lowest order term that restores the system to the classical ground state manifold is the plaquette operator $G_{p}$. It flips four gauge spins around a plaquette plus the corresponding eight matter spins at the corners, leaving the ground state manifold degenerate. The Hamiltonian (2) in this limit is composed of a star and plaquette term:

[TABLE]

$g$ is the energy scale of the plaquette term on the order of ($\Gamma^{8}\widetilde{\Gamma}^{4})/J^{11}$, corresponding to the number of spin flips to lowest order. Because the local gauge symmetry is preserved at all stages, a quantum spin liquid state should emerge in this limit. Again, it is made possible by the combinatorial symmetry which ensures that no spin ordered state is favored. The Hamiltonian (13) is reminiscent of the well-known toric code with its star and plaquette terms Kitaev (2003). The toric code also emerges as a limit of the conventional $\mathbb{Z}_{2}$ lattice gauge theory to lowest order of perturbation theory Fradkin and Susskind (1978); Fradkin (2013); Sachdev (2018).

Numerical study of the gauge-matter Hamiltonian Eq. (2)

Here we provide a numerical comparison of the effective Hamiltonian with the four-spin interaction and the Hamiltonian with the matter and gauge spins. The studies focused on clusters containing either a single star or a single plaquette surrounded by a fixed background, in the presence of the transverse field $\widetilde{\Gamma}$. We have shown in the main text that the gauge-matter Hamiltonian for a single star in the absence of $\widetilde{\Gamma}$

[TABLE]

corresponds to the effective Hamiltonian (9) at low energy, with:

[TABLE]

The analysis passes each of the following three tests in favor of a spin liquid within machine precision.

(i) Fig. 4(a)&(b): For a single star in the presence of the $\widetilde{\Gamma}$ field on the gauge spins, the eigenstates of Eq. (9) supplemented by the transverse field on the gauge spins have a one-to-one correspondence to the lowest energy sector levels of the Hamiltonian (14) supplemented by the transverse field on the gauge spins, with precisely the same degeneracies. This implies that the presence of matter spins in the Hamiltonian (14) does not break the parity symmetry among the gauge spins even with a nonzero $\widetilde{\Gamma}$.

(ii) Fig. 4(c): In the presence of $\widetilde{\Gamma}$, the ground state wavefunction of the full Hamiltonian (14) is an equal weight superposition of the eight $P=+1$ configurations and, with smaller but equal weights, the eight $P=-1$ configurations, as should be the case for the exact star operator in the effective Hamiltonian (9). This indicates that, at the single star level, our “molecule” Hamiltonian (14) replicates the quantum dynamics of an exact star operator, even when quantum fluctuations on both the matter and gauge spins are present.

(iii) Fig. 4(d): the ground state energy of a single plaquette, as well as the gap to the first excited state, is identical for any fixed configuration of external legs surrounding the plaquette. We have confirmed that this holds for all 128 configurations of external legs, implying that the transverse field $\widetilde{\Gamma}$ does not favor any ordering pattern on the lattice. Moreover, the gap to the first excited state is independent of the environment, meaning that the effective plaquette operator generated by $\widetilde{\Gamma}$ does not favor any one spin order either.

The spectrum of the single star of the effective Hamiltonian Eq. (9) in the presence of a transverse field can be obtained straightforwardly, leading to the following 16 eigenvalues (degeneracies): $\pm\sqrt{\lambda^{2}+16\widetilde{\Gamma}^{2}}~{}(\times 1),~{}\pm\sqrt{\lambda^{2}+4\widetilde{\Gamma}^{2}}~{}(\times 4),~{}\pm\lambda~{}(\times 3)$. The constant shift $\gamma$ is not shown for simplicity. The lowest eigenstate is an equal superposition of all 8 spin configurations that satisfy $P=+1$ and, with smaller amplitude that varies with increasing $\widetilde{\Gamma}$, an equal superposition of all 8 spin configurations that satisfy $P=-1$. Symmetry of the ground state of the molecule is necessary, but not sufficient, for symmetry in the lattice.

Diagonalization of the Hamiltonian Eq. (14) in the presence of the transverse field $\widetilde{\Gamma}$ yields states with the exact same degeneracies as those of the effective model on a transverse field presented above. Fig. 4(a) and (b) present the spectrum for both the effective Hamiltonian with the four-spin interaction and the lowest energy sector for the Hamiltonian with the matter and gauge spins.

To compare the weights of the ground state wavefunction on each of the gauge spin configurations, we take the ground state $|\psi\rangle$ of the Hamiltonian with gauge and matter spins and obtain the reduced density matrix by tracing over the matter spins: $\rho_{\bm{\sigma}}={\rm tr}_{\bm{\mu}}|\psi\rangle\langle\psi|$, where ${\bm{\mu}}$ and ${\bm{\sigma}}$ stand for the four matter and gauge spins, respectively. Then we compute the weight of $\rho_{\bm{\sigma}}$ on each of the sixteen (eight with $P=+1$ and eight with $P=-1$) configurations $|{\bm{\sigma}}_{\ell}\rangle$, $\ell=1,2,\ldots,16$:

[TABLE]

These weights are displayed in Fig. 4(c). The data show that the ground state wavefunction contains an equal amplitude superposition of the eight $P=+1$ gauge spin configurations and, with smaller weight, the eight $P=-1$ configurations, mirroring exactly the case for the star operator in the effective Hamiltonian with the four-spin interaction. This indicates that, at the single star level, our “molecule” Hamiltonian (14) replicates the quantum dynamics of an exact star operator Eq. (9), even when quantum fluctuations on both the matter and gauge spins are present via the transverse fields.

To further support our claim that the system does not favor any ordered state when placed on the full lattice, we move to the next level of complexity and focus on a single plaquette. Consider a single plaquette surrounded by 8 external links, as depicted in the inset of Fig. 4(d). Suppose we fix the environment as defined by the gauge spins on the external legs. In the absence of the transverse field $\widetilde{\Gamma}$, the star constraint of positive parity is satisfied, and there are $2^{8-1}=128$ external leg configurations compatible with the constraint.

In the case of the effective Hamiltonian (9), there are only two allowed configurations of the free gauge spins satisfying the star constraints: $|\psi_{1}\rangle$ and $|\psi_{2}\rangle$, and they are related by the plaquette operator: $|\psi_{2}\rangle=\prod_{\square}\sigma^{\rm x}_{i}|\psi_{1}\rangle$. In the presence of a transverse field, the ground state and the first excited state of this plaquette must be the symmetric $|\psi_{S}\rangle$ and antisymmetric $|\psi_{A}\rangle$ superpositions of $|\psi_{1}\rangle$ and $|\psi_{2}\rangle$, respectively; and the energy splitting between these two states is given by the energy scale of the effective plaquette operator.

We find for the “molecule” Hamiltonian Eq. (2) that the same independence on the external leg configuration holds. In Fig. 4(d) shows the ground state and the first excited state energies corresponding to all 128 possible external leg configurations. In the inset of Fig. 4, we show explicitly three examples of external leg configurations. We find that the energies of both the ground state and first excited states are exactly the same for all 128 configurations, within machine precision. This is compelling evidence that the transverse field $\widetilde{\Gamma}$ does not energetically favor any specific ordering pattern on the lattice. In other words, the effective plaquette operator generated by the transverse field $\widetilde{\Gamma}$ does not favor any one spin order in particular, as it should in a spin liquid.

3D toric code

Following the path that we laid out in two dimensions, we can use Hamiltonian (2) to generate the toric code in three dimensions, as well. The main difference is that in 3D, we shall put the matter spins at the center of each plaquette (square face) of the cubic lattice, rather than on the vertex, as shown in Fig. 5(a). This will enable us to construct the exact plaquette operator in the $\sigma^{\rm z}$-basis, first.

The Hamiltonian for a single plaquette is still given by Eq. (14), for which we have shown that the ground state sector takes the form of the product of four $\sigma^{\rm z}$’s around the plaquette, as in Eq. (9). This yields the exact plaquette term of the 3D toric code. Adding a transverse field $\widetilde{\Gamma}$ on the gauge spins and going to the full cubic lattice, we arrive at the effective Hamiltonian on the full lattice:

[TABLE]

When $\widetilde{\Gamma}/{\lambda}$ is small, the leading order term acting within the ground state manifold happens at the sixth order, which involves the product of six $\sigma^{\rm x}$’s on the links emanating from a single vertex. This term is precisely the star operator of the 3D toric code, as shown in Fig. 5(b). Therefore, the same construction can implement the toric code model in both two and three dimensions.

X-cube model

The X-cube Hamiltonian contains, at each vertex, three star terms associated with three intersecting planes. Fig. 6(a) illustrates the labeling of the 6 spins at the edges of a vertex. The three stars correspond, each, to a product of 4 spins: $B_{s}^{xy}=\sigma^{\rm z}_{1}\sigma^{\rm z}_{2}\sigma^{\rm z}_{4}\sigma^{\rm z}_{5}$, $B_{s}^{yz}=\sigma^{\rm z}_{2}\sigma^{\rm z}_{3}\sigma^{\rm z}_{5}\sigma^{\rm z}_{6}$, and $B_{s}^{xz}=\sigma^{\rm z}_{1}\sigma^{\rm z}_{3}\sigma^{\rm z}_{4}\sigma^{\rm z}_{6}$.

We implement each of the three star operators using the same scheme we used in the toric code. We group the 6 gauge spins in (overlapping) sets of 4 spins, forming four-legged stars. The three groups are $(1245),(2356)$ and $(1346)$, matching the groups in the operators $B_{s}^{xy},B_{s}^{yz}$ and $B_{s}^{xz}$ above. For each of the different directions, we need a set of 4 matter spins, thus we require 12 matter spins per vertex or site of the cubic lattice. Without the transverse field $\tilde{\Gamma}$ on the 6 gauge spins, the three directions are decoupled, and it thus follows directly from the construction in the previous section that the low energy effective Hamiltonian takes the form

[TABLE]

where the coefficients $\gamma$ and $\lambda$ are given in terms of $J$ and $\Gamma$ by Eq. (15). The ground state configuration of each of three star operators $B_{s}^{xy},B_{s}^{yz}$ and $B_{s}^{xz}$ has positive parity, $P^{xy}=P^{yz}=P^{xz}=+1$, since $\lambda>0$. (Notice that there is a constraint that the product $P^{xy}\,P^{yz}\,P^{xz}=+1$.)

Paralleling the discussion for the 2D $\mathbb{Z}_{2}$ gauge theory, there is a regime where $|\Gamma|\to\infty$ while keeping $\lambda$ fixed, which opens an infinite gap to the excited sectors. The splitting $2|\lambda|$ between the two parity states within the lowest energy sector remains finite. To access this regime we would fix $\lambda$ and tune $J=\left|\lambda\,\Gamma^{3}/12\right|^{1/4}$ just as in the 2D case.

When $\widetilde{\Gamma}=0$, the spectrum of the gauge-matter Hamiltonian (with six gauge spins and three sets of four matter spins) depends only on the sum of the parities $P^{xy}+P^{yz}+P^{xz}$. This is a manifestation of the combinatorial gauge symmetry we identified in the previous section, and this symmetry is again essential to suppressing any interaction that favors any type of order.

So far we have obtained the exact star operators of a single vertex in the X-cube model. Now consider the full cubic lattice in three dimensions where the gauge spins reside on the links and the matter spins reside on the vertices. We apply a transverse field $\widetilde{\Gamma}$ to the gauge spins.

First, consider the limit $|\Gamma|\rightarrow\infty$, $J=|\lambda\Gamma^{3}/12|^{1/4}$ with $\lambda$ fixed (i.e., projecting down to the ground state sector). In this limit, the Hamiltonian on the entire lattice becomes:

[TABLE]

In the limit where $\widetilde{\Gamma}/\lambda$ is small, the lowest order term in perturbation theory that acts within the ground state subspace is the 12-body interaction around each cube $c$: $A_{c}=\prod_{n\in\partial c}\sigma^{\rm x}_{n}$, which is the cube operator of the X-cube model. Similar to Hamiltonian (10), we expect that there is a range of small $\widetilde{\Gamma}/\lambda$ where the system is in the phase with fractonic topological order. We have thus implemented the full X-cube model on the lattice, using only two-body Ising couplings and a transverse field.

Similar to the implementation of toric code in 2D, we shall now provide numerical evidence indicating that the full Hamiltonian with gauge and matter spins preserves all essential symmetries of the exact star operators in the X-cube model for a wide range of parameters $(\widetilde{\Gamma},\Gamma,J)$. Since there are now three star operators in the X-cube model satisfying the constraint $P^{xy}P^{yz}P^{xz}=+1$, the ground states must have parity +1 for all three star operators, or equivalently, $P^{xy}+P^{yz}+P^{xz}=3$; the excited states are created by flipping the parities of two out of the three stars, i.e. $P^{xy}+P^{yz}+P^{xz}=-1$. In Fig. 6(b)&(c), we plot the eigenenergy spectrum of the ground state sector where $P^{xy}+P^{yz}+P^{xz}=3$, in the presence of a transverse field $\widetilde{\Gamma}$ on the gauge spins. Once again, the levels of the full Hamiltonian [Fig. 6(b)] have a one-to-one correspondence to the levels of the exact star operators [Fig. 6(c)] for nonzero $\widetilde{\Gamma}$, with precisely the same degeneracies. While the levels of the excited states are quite complicated to count directly, transitions to the excited states can be captured again by looking at the weight of the ground state wavefunction on every classical configuration, including excited states with $P^{xy}+P^{yz}+P^{xz}=-1$, just like in Fig. 4 for the 2D toric code. This is shown in Fig. 6(d). It is clear that the ground state wavefunction has equal weights on each of the 16 configurations with $P^{xy}+P^{yz}+P^{xz}=3$ (ground states) and the 48 configurations with $P^{xy}+P^{yz}+P^{xz}=-1$ (excited states), as should be the case for the exact star operators of the X-cube model. While we are unable to perform exact diagonalizations of a single cube surrounded by fixed external legs, as we did in Fig. 4(d) for the toric code in 2D, our numerical results suggest that we indeed replicate a single vertex in the X-cube model even when quantum dynamics is introduced.