Engineering Topological Models with a General-Purpose Symmetry-to-Hamiltonian Approach

Symmetry is at the heart of modern physics. Phases of matter are classified by symmetry breaking, topological phases are characterized by non-local symmetries, and point group symmetries are critical to our understanding of crystalline materials. Symmetries could then be used as a criterion to engineer quantum systems with targeted properties. Toward that end, we have developed a novel approach, the symmetric Hamiltonian construction (SHC), that takes as input symmetries, specified by integrals of motion or discrete symmetry transformations, and produces as output all local Hamiltonians consistent with these symmetries (see github.com/ClarkResearchGroup/qosy for our open-source code). This approach builds on the slow operator method [PRE 92, 012128]. We use our new approach to construct new Hamiltonians for topological phases of matter. Topological phases of matter are exotic quantum phases with potential applications in quantum computation. In this work, we focus on two types of topological phases of matter: superconductors with Majorana zero modes and $Z_2$ quantum spin liquids. In our first application of the SHC approach, we analytically construct a large and highly tunable class of superconducting Hamiltonians with Majorana zero modes with a given targeted spatial distribution. This result lays the foundation for potential new experimental routes to realizing Majorana fermions. In our second application, we find new $Z_2$ spin liquid Hamiltonians on the square and kagome lattices. These new Hamiltonians are not sums of commuting operators nor frustration-free and, when perturbed appropriately (in a way that preserves their $Z_2$ spin liquid behavior), exhibit level-spacing statistics that suggest non-integrability. This result demonstrates how our approach can automatically generate new spin liquid Hamiltonians with interesting properties not often seen in solvable models.


I. INTRODUCTION
Symmetry is central to our understanding of the phases of matter seen in nature. Many phase transitions, such as those of liquids, magnets, or superconductors, can be described by the spontaneous breaking of symmetry according to Landau's theory [1]. Moreover, the existence of space group and point group symmetries in crystalline phases of matter highly influences the formation of order in these systems. Even topological phases of matter, exotic quantum phases that are notable for their lack of order, have non-trivial topological symmetries that give rise to their exotic properties.
We propose a new approach for studying quantum phases of matter based on symmetries. Our new approach, the symmetric Hamiltonian construction (SHC), is an algorithm that takes as input symmetries and produces as output Hamiltonians that obey those symmetries. The SHC is an example of an inverse method, a method for generating models from data. Inverse methods are widely used throughout machine learning, such as in deep learning [2]. In physics, they have been used in classical systems to design interaction potentials that stabilize crystalline and magnetic order [3][4][5][6] and in quantum systems to design or reconstruct Hamiltonians from eigenstates or density matrices [7][8][9][10][11] as well as to build single-body Hamiltonians compatible with a given symmetry group [12]. The SHC algorithm extends ideas developed in the slow operator method [13] and is quite general. The symmetries provided as input can be either integrals of motion, which can generate continuous symmetries, or discrete symmetry transformations. Example symmetries include particle number conservation, SU (2) symmetry, and point group symmetries, as well as more exotic topological symmetries such as those that we consider in this work. The Hamiltonians produced as output can be interacting as well as non-interacting; and can be made to commute or anticommute with the input symmetry operators. Our numerical implementation of the SHC is publicly available as the QOSY: Quantum Operators from Symmetry Python package [14].
In this work, we use the SHC to construct new Hamiltonians for two topological systems: superconductors with Majorana zero modes and Z 2 quantum spin liquids. We engineer superconducting Hamiltonians that commute with Majorana zero mode operators whose spatial distributions are specified as input. Separately, we construct multiple new interacting Z 2 quantum spin liquid Hamiltonians on different lattice geometries that commute with topologically non-trivial Wilson loop operators.
Majorana zero modes are integrals of motion that occur in certain superconductors that exhibit the statistics of Majorana fermions [15], an exotic class of fermions that are their own antiparticles [16]. In addition to their importance to fundamental physics, Majorana fermions, which are non-Abelian anyons, have potential applications as the building blocks for qubits in fault-tolerant quantum computers [15,17]. Many experiments have attempted to observe Majorana fermions by engineering particular superconducting Hamiltonians, such as the Kitaev chain [18], that are theoretically known to host Ma-jorana zero modes [19]. To expedite the experimental search for Majorana fermions, it is desirable to expand the small list of superconducting Hamiltonians known to possess Majorana zero modes. For this reason, we apply the SHC to design new examples of such Hamiltonians. We successfully construct a large family of local, highly tunable superconducting Hamiltonians that commute with Majorana zero modes that can be distributed arbitrarily in space. Many of these Hamiltonians have the potential to be realized in experiment.
Quantum spin liquids are exotic magnets in which spins do not order even at zero temperature due to quantum fluctuations [20]. Gapped quantum spin liquids are topologically ordered, meaning that they exhibit anyonic quasiparticles, ground state degeneracy that depends on the topology of the underlying lattice, and non-local symmetries. In particular, gapped Z 2 spin liquids host Abelian anyons and a particular set of symmetries known as Wilson loops, non-local loop operators with non-trivial topological properties. Despite significant study, there are few known model Hamiltonians that exhibit the physics of Z 2 quantum spin liquids. Some exactly solvable Hamiltonians, such as the toric code [21] and related models [21][22][23][24][25][26], are often sums of commuting operators while other candidate Z 2 spin liquid Hamiltonians [27,28] are difficult to solve numerically.
Using the SHC, we find new families of Z 2 spin liquid Hamiltonians for spins on the square and kagome lattices. The Hamiltonians that we discover are not sums of commuting operators, are not frustration-free, and can possess local and non-local integrals of motion. We find that these Hamiltonians, perturbed in an appropriate way, exhibit GOE level-spacing statistics in particular quantum number sectors, suggesting that they could be non-integrable. These models provide new, interesting examples of Z 2 topological order in spin systems.
The remainder of this paper is organized as follows. Section II describes the SHC method. Section III details our construction of superconducting Hamiltonians that commute with Majorana zero modes. Section IV describes the new Z 2 spin liquids that we found by searching for Hamiltonians that commute with Wilson loops. We conclude in Section V.

II. THE SYMMETRIC HAMILTONIAN CONSTRUCTION METHOD
In this section, we describe the symmetric Hamiltonian construction (SHC) procedure for generating Hamiltonians with desired symmetries (see Fig. 1), which includes integrals of motion (Sec. II A) and symmetry transformations (Sec. II B). In Sec. II C, we describe how these calculations can be interpreted as finding the "ground states" of super-operator "Hamiltonians." Our numerical implementation of these methods is publicly available as the QOSY: Quantum Operators from Symmetry Python package [14], whose features are . We are interested in a small subspace of this space that consists of local Hamiltonians (shown in blue). Local Hamiltonians with particular symmetries, such as Hamiltonians that commute with the integral of motionÂ or anti-commute with the antisymmetryB or are invariant under the symmetry transformationÛ, are also vector spaces (shown in green, purple, and orange, respectively). The symmetric Hamiltonian construction takes as input a list of symmetries, such asÂ,B, andÛ, and produces as output the space of local symmetric Hamiltonians that obey all of these symmetries (shown in white).
briefly discussed in Appendix F.

A. Constructing Hamiltonians with desired integrals of motion
Here we present our method for constructing Hamiltonians with a desired integral of motion, i.e., a Hermitian operatorÔ that commutes with the Hamiltonian [Ĥ,Ô] = 0. This method takes in as input an integral of motionÔ and produces as output HamiltoniansĤ that commute withÔ. Our inverse method for finding these symmetric Hamiltonians is based on the slow operator method [13], a forward method for finding integrals of motion from Hamiltonians.
In the inverse method, our goal is to find a Hamiltonian that minimizes the norm of the commutator ε ≡ [Ĥ,Ô] 2 (1) subject to the constraint that Ĥ = 1 (which can always be achieved with a trivial normalization), where Ô 2 ≡ tr Ô †Ô /tr Î is the Frobenius norm andÎ is the identity operator.
In the slow operator forward method, the goal is to find integrals of motion that minimize Eq. (1) for a given Hamiltonian. Numerically, this has been done using exact diagonalization [13,[29][30][31], matrix product operators [31,32], and other tensor networks [33]. In these contexts, the Hamiltonians studied were non-integrable models and the minimal commutator norms ε discovered were often small, but not exactly zero, making the optimized operators approximate integrals of motion. Generically, we expect the inverse problem to be easier than the forward one. This is because integrals of motion can be highly non-local while physical Hamiltonians should be local. Since local operators make up a much smaller space of operators, it is much easier to search for local Hamiltonians than non-local integrals of motion. In Section IV, for particularÔ, we are able to efficiently find Hamiltonians for which ε ≈ 10 −16 in finite-size systems.
Instead of representing the operatorsĤ andÔ in Eq. (1) as matrices or tensor networks, we find it useful to expand both of these operators in a basis of operator stringsŜ 1 , . . . ,Ŝ d that span a d-dimensional space of Hermitian operators. In particular, we consider three different types of operator strings, Pauli strings, Fermion strings, and Majorana strings, to represent all possible spin-1/2 and fermionic operators.
Finally, for fermions, we also define a 4 n -dimensional basis of Majorana strings: where , and σ a ∈ {0, 1} is chosen to makeŜ a Hermitian. The opera- iĉ i is a fermion parity operator. A discussion of the properties of these three operator string bases is provided in Appendices B, C, and D.
We perform our SHC calculations directly in a basis of operator strings. The Hamiltonian (unknown) and integral of motion (known) can be written in terms of operator strings asĤ = a J aŜa andÔ = b g bŜb , respectively, where J a (unknown) and g b (known) are real coupling constants. Our goal is to find J a such that the commutator norm is zero: ε = 0. The commutator between two operator strings can be expanded in terms of other operator strings (5) where f c ab = −f c ba are (basis-dependent) structure constants. Importantly, for bases of Pauli strings and Majorana strings, the structure constants f c ab are highly sparse and easy to compute algebraically without representing the operator strings as matrices or tensor networks (see Appendices B and D for details). Computing the structure constants for Fermion strings, however, is not efficient, so we instead perform computations in the Majorana string basis and convert back and forth to the Fermion string basis as needed.
Using the structure constants, we define the Liouvillian matrix (LÔ) ca ≡ b g b f c ab , which describes how operator strings commute with the known integral of motion: Finally, using the Liouvillian matrix, we can define the commutant matrix CÔ [29], the central quantity that we will work with in the inverse method [34]: The commutator norm can then be written as a quadratic form involving the commutant matrix where J is the vector of coupling constants of the Hamiltonian. The commutant matrix CÔ is Hermitian and positive semi-definite, making its eigenvalues real and non-negative. A normalized eigenvector J of CÔ with eigenvalue ε corresponds to a HamiltonianĤ = a J aŜa whose commutator norm withÔ is [Ĥ,Ô] 2 = ε. This indicates that the operatorsĤ andÔ exactly commute when J is a null vector, i.e., an eigenvector with zero eigenvalue (ε = 0), of the commutant matrix. Therefore, a null vector J corresponds to a HamiltonianĤ withÔ as an integral of motion. Moreover, since any linear combination of null vectors is also a null vector, we see that the null space of CÔ corresponds to an entire vector space of Hamiltonians with the desired integral of motionÔ. We now can see that the inverse method for finding Hamiltonians with a desired integral of motion amounts to finding the null space, or zero modes, of the commutant matrix. As we will discuss in Sec. II C, computing the null space of the commutant matrix can be done in a number of ways. The simplest way to do this is to explicitly construct the commutant matrix and diagonalize it numerically for a finite-dimensional basis of operator strings. This is feasible when searching for Hamiltonians since they can be represented as local operators, which can be spanned by relatively low-dimensional spaces of local operator strings.
For a fixed d-dimensional space of operators, it is interesting to consider the possible dimensionality of the commutant matrix's null space. If the null space of CÔ is one-dimensional, then there is a unique Hamiltonian H in the chosen space that commutes withÔ. If the null space dimension is greater than one, then there are many HamiltoniansĤ 1 ,Ĥ 2 , . . . in that space that, in any linear combination, commute withÔ. If the null space dimension is zero, then there is no Hamiltonian that exactly commutes withÔ in the chosen space. Nonetheless, since the eigenvalues of CÔ correspond to commutator norms, the smallest eigenvalue eigenvector of CÔ corresponds to the Hamiltonian in the space that is "closest" to commuting withÔ [35].
We also generalize this inverse method to construct Hamiltonians with desired antisymmetries, i.e., Hamilto-niansĤ that anti-commute with a desired operatorÔ. To find Hamiltonians that anti-commute withÔ, we look for Hamiltonians that minimize the anti-commutator normε whereCÔ =L †ÔLÔ is the anti-commutant matrix, abŜ c . Similar to the method described above, finding Hamiltonians that anti-commute withÔ amounts to finding the null space of the anti-commutant matrixCÔ.
Finally, we note that the inverse method described in this section is directly related to the recently developed eigenstate-to-Hamiltonian construction (EHC) algorithm [7][8][9] for constructing Hamiltonians from eigenstates. In fact, as pointed out in Ref. 8, if the integral of motion corresponds to a pure state density matrix,Ô = |ψ ψ|, then the commutant matrix of Eq. (6) is proportional to the quantum covariance matrix ψ|Ŝ aŜb |ψ − ψ|Ŝ a |ψ ψ|Ŝ b |ψ , which is the central quantity computed in the EHC algorithm. Moreover, if the integral of motion is a mixed state density matrix, then the SHC method is closely related to the method described in Ref. 10  Next, we detail how to construct Hamiltonians that are invariant under desired discrete symmetry transformations, i.e., unitary operatorsÛ g associated with a finite group g ∈ G that leave the Hamiltonian invariant U gĤÛ −1 g =Ĥ. Since [Ĥ,Û g ] = 0, in principle the method from Sec. II A could be applied. However, generically,Û g will be a sum of many operator strings, making such calculations intractable in the basis of operator strings. Instead, in this section, we describe two alternative efficient approaches that take as input a finite symmetry group G made of symmetry transformations, such as space group, time-reversal, and charge-conjugation symmetries, and produce as output Hamiltonians invariant under the action of these transformations.
In the first approach, we directly symmetrize our basis of operator strings to produce a new basis of symmetric operators [12] where g ·Ŝ a is the group action of the element g on the operator stringŜ a , which maps that operator string to a linear combination of other operator strings. This is similar to a symmetrization procedure performed by Ref. 12. For finite groups of order |G| and bases of operator strings of dimension d, we can enumerate over all operator strings and perform this calculation in time d|G|. When performing such calculations, one needs to take care to ignore non-symmetrizable operator strings [36] and not include linearly dependent symmetrized operators. In this new symmetrized basis, any linear combi-nationĤ = a J aŜ a of the symmetrized operators is a Hamiltonian that is invariant under the transformations in the symmetry group G.
In the second approach, we analyze the spectrum of the representations of the generators of the symmetry group G in the space of operator strings. In particular, consider an element g ∈ G that can be represented by a unitary (or anti-unitary) operatorÛ g acting on the usual Hilbert space of states. The action of the symmetry transformation g on an operator string is given by conjugation withÛ g and can be expanded in terms of other operator strings: The matrix D g is the representation of the symmetry transformation g on the space of operator strings and for common symmetry operators, such as space group symmetries, particle-hole symmetry, and time-reversal symmetry, can be straight-forwardly computed as we discuss in Appendix E. From Eq. (10), we see that a Hamiltonian H = a J aŜa transforms under the symmetry as [37] g ·Ĥ =Û gĤÛ From Eq. (11), we see that g ·Ĥ = ±Ĥ when J is an eigenvector of D g with eigenvalue ±1. The (+1)eigenvalue eigenvectors correspond to the coupling constants of Hamiltonians with g as a symmetry so that [Ĥ,Û g ] = 0. Likewise, the (−1)-eigenvalue eigenvectors correspond to Hamiltonians with g as an antisymmetry so that {Ĥ,Û g } = 0. Therefore, we see that we can find Hamiltonians that are invariant under a desired symmetry transformation g by finding the (+1)-eigenvectors of the representation D g , and likewise for antisymmetries the (−1)-eigenvectors of D g . To find Hamiltonians that are invariant under all of the symmetry transformations in the symmetry group G, one can compute the intersection of the (+1)-eigenspaces of the representations of the group's generators.
C. Constructing Hamiltonians with desired symmetries by finding the ground states of superoperators In the methods described above in Secs. II A and II B, we perform calculations in a Hilbert space of Hermitian operators. The commutant matrix CÔ and representation D g can be interpreted as superoperators acting on operators in this space. From this perspective, we can frame the inverse problem of constructing Hamiltonians with a set of desired symmetries and antisymmetries as finding the ground state of a particular Hermitian, positive semi-definite "superoperator Hamiltonian." In particular, if we desire to construct Hamiltonians that commute with integrals of motionÔ 1 , . . . ,Ô N1 , anticommute withÔ 1 , . . . ,Ô N2 , are invariant under symmetry transformations g 1 , . . . , g M1 and anti-invariant under g 1 , . . . , g M2 , then we can do so by finding the "ground states" of the superoperator Hamiltonian which is Hermitian and positive semi-definite by construction, where I is the identity superoperator. If we find ground states of H that have zero "energy," then we have found Hamiltonians that obey all of the desired symmetries at once. If the ground state energy is nonzero, then the energy indicates how much the discovered Hamiltonian fails to commute (or anti-commute) with the given symmetries. The simplest way to find the ground states of the superoperator Hamiltonian H is to write it as a matrix in a basis of operators and perform exact diagonalization on the matrix, e.g., using the Lanczos algorithm. Consider performing such a calculation in the Pauli string basis. Rather than working in the full 4 n -dimensional space of operators for a system of n spins, it is convenient to consider a much smaller basis of range-R k-local Pauli strings. Range-R k-local Pauli strings are Pauli strings made from a product of k (non-Identity) Pauli matrices on sites separated spatially by at most the maximum distance R. For example, the space of range-2 3-local Pauli strings on a 1D chain includes all possible spin Hamiltonians with three-site interactions between nearest and next nearest neighbor sites. In Sections III and IV, we obtain our results by exactly diagonalizing H in range-R k-local bases of Majorana strings and Pauli strings.
Finally, we note that many other well-developed methods can be used to find the ground states of H. For example, one can represent the superoperator H as a matrix product operator [32,38] and perform DMRG to find its ground state(s) or use methods such as variational Monte Carlo or other forms of quantum Monte Carlo to do so.
In attempting this, one should keep in mind that, even though the notion of locality might be kept in H, there is no known guarantee for its gapped/gapless nature, which might hinder the applicability of these methods.

III. HAMILTONIANS WITH ZERO MODES
In this section, we use the SHC to analytically design non-interacting and interacting Hamiltonians that commute with a desired pair of zero modes that can be distributed arbitrarily in space. First, we present some theoretical background on zero modes and Majorana zero modes (MZMs). Then, we describe our general framework for constructing these Hamiltonians out of "bond operator" building blocks, which are two-site operators that involve fermionic hopping, pairing, and chemical potentials. Finally, we provide examples of how this procedure can produce various Hamiltonians that commute with zero modes. We start by finding p-wave superconducting Hamiltonians that exactly commute with either exponentially-decaying MZMs, Gaussiandistributed MZMs, or zero modes with complicated spatial distributions. Then we give examples of s-wave superconducting Hamiltonians and interacting Hamiltonians that commute with MZMs. Background. A zero modeγ is a Hermitian operator that [39][40][41] whereĤ is a Hamiltonian that conserves fermion parity [Ĥ, (−1)N ] = 0, where (−1)N ≡ j (Î − 2n j ). Property 1 indicates that a zero mode is a symmetry of the system; property 2 indicates that its eigenvalues are ±1; and property 3 indicates that each state |ψ + with definite fermion parity η = ±1 comes paired with an orthogonal state |ψ − ≡γ|ψ + with opposite parity −η. For fermion-parity-conserving Hamiltonians, the opposite parity states |ψ + and |ψ − are degenerate energy eigenstates.
Generically, zero modes come in pairs. For M pairs of zero modes, the operatorsγ (1) , . . . ,γ (2M ) all commute with the Hamiltonian and satisfy the anticommutation relations These 2M zero modes lead to a 2 M -fold degeneracy in each of the eigenstates of the Hamiltonian. By pairing zero modes into complex fermionsf m ≡ (γ (2m−1) + iγ (2m) )/2 for m = 1, . . . , M , we can see that the number operatorsf † mfm commute with the Hamiltonian and each other. Therefore, these operators are simultaneously diagonalizable and the Hamiltonian eigenstates can be labeled by the occupation numbers 0, 1 for each of the M number operators [17]. From properties 1-3, we can deduce that these 2 M many-body states are degenerate in energy and that thef † m operators create singlequasiparticle modes of zero energy. Finally, we mention that Majorana zero modes (MZMs) are zero modes that are spatially localized and well-separated from one another [15]. These properties allow such zero modes to exhibit the non-Abelian statistics of Ising anyons [17].
In this work, we will consider zero modes of the form where α are real parameters that specify the "amplitudes" of the zero modes on site j (or, more generically, orbital j). To ensure thatγ (1) andγ (2) are zero modes with the properties mentioned above, the parameters are constrained such that (γ (1) ) 2 = (γ (2) ) 2 = I and {γ (1) ,γ (2) } = 0. This implies that the zero modes are orthonormal, Framework for designing zero mode Hamiltonians. We present a novel framework, illustrated in Fig. 2, for designing Hamiltonians that commute with a desired pair of zero modesγ (1) andγ (2) that can be distributed arbitrarily in space. Given a spatial distribution of the two zero modes specified by the values of α that commute with the two zero modes and are built from Hermitian bond operatorŝ which act only on sites i and j. The J ij are real parameters that independently scale each bond operator and are arbitrary up to the constraint that they must be non-zero on a connected graph of bonds {(i, j)}. This constraint guarantees that the Hamiltonians do not commute with any other zero modes that are linear combinations of the Majorana fermionsâ j andb j other thanγ (1) andγ (2) [42]. The family of Hamiltonians in Eq. (16) includes Hamiltonians on various lattices and graphs such as square lattices, kagome lattices, tetrahedral lattices, trees, and aperiodic tilings. We were able to determine the analytic form of the bond operatorsĥ ij and their properties by computing the commutant matrices Cγ(1), Cγ (2), and Cĥ ij , as described in Appendix G. The parameters of the bond operatorŝ h ij depend on the zero mode amplitudes α For the chemical potential parametersμ j , we add an (ij) superscript to make clear that these chemical potentials on sites i and j are associated with the bond operatorĥ ij and not another bond operator, such asĥ ik orĥ kj . Also, importantly, we choose an ordering convention for our fermionic operators and require that i < j for each bond operatorĥ ij to be consistent with our convention.
Each bond operatorĥ ij individually commutes witĥ γ (1) andγ (2) , though overlapping bond operators do not commute with one another: [ĥ ij ,ĥ jk ] = 0. This makes theĤ ZM Hamiltonians analogous to frustrationfree Hamiltonians, which by definition have ground states that are simultaneously the ground states of each local term in the Hamiltonian [43,44].
We can rewrite Eq. (16) aŝ whose parameters take on the values and 1.

2.
FIG. 2. Steps for constructing zero mode Hamiltonians on an arbitrary graph. 1. Specify the spatial distributions of the zero modesγ (1) andγ (2) by choosing their amplitudes α j , respectively, on the vertices j of the graph. 2. Draw a set of edges between the vertices until the vertices and edges form a connected graph. The resulting graph represents a family of Hamiltonians of the formĤZM = ij Jijĥij, where Jij = 0 on the (i, j) edges andĥij are bond operators specified by Eq. (17). These Hamiltonians exactly commute withγ (1) andγ (2) and no other zero modes of the form of Eq. (14).
The chemical potential on a site j is the sum of the chemical potentials contributed by each bond operator. The need for two sums is due to our choice of convention that i < j for each bond operatorĥ ij .
To simplify the remaining discussion in our examples, we consider the restricted class of zero modes with β where we relabeled α (1) j → α j and β (2) j → β j . One reason to consider this class of zero modes is that it contains the zero modes of the Kitaev chain [18]. Another reason is that, upon normalization, these zero modes automatically satisfy Eq. (15).
For the zero modes of Eq. (22), the bond operator can be simplified tô Example: Exponentially decaying Majorana zero modes. As a first example, we discuss how to use our framework to construct p-wave 1D chain Hamiltonians with exponentially decaying Majorana zero modes. We will see that the models we construct this way are closely related to the Kitaev chain, which has nearly exponentially decaying MZMs [18]. As input, we choose two MZMs exponentially localized at each edge of a 1D chain of L spinless fermions: where ξ > 0 is a correlation length in units of the lattice spacing. For a bond between the sites i and j separated by a distance d = |i − j|, the parameters of the bond operator Eq. (24) becomẽ The bond operators with these specific parameters define a large family of HamiltoniansĤ ZM = ij J ijĥij that exactly commute with the zero modes of Eq. (25).
First, we will focus on a simple, interesting subspace of these zero mode Hamiltonians. In particular, we consider those Hamiltonians constrained to have constant nearest neighbor (d = 1) hopping −t on a 1D chain. We implement this constraint by choosing J ij = −tδ i,j−1 /t ij for each bond (i, j). Given this constraint and our input zero modes, we find a unique Hamiltonian that commutes with the desired zero modeŝ where Note that this Hamiltonian almost matches the Kitaev chain [18] with pairing ∆/t = tanh(1/ξ) and chemical potential µ/t = −2/ cosh(1/ξ), though slightly differs at the edges j = 1, L of the chain:Ĥ Another interesting subspace of this large family of Hamiltonians is the space of Hamiltonians with d-th neighbor bonds for 1 < d < L/2. For this subspace, we choose the constraint that J ij = −tδ i,j−d /t ij . Under this constraint, we can use the bond operators to construct the following Hamiltonianŝ where Interestingly, all Hamiltonians of Eq. (28), in any linear combination, exactly commute with the exponentially decaying MZMs of Eq. (25). However, note that only the nearest-neighbor d = 1 Hamiltonian connects together all of the sites in the 1D chain, while the d > 1 Hamiltonians form disconnected graphs. Therefore, Hamiltonians of the form exp with J 1 = 0 commute with exactly two zero modes, while those with J 1 = 0 potentially commute with more than two.
Example: Gaussian-distributed Majorana zero modes. Next, we construct Hamiltonians that commute with MZMs that are spatially localized as Gaussians of width σ centered at positions x 1 and x 2 . We provide as input zero modes with the amplitudes where we replaced the site label j with its spatial coordinate x in a lattice so that α j , β j → α x , β x . For concreteness, we will consider a 1D chain lattice and 2D square lattice, but the same construction applies to arbitrary lattices in any dimension.
The α x , β x parameters determine the coefficients of the bond operatorĥ ij →ĥ x,x+δ that connects the site i at position x to site j at position x + δ. For Gaussiandistributed zero modes, the parameters of the bond operatorĥ x,x+δ of Eq. (23) satisfỹ An interesting property to note is that∆ x,x+δ /t x,x+δ is actually independent of position x and only depends on the displacement of the MZMs x 2 − x 1 and the direction of the bond δ.
By arranging these bond operators uniformly onto the nearest neighbor bonds of a 1D chain or 2D square lattice, i.e., choosing J x,y = −tδ x,y−δ /t x,y , we construct the following two Hamiltonians that commute with the two

Gaussian-distributed zero modeŝ
where the pairings and chemical potentials can be generated from Eq. (30) by applying Eqs. (20)- (21). To illustrate, we show the chemical potential distributions for a 100-site 1D chain and a 100 × 100 2D square lattice with Gaussian Majorana zero modes in Figs. 3 and 4, respectively. Again, there is nothing special about chains, square lattices, 1D, or 2D. The same Hamiltonian as Eq. (32), but in different dimensions and with different lattice vectors δ, would also commute exactly with two Gaussian-distributed MZMs. Even more generally, a similar construction for any arbitrarily connected graph is straightforward, as discussed above.
Example: Zero modes with complicated spatial distributions. Using theĥ ij bond operators, we can also design Hamiltonians with zero modes that have highly non-trivial spatial distributions.
To illustrate, we provide as input to our framework a pair of two spatially separated "Majorana" MZMs shaped according to a portrait of Ettore Majorana, shown in Fig. 5(a). As we did above, we also lay down bond operators uniformly onto a 2D square lattice so that the nearest-neighbor hopping between sites is −t. The Hamiltonian that we find, whose parameters are shown in Fig. 5(b)-(d), is both complicated and simple. It is complicated because of the non-trivial spatial distributions of the pairing and chemical potential terms, but it is also simple because of its locality and non-interacting nature.
We also construct Hamiltonians that commute with exotic zero modes with non-trivial spatial distributions that possess some but not all of the properties of MZMs. Recall that a MZM is (1) spatially localized into a single location and (2) well separated from other MZMs. Here we consider two examples of zero modes that break the first of these two properties. As we did above, we consider a 2D square lattice geometry for our Hamiltonians and require constant hopping between sites. In our first example, shown in Fig. 6(a), we provide as input a pair of zero modes that are spatially localized into two locations but still well-separated from one another. The parameters of the resulting Hamiltonian are shown in Fig. 6(b)-(c). In our second example, shown in Fig. 6(d), we provide as input a pair of well-separated zero modes in which one of the modes is spatially delocalized into a ring surrounding the other zero mode. The parameters of the Hamiltonian that commutes with the pair of zero modes are shown in Fig. 6(e)-(f). In both cases, the chemical potential distributions and pairing distributions required to produce these zero modes are quite non-trivial.
Note that for the Hamiltonians depicted in Figs. 5 and 6, we use the same notation (replacing the labels i, j → x, x + δ) as we did for the Gaussian zero modes, but the µ x , ∆ x,x+δ parameters are from the more general Eq. (24) instead of the special-case Eq. (30). Example: s-wave superconducting Hamiltonians with Majorana zero modes. While we have restricted our attention to spinless fermions with p-wave superconductivity, it is also possible to use our framework to construct spinful Majorana zero mode Hamiltonians with different superconducting order parameters. However, building these models can come at the expense of breaking symmetries, such as spinful time-reversal symmetry and spin conservation in the z-direction; or the complication of employing spin-orbit coupling or applying a local magnetic field. To add spins to our models, we substitute our labels with i, j → iσ, jσ where i, j now correspond to sites and σ, σ ∈ {↑, ↓} to spins.
Using our framework, we construct a 1D s-wave superconducting Hamiltonian with exponentially localized MZMs on the edges in two steps. (1) First, we construct two disconnected number-conserving Hamiltonians, one for the spin-up fermions and one for the spin- down fermions, that commute with four zero modes (two at each edge). (2) Second, we add a coupling between the up and down spins that does not commute with two of the four zero modes, leaving only one zero mode on each edge. Below we describe the two steps in detail.
(1) Consider the spin-up fermions. We first want to construct a Hamiltonian that commutes with a pair of spin-up exponentially decaying edges modes, both of which decay at rate 0 < r < 1 from the left edge: γ For these two zero modes, the corresponding bond operator between neighboring sites in the chain isĥ j↑,j+1↑ = (ĉ † j↑ĉ j+1,↑ + H.c.) − rn j↑ − r −1n j+1↑ . Note that because the two zero modes are decaying in the same direction this operator has no superconducting pairing, only hopping and chemical potential. Now, consider the spin-down fermions. We would like to construct another Hamiltonian that commutes with two spin-down zero modes decaying exponentially at the same rate from the right edge:γ For these zero modes, the bond operator between neighboring sites isĥ j↓,j+1↓ = (ĉ † j↓ĉ j+1,↓ +H.c.)−r −1n j↓ −rn j+1↓ . We can add the spinup and spin-down bond operators together to construct a 1D chain Hamiltonian −t L−1 j=1 σĥ jσ,j+1σ that commutes with the four zero modesγ ↓ .
(2) To ensure that only one zero mode persists at each edge, we add a perturbing bond operatorĥ j↑,j↓ = ib j↑âj↓ = −ĉ † j↑ĉ j↓ −ĉ † j↑ĉ † j↓ +H.c. that commutes withγ Altogether, the Hamiltonian that we construct iŝ where t and ∆ s are real free parameters,n j ≡n j↑ +n j↓ , µ = −t(r + r −1 ), and the edge term iŝ which involves a chemical potential and magnetic field. This Hamiltonian is an s-wave superconductor that commutes with only the two desired Majorana zero modesγ (1) ↑ andγ (2) ↓ . In this case, this Hamiltonian breaks timereversal symmetry, does not conserve z-magnetization, involves spin-orbit coupling of the same strength as the pairing, and requires a finely tuned magnetic field at the edge. This edge magnetic field, however, we do not expect to be essential and could possibly be removed by slightly modifying the spatial distributions of the input zero modes. Interacting Hamiltonians with zero modes. Finally, we mention how to construct interacting Hamiltonians that commute with particular zero modes. The main fact to note is that if Hermitian operatorsÂ and B commute with a zero modeγ (1) , then so do the Hermitian operators i[Â,B] and {Â,B} if they are non-zero. For example, for bond operatorsĥ ij andĥ jk that commute withγ (1) andγ (2) , the operator {ĥ ij ,ĥ jk } = 0 also commutes withγ (1) andγ (2) . Using this observation, we can construct the class of fermion-parity-conserving interacting Hamiltonianŝ where the coefficients c ij , d ijkl , . . . form a connected graph. These Hamiltonians often contain complicated interacting terms such asn inj ,n iĉ † As an example, consider the s-wave Hamiltonian Eq. (33) that we constructed above. To this Hamiltonian we can add an interacting term between neighboring sites of the form j {ĥ j↑,j↓ ,ĥ j+1↑,j+1↓ }/2 = − jb j↑âj↓bj+1↑âj+1↓ and still have the resulting Hamiltonian commute with exactly the same two zero modes. When written in terms of complex fermionsĉ j andĉ † j , this term is a sum of eight quartic fermionic operators, many of which do not conserve particle number.

IV. Z2 QUANTUM SPIN LIQUID HAMILTONIANS
In this section, we use the SHC to numerically construct large classes of new Z 2 quantum spin liquid Hamiltonians on the square and kagome lattices. We discover many Hamiltonians that commute with the Wilson loops shown in Fig. 7. All of these Hamiltonians have at least four-fold degenerate ground states. We perform exact diagonalization on these Hamiltonians and determine that many have exactly four-fold ground state degeneracy. For many of these Hamiltonians with four-fold ground states, we compute the modular S-matrix and find that it exactly matches the modular S-matrix of Z 2 spin liquids. Generically, the Z 2 spin liquid Hamiltonians that we find are not sums of commuting projectors, nor frustrationfree. For some of these Hamiltonians, the Wilson loops are "rigid", i.e., they are only of a fixed length. In other Hamiltonians we find that some Wilson loops can be deformed, or extended to arbitrary length, like Wilson loops in the toric code.
The Hamiltonians with deformable Wilson loops possess many local integrals of motion, while those with rigid loops do not. For none of these Hamiltonians do the integrals of motion that we identify form a complete mutually commuting set that fully specifies all of the eigenstates. While it is possible that there are integrals of motion we do not know about, we are able to explicitly rule out some types of integrals of motion. An exhaustive numerical search rules out the existence of any additional truly local integrals of motion up to some range. A novel eigenstate clustering approach, discussed in Appendix H, rules out the existence of a complete set of integrals of motion that commutes with a class of (Wilson-loop-preserving) perturbations of the Hamiltonian. Finally, we consider the level-spacing statistics of the Z 2 spin liquid Hamiltonians under these perturbations. For the energy levels of these perturbed Hamiltonians in individual quantum number sectors, we find GOE level-spacing statistics suggesting non-integrable behavior.
Background. Wilson loops arise as integrals of motion in Z 2 quantum spin liquids. Consider a square lattice of spins wrapped into a torus so that there are periodic boundaries in both directions. On the torus, we can form two topologically inequivalent non-contractible loops Lx and Lŷ that span the entire system. Consider four nonlocal Wilson loop operatorsX Lx ,X Lŷ ,Ẑ Lx andẐ Lŷ with non-trivial support along these two loops. Suppose that these Wilson loop operators are integrals of motion of a HamiltonianĤ, that they square to identity (so that their eigenvalues are ±1), and that they obey the following set of commutation and anti-commutation relations The existence of Wilson loops satisfying these properties implies at least a four-fold degeneracy in each of the energy eigenstates of the HamiltonianĤ.
In this work, we will focus on a particularly simple set of Wilson loops of the form where Lx and Lŷ are two straight-line topologically distinct non-contractible loops across the torus that wind in the horizontal and vertical directions, respectively (see Fig. 7(a)). It can be verified that the Wilson loops of Eq. (36) square to identity and satisfy the properties of Eq. (35). These particular Wilson loops are integrals of motion of the toric code, a well known solvable Z 2 spin liquid Hamiltonian that is a sum of commuting terms [23]. Note that there are alsoŶ L ∝X LẐL Wilson loops, though they are not independent from the ones in Eq. (36). Also, there are actually many straight-line Wilson loop operators whose loops are parallel to Lx and Lŷ, though shifted by jŷ or kx. We will refer to these shifted loops as L (j) x and L (k) y for j, k = 1, . . . , L. In addition to these Wilson loops on the square lattice, we also consider straight-line Wilson loops of the same form as Eq. (36) on the periodic kagome lattice. We refer to the kagome Wilson loops asẐ La 1 ,X La 1 ,Ẑ La 2 ,X La 2 , where a 1 = (1, 0) and a 2 = (1/2, √ 3/2) are the kagome lattice vectors (see Fig. 7(c)).
Finally, we would like to mention the form of the toric code on the kagome lattice, as it will be relevant to our later discussion. The toric code model can be defined on many lattices, such as the honeycomb lattice. While it is customary for the spins of the toric code to be on the links of the lattice, in this work we always consider spins to be on the sites. When the spins are defined on sites instead of links, the honeycomb toric code gets mapped to the kagome lattice. This model, which we will refer to as the kagome toric code model, iŝ where theX = j∈ σ x j andẐ = j∈ σ z j operators are three-spin and six-spin interactions defined on the triangles and hexagons, respectively, of the lattice. The model has the same essential features as the square lattice toric code: it is a sum of commuting terms, theX ,Ẑ are local integrals of motion, and the model commutes with the straight-line Wilson loopŝ Z La 1 ,X La 1 ,Ẑ La 2 ,X La 2 .
A. Z2 spin liquid Hamiltonians on the square lattice SHC numerics. To generate new Z 2 spin liquid Hamiltonians on the square lattice, we provided a list of desired symmetries as input to SHC: (1) the four straight-line Wilson loop operatorsX Lx ,Ẑ Lx ,X Lŷ ,Ẑ Lŷ (see Fig. 7(a)); and (2) the symmetry group of the square lattice, generated by translations of lattice vectorsx,ŷ, 90 • -rotation, and reflection about the side of a square. Note that using theX,Ẑ pair of Wilson loops is an arbitrary choice and we could have instead used theX,Ŷ orŶ ,Ẑ pair. Given this input, one might expect that the SHC would produce the toric code as output, since it is an example of a square lattice Hamiltonian that commutes with Wilson loops of this form. However, it does not because the toric code's translational symmetry is generated by translations of 2x and 2ŷ instead ofx and y.
We performed the SHC calculations on a finite-size N = 8 × 8 = 64 site lattice. We used a basis of range-R k-local Pauli strings, where R = 2 and k ∈ {1, 2, 3, 4, 5}; this includes up to five-spin interactions on nearest, nextnearest, and next-next nearest neighbor sites on the square lattice. Before symmetrization by spatial symmetries, this basis was 67, 584-dimensional. After symmetrization, the basis was reduced to a 234-dimensional basis of spatially-symmetric Hamiltonians. In this 234dimensional basis, we numerically constructed the commutant matrices CX Lx , CẐ Lx , . . . for the four Wilson loops and found three vectors that were null vectors of all of these matrices. These three vectors correspond to the coupling constants of three Hamiltonians with all of the desired symmetries. This three-dimensional space of symmetric Hamiltonians takes the form whereX = j∈ σ x j ,Ŷ = j∈ σ y j , andẐ = j∈ σ z j are four-spin interactions on the nearest neighbor squares of the lattice and ξ x , ξ y , ξ z are arbitrary real constants. These Hamiltonians are depicted in Fig. 7(b). Numerical checks of Z 2 order. While the Hamiltonians of Eq. (38) commute with the Wilson loop operators, it is not guaranteed that they are Z 2 spin liquids. In this section, we numerically tested particular Hamiltonians in this space to check if they have Z 2 topological order. By construction, these Hamiltonians commute with Wilson loops satisfying Eq. (35) and so are guaranteed to have eigenstates with degeneracies that are multiples of four. However, it is possible that they have greater than four-fold degeneracy, either due to the existence of accidental degeneracy or additional symmetries that we did not require.
For eachĤ square that we tested, we used exact diagonalization (ED) to determine if the ground state was exactly four-fold degenerate. If it was, we then calculated the modular S-matrix, a quantity that encodes the properties of anyons in a topologically ordered system [17,45]. The modular S-matrix is an overlap matrix S ij ≡ Ξx i |Ξŷ j between minimally-entangled-states (MES) |Ξx i and |Ξŷ i for i = 1, . . . , 4 across loops Lx and Lŷ of the torus. The MES are particular linear combinations of the four degenerate ground states in the system that minimizes the Renyi entanglement entropy along the cuts defined by loops Lx and Lŷ. For the toric code, the MES are e and m flux eigenstates along the two loops, i.e., they are eigenstates of the Wilson loop operatorsẐ L ,X L [45]. Assuming this holds for theĤ square Hamiltonians as well, we computed the MES for our Hamiltonians by finding the flux eigenstates ofẐ L δ ,X L δ for δ =x,ŷ from the four-fold degenerate ground states. In particular, for δ =x andŷ, we built the set of four MES |Ξ δ 1 , . . . , |Ξ δ 4 by computing the (unique) ground state ofĤ square − κ eẐL δ − κ mXL δ for (κ e , κ m ) = (+1, +1), (+1, −1), (−1, +1), (−1, −1), respectively. We numerically verified for the ξ x = ξ y = ξ z = −1 Hamiltonian on the 4 × 4 lattice that the MES that we used did in fact minimize the Renyi entanglement entropy across the Lx and Lŷ cuts of the torus, as is expected [45].
Using ED, we checked 27 particularĤ square Hamiltonians on a 4 × 4 square lattice to determine if they were Z 2 spin liquids. We considered the 3 3 = 27 possible Hamiltonians with ξ α set to −1, 0, +1 for each α = x, y, z. Six of these Hamiltonians (up to cyclic permutations of ξ x , ξ y , ξ z ), listed in Table I, had exactly four degenerate ground states. For these six Hamiltonians, we computed the modular S-matrix and found it to be [46] up to numerical precision. This S-matrix corresponds to Z 2 topological order [45]. Symmetries of the discovered Hamiltonians. Here we consider the symmetries of the Hamiltonians that we found. By construction, these models possess straightline Wilson loops as integrals of motion and obey the spatial symmetries of the square lattice. However, it is possible that these Hamiltonians possess additional symmetries that we did not require. We numerically determined that theĤ square Hamiltonians listed in Table I do not possess highly local integrals of motions. We verified this numerically by using the slow operator forward method, i.e., by computing the commutant matrices CĤ square for each of the Hamiltonians in Table I with a local basis of operators. For the Hamiltonians tested, we found no integrals of motion with up to 9-site terms on a local 3 × 3 square cluster of sites in the lattice. The fact that these Hamiltonians do not possess such local integrals of motion suggests that the Wilson loops of these models are "rigid," i.e., cannot be locally deformed.
Even without any apparent local integrals of motion, we are able to identify a set of mutually commuting integrals of motion built from the straight-line Wilson loops. Consider an L x × L y square lattice. Without local integrals of motion, theẐ L (j) x ,Ẑ L (k) y Wilson loops for j = 1, . . . , L y ; k = 1, . . . , L x − 1 are all independent conserved quantities [47]. The same is true for the X Wilson loops, but the X Wilson loops do not commute with all of the Z Wilson loops (see Eq. (35)). Nonetheless, prod-ucts of two X loops do commute with all of the Z loops. Therefore, the set of operators form a set of 2L x + 2L y − 4 mutually commuting operators. Generically, the square latticeĤ square Hamiltonians of Eq. (38) also possess a global integral of motion that is a sum of Pauli strings. These Hamiltonians can be written as a sum of two commuting operatorsĤ square =Â +B, whereÂ and A are the "black squares" and B are the "white squares" of a black-white checkboard pattern laid down on the square lattice. Since [Â,B] = 0, we can consider one of these operators, sayB, as a global integral of motion, so that [Ĥ square ,B] = 0. TheB operator also commutes with the integrals of motion listed in Eq. (40). Note that theĤ square Hamiltonians of Eq. (38) are not sums of commuting terms nor frustration-free, making them difficult to solve analytically. The set of mutually commuting integrals of motion of Eqs. (40) and (42) is not enough to fully diagonalize theĤ square Hamiltonians, but can be used to block diagonalize them into quantum number sectors, allowing us to more effectively study these models numerically. Level-spacing statistics. Typically, the level-spacing statistics of a quantum Hamiltonian are either distributed according to the Gaussian orthogonal ensemble (GOE) distribution or the Poisson distribution. If the level-spacing statistics are GOE distributed, then that is good evidence that the system is non-integrable. For Hamiltonians with many integrals of motion, such as the ones we found, the level-spacing statistics will generally appear Poisson when considering energies spread out over multiple quantum number sectors of the integrals of motion. However, it is possible for the statistics to be GOE in particular quantum number sectors.
Using ED, we numerically examined the levelspacing statistics of particularĤ square Hamiltonians, i.e., the statistics of the level-spacing ratio r n = min(δ n , δ n−1 )/ max(δ n , δ n−1 ) where δ n ≡ E n − E n−1 and E n are the sorted energy eigenvalues of the Hamiltonian. For the Poisson distribution the average level-spacing ratio is expected to be r Poisson = 0.3863, while for the GOE distribution it should be r GOE = 0.5307. After accounting for the quantum number sectors we know about -i.e., the ones listed in Eqs. (40) and (42) -we observed significant degeneracy in the spectrum, leaving only a small number of unique energy levels on the 8 × 4 square lattice that we considered. To break this degeneracy, we perturbed the sixth Hamiltonian of Table I  where δĤ = h Ẑ and h are random numbers drawn from the uniform distribution between −1 and 1 for each square and is the disorder strength. This particular perturbation breaks the spatial symmetries of the square lattice, but preserves the Wilson loop integrals of motion in Eq. (40) (and a modified global integral of motion Eq. (42)) and generically destroys the eigenstate degeneracy within the known quantum number sectors. After this perturbation, there are still eigenstates in a given quantum number sector that do not couple with each other. We cluster these eigenstates by grouping together the connected set of states that couple through the perturbation to the Hamiltonian (see Appendix H). This coupling still leaves multiple eigenstates per cluster and suggests the existence of some "hidden" integrals of motion we have not explicitly identified. For 10 random realizations of Eq. (43) for from 0 to 6, we computed the average level-spacing ratio within these clusters. As shown in Fig. 9, as increases to 2 the average levelspacing ratio r approaches the GOE value. However, it does decrease slightly below that value for larger , as shown in Fig. A1. For the 8 × 4 square lattice, the eigenstate clusters we found were quite small, typically containing either 28 or 35 states. We also tested 10 random realizations of these perturbed Hamiltonians with disorder strengths = 0.25, 0.5, 0.75, 1 on a 4 × 4 square lattice and observed that they always possessed exactly four-fold degenerate ground states, suggesting that the perturbed models are also Z 2 spin liquids.  (2) the symmetry group of the kagome lattice generated by translations of lattice vectors a 1 and a 2 , 60 • -rotation, and reflection. We performed our SHC calculations on the finite-size N = 48 site symmetric cluster shown in Fig. 7(c). To construct Hamiltonians, we used a basis of range-R 3local and 6-local Pauli strings, where R = 2/ √ 3. Before symmetrization by spatial symmetries, this basis was 31, 536-dimensional. After symmetrization, the basis was reduced to a 220-dimensional space. In this large space of local Hamiltonians, we found the following sixdimensional space of symmetric Hamiltonians that ex- actly obey all of the desired symmetrieŝ where η α , χ α for α = x, y, z are arbitrary real constants, X = j∈ σ x j ,Ŷ = j∈ σ y j , . . . , and the summations are over all of the triangles (both upward and downward facing) and hexagons in the kagome lattice. These Hamiltonians are depicted in Fig. 7(d).
We also note that we performed the same SHC calculations with a basis of range-R k-local Pauli strings with R = 2/ √ 3 and k ∈ {1, 2, 3, 4, 5, 6} and found six additional symmetric Hamiltonians, which involve four and five-site interactions. These Hamiltonians can be written in terms of products of triangle operators: where the summations are over nearest-neighbor triangles and that overlap at a single site. While we include these Hamiltonians for completeness, we will not examine their properties in the discussion below. We instead will focus on the Hamiltonians of Eq. (44). Numerical checks of Z 2 order. Using ED on finitesize kagome lattices, we found Hamiltonians of the form of Eq. (44) that exhibit Z 2 topological order. On a 3×2× 3 = 18 site kagome lattice, we consideredĤ kagome Hamiltonians with all possible 3 3 = 27 combinations of η α = −1, 0, +1 for α = x, y, z with (χ x , χ y , χ z ) = (0, 0, −1). We found that 24 of these 27 Hamiltonians had four-fold degenerate ground states and exactly the Z 2 modular Smatrix of Eq. (39), which we computed in the same way as described in the previous section. The three Hamiltonians that did not have these properties were the effectively classical (η x , η y , η z ) = (0, 0, −1), (0, 0, 0), (0, 0, +1) Hamiltonians. Symmetries of discovered Hamiltonians. By construction, all of theĤ kagome Hamiltonians of Eq. (44) obey the symmetries of the kagome lattice and commute with the straight-line Wilson loops. Yet, generiĉ H kagome Hamiltonians do not possess any highly local integrals of motion. We checked this numerically using the slow operator forward method. In particular, we computed the commutant matrix CĤ kagome for 100 random H kagome Hamiltonians, with η α , χ α sampled uniformly from [−1, 1], and found that they had no null vectors for local bases of Pauli strings on the triangles, bowties, and hexagons of a 48-site kagome lattice. For genericĤ kagome Hamiltonians, these calculations suggest that there are no local integrals of motion and that the Wilson loops are rigid like they were for the square lattice Hamiltonians. For these generic Hamiltonians, we can identify the following set of mutually commuting integrals of motion However, as we discuss below, particular subspaces of these Hamiltonians possess different sets of mutually commuting integrals of motion, which do include local integrals of motion and deformable Wilson loops. For example, the Hamiltonians in the four-dimensional subspace of Hamiltonians where χ z = 0, commute withẐ for all hexagons. For these models, we can therefore build a large set of mutually commuting operators Moreover, for these models, the Wilson loops are partially deformable: no Wilson loops can be deformed around triangles, but the Z Wilson loops can be deformed around hexagons. Similar subspaces with different hexagon operators, e.g., with χ x = 0, χ y = χ z = 0, also exist and have the same properties.
Likewise, Hamiltonians of the form where η z = 0, commute withẐ for all triangles. For these models, we can build an even larger set of mutually commuting operators since there are more triangles than hexagons in the kagome lattice. For these models, the Z Wilson loops can be deformed around triangles, but no Wilson loops can be deformed around hexagons. Interestingly, the kagome lattice toric code of Eq. (37) is a particularly special point, (η x , η y , η x , χ x , χ y , χ z ) = (−1, 0, 0, 0, 0, −1), in the space of Hamiltonians. It is special in that both X triangles and Z hexagons are local integrals of motion, making X Wilson loops deformable about triangles and Z Wilson loops deformable around hexagons (see Fig. 8). Since the toric code is an exactly solvable model, it is interesting to consider perturbing it to better understand the other models in this space of Hamiltonians, which have different symmetry properties. Fig. 10 shows the full spectra of three families of Hamiltonians that all include the toric code as a special point. The Hamiltonians of Fig. 10(a) are of the type described in Eq. (46), the Hamiltonians of Fig. 10(b) are of the type described in Eq. (48), and the Hamiltonians of Fig. 10(c) are genericĤ kagome Hamiltonians of the type described in Eq. (44). Notably, the ground state for each of the Hamiltonians shown stayed exactly four-fold degenerate, though the gap to the first-excited states did decrease with large enough perturbations.
Generically, the kagome latticeĤ kagome Hamiltonians of Eq. (44), like the square lattice models, also possess a global integral of motion that is a sum of Pauli strings. A genericĤ kagome Hamiltonian is a sum of two commuting terms,Ĥ kagome =Ĉ +D, wherê are the terms only on the triangles or hexagons of the kagome lattice, respectively. We can take one of these operators, sayD, as a global integral of motion, since [Ĥ kagome ,D] = 0. TheĤ kagome Hamiltonians of Eq. (44) are not sums of commuting terms nor frustration-free. Generically, the set of mutually commuting integrals of motion that we found do not fully specify the degrees of freedom of the model and cannot be used to fully diagonalize thê H kagome Hamiltonians, except at the special points corresponding to the toric code model. However, the set of mutually commuting integrals of motion Eq. (45) (or The full spectra of the Hamiltonians (a) H T C,kagome + α1 Ẑ , (b)Ĥ T C,kagome + α2 X , and (c)Ĥ T C,kagome + α3( Ẑ + X ) on a 3 × 2 × 3 kagome lattice. The width of each horizontal line corresponds to the degeneracy of the energy eigenstates, plotted on a log scale. The smallest degeneracy is four-fold, which occurs for all of the ground states of these models, and the largest degeneracy is about 80, 000-fold, which occurs for the eigenstates in the center of the toric code's spectrum (α1 = α2 = α3 = 0).
Eq. (47) or (49)) and Eq. (51) can be used to block diagonalize the Hamiltonians into quantum number sectors, which allows us to more effectively study these models numerically. Level-spacing statistics. We examined the levelspacing statistics of the Hamiltonians we discovered with the SHC to see if they obeyed GOE or Poisson statistics in their quantum number sectors. Initially, we computed the level-spacing statistics of the Hamiltonians in Fig. 10 for 3 × 2 × 3 = 18 and 3 × 3 × 3 = 27 site kagome lattices in the quantum number sectors of Eq. (45) (or Eq. (47) or (49)) and Eq. (51). Like we observed for the square lattice Hamiltonians, the kagome lattice Hamiltonians possess large eigenstate degeneracy. We again remove these degeneracies by considering perturbed Hamiltonians of the form where δĤ = h Ẑ and h are random numbers drawn from the uniform distribution between −1 and 1 for each triangle and is the disorder strength. These particular perturbations break the spatial symmetries of the kagome lattice, but still preserve the integrals of motion described in the previous section.
On 18-site lattices, for both Hamiltonians Eqs. (52) and (53), the δĤ perturbation was sufficient to produce GOE statistics at intermediate for the 256 energy levels in the considered quantum number sectors, as shown in Fig. 9. On the 27-site lattice (for which we only considered Hamiltonian Eq. (52)), this perturbation was not sufficient to produce GOE statistics; instead we found that there were additional unidentified "hidden" quantum number sectors that were affecting the level-spacing results.
Using the eigenstate clustering approach described in Appendix H, we identified (within one sector) the set of eigenstates that correspond to a hidden quantum number sector without identifying their corresponding integrals of motion. This procedure involved looking at the set of energy eigenstates coupled through the perturbation δĤ . After accounting for the hidden quantum numbers, our originally 8192-dimensional sector was reduced to four 2048-dimensional hidden sectors, suggesting that for the 27-site lattice there are likely two missing binary integrals of motion that we were unable to directly identify. Using the energy levels in the hidden sectors, we computed the average level-spacing ratio and observed that it indeed approached the GOE value as increased, in fact doing so much more rapidly than for the 18-site lattice (see Fig. 9). This seems to suggest that the non-GOE behavior for low could be a finite-size effect. We also numerically determined the ground state degeneracy of 10 random realizations of the Hamiltonians of Eqs. (52) and (53) with disorder strengths = 0.25, 0.5, 0.75, 1 on an 18-site kagome lattice. For all random Hamiltonians considered, the ground states were exactly four-fold degenerate, suggesting that these perturbed models are still Z 2 spin liquids.

V. DISCUSSION AND CONCLUSIONS
We have introduced a new approach, the symmetric Hamiltonian construction (SHC), for constructing Hamiltonians with desired symmetries. This method is general and can be used to construct families of Hamiltonians that commute (or anti-commute) with desired integrals of motion, such as zero mode and Wilson loop operators, and are invariant (or anti-invariant) under discrete symmetry transformations, such as point-group symmetries. In this work, we applied the SHC approach to design new topologically ordered Hamiltonians by providing as input to the method integrals of motion with topological properties and spatial symmetries. We analytically determined large families of superconducting Hamiltonians with Majorana zero modes (MZMs) and numerically found new Z 2 spin liquid Hamiltonians on the square and kagome lattices.
Using the SHC, we developed a general framework for designing Hamiltonians that commute with a pair of zero mode operators. In this framework, we provide as input the spatial distribution of two zero modes and as output produce Hamiltonians that commute exactly with those two zero modes and no others. These Hamiltonians can be put onto arbitrary lattice geometries, e.g., square, kagome, tetrahedral, or quasicrystal lattices, or even arbitrary graphs. Some examples of Hamiltonians that we designed with this framework are: a one-dimensional swave superconducting Hamiltonian that commutes with two exponentially-localized MZMs, a two-dimensional pwave superconducting Hamiltonian that commutes with Gaussian-localized MZMs, and two-dimensional Hamiltonians that commute with exotic semi-localized zero modes.
Using the SHC numerically, we discovered new classes of Z 2 spin liquid Hamiltonians on the square and kagome lattices whose properties differ from known solvable models. The Hamiltonians are not sums of commuting projectors nor are they frustration-free. Generically, these Hamiltonians possess many integrals of motion, though not enough to fully diagonalize the Hamiltonian. We observed that particular perturbations of these Hamiltonians, which break spatial symmetries but preserve the integrals of motion, exhibit GOE level-spacing statistics in their quantum number sectors. Many of the Hamiltonians that we found, even with these perturbations, possess numerically exact four-fold ground state degeneracy in finite-size systems.
There are many future directions for using the SHC to study topological order. Our general framework for designing zero mode Hamiltonians opens the door to the design of new experiments to search for MZMs. Using our framework, it is now possible for experimentalists to design a Hamiltonian with MZMs that best fits their experimental constraints, rather than focusing on a particular idealized model such as the Kitaev chain. To realize these MZM Hamiltonians in the lab, these experiments need some degree of control over the spatial distribution of some combination of chemical potentials, magnetic fields, superconducting pairings, or hoppings, in fermionic superconducting systems. Such control can potentially be realized in existing experiments, such as in Josephson junction arrays [48], superconductor-semiconductor heterostructures [49], twisted bilayer graphene [50], and strain-modulated superconductors [51,52]. Our framework also potentially allow theorists to design new model Hamiltonians with new exotic zero mode physics, such as zero modes that realize non-Ising non-Abelian anyons or zero modes displaying the physics of higher-order topological insulators and superconductors [53,54]. One can also use similar techniques as we used to discover Z 2 spin liquids in order to find other topologically ordered Hamiltonians. For example, by providing as input different topological symmetry operators, such as different types of Wilson loops, one can attempt to discover new model Hamiltonians with double semion or Fibonacci anyons.
Broadly speaking, the SHC is a tool that can be used to systematically study Hamiltonians with symmetries of interest. For example, it can be used to find all local Hamiltonians consistent with particular crystallographic symmetries. Once identified, these models can then be studied numerically or analytically to better understand their behavior. The SHC method can also be used to generate (potentially interacting) realizations of Hamiltonians from well-known symmetry classifications, such as the "ten-fold way" classification of non-interacting topological insulators and superconductors [55]. It could also be used to engineer Hamiltonians with properties that make them easier to study. Certain classes of interacting fermionic Hamiltonians with specific symmetries are known to be sign-problem-free, a property that allows for their efficient numerical simulation using methods such as quantum Monte Carlo [56][57][58]. The SHC method could be used to generate particular realizations of Hamiltonians with such symmetries, and thereby provide many new sign-problem-free Hamiltonians for numerical study. The availability of new numerically solvable models could provide valuable insights into strongly correlated systems.

VI. ACKNOWLEDGEMENTS
We acknowledge useful discussions with David Huse, Barry Bradlyn, Roger Mong, Curt Von Keyserlingk, David Pekker, Dmitrii Kochkov, Ryan Levy, Jahan Claes, and Di Luo. We also acknowledge useful conversations with Shivaji Sondhi, including a suggestion to consider our spin-liquid's level statistics. This work was supported by DOE de-sc0020165. BV acknowledges support from the Google AI Quantum team.

Appendix A: Derivation of commutant matrix
In this section, we derive Eq. (7), which expresses the commutator norm between HamiltonianĤ and integral of motionÔ in terms of the commutant matrix.
Consider a space of Hermitian operators spanned by a d-dimensional basis of Hermitian operatorsŜ a for a = 1, . . . , d. The commutator of two basis vectors in this space satisfies [Ŝ a ,Ŝ b ] = −[Ŝ b ,Ŝ a ] = −[Ŝ a ,Ŝ b ] † and is therefore anti-Hermitian. An anti-Hermitian operator can be represented as an imaginary number times a Hermitian operator. So if this space of operators is large enough, then the operator [Ŝ a ,Ŝ b ] can be represented as where the expansion coefficients f c ab are purely imaginary. Note that if this equation holds for all a, b in the space of operators we are considering, then the space is a Lie algebra with the commutator as the Lie bracket and f c ab as the (basis-dependent) structure constants of the Lie algebra.
Consider two operatorsĤ andÔ expressed in our basis so thatĤ = d a=1 J aŜa andÔ = d b=1 g bŜb with real J a , g b . Their commutator is where CÔ ≡ LÔ † OLÔ is the commutant matrix. The matrix O cc = tr Ŝ † cŜc /tr Î = Ŝ c ,Ŝ c is an overlap (or Gram) matrix, which is Hermitian and positive semi-definite, and is calculated by taking inner products between all operator strings. Since we require that thê S a form a basis, they are linearly independent, making the overlap matrix O non-singular and therefore positive definite. Since O is Hermitian and positive definite, LÔ † OLÔ is Hermitian and positive semi-definite. For orthonormal bases, O cc = δ cc and the commutant matrix simplifies to CÔ = LÔ † LÔ.

Appendix B: Properties of Pauli string basis
In this section, we describe some properties of the Pauli string basis, a complete basis for the space of Hermitian operators that can act on n qubits, and discuss how to compute the structure constants f c ab for this basis. The Pauli string operators, which are tensor products of the identity matrix and Pauli matrices, are defined in Eq. (2). These operators are Hermitian and unitary. These properties imply that the operators square to identity and have eigenvalues ±1. Moreover, except for the identity operatorÎ, the operators in this basis are traceless, implying that they have an equal number of +1 and −1 eigenvalues. Therefore, the Pauli stringsŜ a are highly degenerate, with two 2 n−1 -dimensional degenerate subspaces.
The term with the identity operator, which is non-zero only if (s 1 , . . . , s n ) = (t 1 , . . . , t n ), i.e., only ifŜ a =Ŝ b , is the only operator with non-zero trace in the expansion of the product. Therefore Ŝ a ,Ŝ b = tr Ŝ aŜb /tr Î = δ ab and the basis is orthonormal.
Next, we describe in more detail how to compute the product, and thereby the structure constants, of Pauli strings. The commutator of the two strings in Eq. (B1) is To compute the structure constants f c ab in the final line for a particular (a, b)-pair in our basis, we need to apply Eq. (B2). To precisely illustrate this calculation, we need to introduce some notation. We say that two Pauli strings "agree" on p sites i 1 , . . . , i p when the operators on those sites match: s i1 = t i1 , . . . , s ip = t ip . They "trivially disagree" on q sites j 1 , . . . , j q when the operators on those sites do not match, but one of the two operators on the site is identity: s j1 = t j1 , . . . , s jq = t jq and s jm = 0 or t jm = 0 for all m = 1, . . . , q. Finally, the two strings "non-trivially disagree" on r = n − p − q sites k 1 , . . . , k r when the operators do not match and are both non-identity operators: s k1 = t k1 , . . . , s kr = t kr and s k1 , t k1 , . . . , s kr , t kr ∈ {1, 2, 3}. According to Eq. (B2), the sites that agree become identity operators, while the sites that disagree, both trivially and non-trivially, become Pauli matrices. We then see that From this result, we see that if r is even, then the two terms are purely real and cancel, leading to [Ŝ a ,Ŝ b ] = 0 and f c ab = 0 for all c. If r is odd, then the two terms are purely imaginary and add, resulting in a single c for which f c ab = 0. Therefore, the structure constant associated withŜ a andŜ b is The structure constants in Eq. (B4) can be computed efficiently, in time O(q + r), for each pair of operator stringsŜ a andŜ b . Therefore, for a d-dimensional basis of k-local Pauli strings, all of the structure constants for the basis can be computed in time O(kd 2 ).

Appendix C: Properties of fermion string basis
In this section, we describe some properties of the fermion string basis, a complete basis for the space of Hermitian operators that can act on n fermions. A similarly defined basis was considered in Ref. 59.
Fermion strings, expressed in terms of the standard fermionic creation and anhillation operatorsĉ † i andĉ i , are defined in Eq. (3). Unlike the Pauli strings, fermion strings are neither unitary, traceless, nor orthonormal according to the Frobenius inner product.
As shown in Eq. (3), we classify the fermion string operators into three types.
Altogether, from our counting, we see that the {Ŝ a } fermion string basis consists of 2 n + 2 × 2 n−1 (2 n − 1) = 4 n linearly independent Hermitian operators and therefore spans the entire space of Hermitian operators.
Finally, we note that the product of two fermion strings,Ŝ aŜb , and therefore the commutator, is nontrivial to calculate. Rather than working out the commutator and structure constants in the fermion string basis, we compute the structure constants in the Majorana string basis, which spans the same space of operators. We can then convert to and from the fermion string basis as needed by applying an invertible basis transformation, as discussed in the next section.

Appendix D: Properties of Majorana string basis
In this section, we describe some properties of the Majorana string basis -another complete basis for the space of Hermitian operators that can act on n fermions. We also discuss how to compute the structure constants f c ab for this basis and how to convert from the Majorana string basis to the fermion string basis. The Majorana string basis, while more difficult to interpret physically than the fermion string basis, is more amenable to the computation of structure constants, making it useful for our methods.
and so the properties of Majorana strings can be thought of as being inherited from the Pauli strings.
To understand the properties of the Majorana string basis, it is important to understand the algebraic properties of theâ i ,b i ,d i operators. Operators with different labels i = j, satisfy the following commutation and anticommutation relations: which one can derive from the canonical fermionic anticommutation relations {ĉ i ,ĉ † j } = δ ij and {ĉ i ,ĉ j } = 0. One can also show that operators with identical labels i = j obey the relations: Next, we clarify the imaginary prefactor in Eq. (4). Note that, unlike the Pauli string basis in Eq. (2), the Majorana string basis in Eq. (4) involves products, not tensor products, of operators. This is because the ordering of the operators are important in our calculations due to the anti-commutation relations of theâ i ,b i operators shown in Eq. (D1). We now illustrate the effect of the operator ordering with an example. Suppose that, of the n sites of a Majorana stringŜ a , m AB of the sites haveâ i or b i operators on them so that m AB = n i=1 (δ ti,1 + δ ti,2 ). For example, the Majorana stringŜ a = i σaâ 2â3b4â5 has m AB = 4. For thisŜ a to be Hermitian, we require thatŜ † a = (−i) σaâ 5b4â3â2 =Ŝ a upon the reordering of the anti-commutingâ i ,b i operators. In general, reversing the order of these operators can be done with m AB 2 = m AB (m AB − 1)/2 swaps, which multiplies the operator by a sign (−1) m AB (m AB −1)/2 . Therefore, we use the convention that σ a = m AB (m AB − 1)/2 mod 2 to makeŜ a Hermitian. Note, from Eq. (D1), that the identity operator and parity operatorsd i commute with operators on different sites, so they do not contribute signs like theâ i andb i operators do. Now, we discuss how to compute a product of Majorana strings, which will demonstrate the orthonormality of the Majorana string basis. Consider a pair of length n Majorana stringsŜ a = i σaτ s1 1 · · ·τ sn n S b = i σ bτ t1 1 · · ·τ tn n where the σ a , σ b ∈ {0, 1} andτ ti i are as defined above and s i , t i ∈ {0, 1, 2, 3}. The product of these two operators iŝ where s ab = ±1 is a sign picked up from reordering thê a i andb i operators. (In practice, we compute the sign s ab by sorting theâ i andb i operators in theŜ aŜb string with a stable sorting algorithm and counting the number of swaps performed. An odd number of swaps leads to a minus sign.) Sinceâ i ,b i ,d i , and strings of these operators are traceless, we see that the only possible term with non-zero trace in the final line is the identity operator, which occurs whenŜ a =Ŝ b . Therefore, just like the Pauli string basis, the Majorana string basis is orthonormal, satisfying Ŝ a ,Ŝ b = tr Ŝ aŜb /tr Î = δ ab . Next, we discuss how to compute the structure constants. From Eq. (D4), we see that the commutator of two Majorana strings is After the reordering of the Majorana operators, the calculation of f c ab andŜ c parallels the one for the Pauli string basis. To compute the structure constants f c ab in the final line of Eq. (D5) for a particular (a, b)-pair in our basis, we need to apply Eq. (D3). To describe this calculation, we use the same notation as we used for the Pauli strings. We say that two Majoranas strings "agree" on p sites i 1 , . . . , i p when the operators on those sites match: s i1 = t i1 , . . . , s ip = t ip . They "trivially disagree" on q sites j 1 , . . . , j q when the operators on those sites do not match, but one of the two operators on the site is identity: s j1 = t j1 , . . . , s jq = t jq and s jm = 0 or t jm = 0 for all m = 1, . . . , q. Finally, the two strings "non-trivially disagree" on r = n − p − q sites k 1 , . . . , k r when the operators do not match and are both non-identity operators: s k1 = t k1 , . . . , s kr = t kr and s k1 , t k1 , . . . , s kr , t kr ∈ {1, 2, 3}. According to Eq. (D3), the sites that agree become identity operators, while the sites that disagree, both trivially and non-trivially, be-comeâ i ,b i ord i operators. We then see that where the notation · indicates that the bracketedτ ti Finally, we discuss the conversion of Majorana strings to linear combinations of fermion strings. This conversion is done by applying the definitionsâ Inserting these relations into Eq. (4) and expanding, we see that we obtain 2 k terms made of products ofĉ i andĉ † i , where k is the number of non-identity terms in the Majorana string. These terms can cancel and can be combined to form non-diagonal fermion strings, which involve Hermitian conjugates. To correctly convert these terms to the fermion strings of Eq. (3), we need to normal order our expanded operators and reorder them so that they follow our label ordering convention. We implement this process algorithmically to build up a basis transformation matrix B ab , which is invertible but not unitary, since the fermion string operators are not orthonormal. The construction of the B matrix takes time O(d min(2 k , d)).

Appendix E: Representations in the operator string basis
Here we discuss how to compute the representations D g of a symmetry transformation g ∈ G in the operator string basis for a few common discrete symmetries.
For fermionic systems, spatial unitary symmetry transformations can be represented by their action on the fermionic creation and anhillation operatorsĉ † where i, j label lattice site degrees of freedom and U ji is a unitary matrix. For example, for reflection and rotation symmetries, the matrix U ji is a permutation matrix that specifies how lattice site labels are permuted under the transformation. Eq. (E1) and its generalizations provide us with a rule for how to transform fermionic operator stringsŜ a . For Majorana fermions, Eq. (E1) can be reexpressed aŝ Now we can see that, for a Majorana string operator made of many Majorana fermion operators, these rules specify how the string transforms. For example, the string iâ ibj transforms aŝ For a space group symmetry, this is particularly simple and the U ki , U lj , U (kl),(ij) matrices are all permutation matrices. For non-spatial symmetry transformations, such as charge-conjugation or time-reversal symmetry, the transformations also involve changes in sign in addition to permutations (see Tab. A1). For Pauli strings, spatial symmetry transformations work the same way as for Majorana strings: they simply permute the labels of the Pauli matrices. The time-reversal symmetry transformation, on the other hand, involves an additional sign for every Pauli matrix, Here we derive a large family of Hamiltonians that commute with desired zero modes. First, we look for two-site Hamiltonians that commute with a single zero mode. Then, we look for two-site Hamiltonians that commute with a pair of zero modes. Finally, we discuss how these two-site Hamiltonians can be used to construct manybody Hamiltonians with desired zero modes.
Suppose that we wish to find Hamiltonians that commute with a single zero mode of the formγ = j (α jâj + β jbj ). On the two sites i and j, this zero mode has supportγ ij = α iâi + β ibi + α jâj + β jbj . To find two-site fermion-parity-conserving Hamiltonians, we construct the commutant matrix Cγ ij in the 7-dimensional basis spanned by the Majorana stringŝ In this basis, the 7 × 7 commutant matrix is where the lower triangle of this matrix is specified by the upper triangle since it is symmetric. We find that this matrix has three null vectors and four degenerate eigenvectors with positive eigenvalue 4(α 2 For α i , β i , α j , β j = 0, the three null vectors correspond to the following three operators that commute withγ ij j + H.c.. In summary, when we looked for two-site Hamiltonians that commute with one zero mode, we found three such Hamiltonians. Now suppose that we wish to find Hamiltonians with two particular zero modes,γ (1) = j (α jb j ) andγ (2) = j (α (2) jâ j + β (2) jb j ), with supportγ (1) ij and γ (2) ij on sites i, j. To find the two-site Hamiltonians that commute with both of these zero modes, we examined the null space of the sum of their commutant matrices Cγ(1) ij + Cγ (2) ij in the same 7-dimensional basis as before.
We found a one-dimensional null space of this matrix, which corresponds to a unique two-site Hamiltonian that commutes with both of the desired zero modes. This unique Hamiltonian, converted into complex fermions, iŝ Note that if we consider the zero mode coefficients as vectors γ (1) ≡ (α j ) T and γ (2) ≡ (α (2) i , β (2) i , α (2) j , β (2) j ) T , then each of the parameters of the Hamiltonian are indefinite quadratic forms, e.g.,t R ij = γ (1) T Aγ (2) for a particular antisymmetric matrix A.
Next, we checked if the operatorĥ ij commuted with zero modes other thanγ (1) andγ (2) . We performed this check by computing the commutant matrix Cĥ ij in the four-dimensional basis spanned by the Majorana fermionsâ i ,b i ,â j ,b j . For the bond operator specified by Eqs. (G2) and (G3), the commutant matrix has a null space that is exactly two-dimensional and spanned precisely by theγ (1) ij andγ (2) ij zero modes. The two remaining eigenstates of the 4 × 4 Cĥ ij matrix are degenerate with eigenvalue 4∆ε ij , where This eigenvalue gap is non-negative (∆ε ij ≥ 0) by the Cauchy-Schwarz inequality and is positive as long aŝ γ (1) ij ∝γ (2) ij . It is largest when the two zero modes are orthogonal on sites i and j.
Let us examine what happens when we attempt to construct an N -site Hamiltonian made of only the bond op-eratorsĤ In this procedure, imagine that we have beforehand decided on a desired pair of zero modes,γ (1) ,γ (2) , so that we have specified the α k , β (2) k parameters for all k = 1, . . . , N . This in turn specifies all of theĥ ij operators. To avoid the case where ∆ε ij = 0, we also require thatγ (1) ij ∝γ (2) ij for the (i, j) pairs that we consider below.
Suppose that we start building our Hamiltonian from the zero operator, so that J ij = 0 for all i, j. In this case, there are 2N zero modes that (trivially) commute witĥ H ZM = 0:â k ,b k for k = 1, . . . , N . Note that each site has two zero modes and each bond has four. Next, suppose that we turn on a bond J lm = 0 for a particular pair of sites l, m. On sites l, m, the bond operator gaps out two of the four zero modes, so that only two zero modeŝ γ (1) lm ,γ (2) lm on sites l, m commute with the Hamiltonian. After laying down the first bond, there are 2N − 2 zero modes. Moreover, the zero modes on sites l, m are now constrained to locally matchγ (1) andγ (2) . Now, consider iterating the procedure and laying down one bond at a time. If we connect the bonds together, e.g., so that J lm , J mn , J np = 0, then each bond we lay down eliminates two zero modes from the system and acts as a constraint that forces the zero modes on those sites to match our desired zero modes. If we think of the non-zero J ij as being the edges of a graph, then we can see that we are building a connected component into the graph and that the only zero modes that commute with the Hamiltonian on that connected component are constrained to matcĥ γ (1) andγ (2) . In general, if we build theĤ ZM Hamiltonian to have N C connected components, then there will be 2N C zero modes, each of which are "pieces" of thê γ (1) andγ (2) zero modes. If N C = 1, i.e., the graph is connected, then the only two zero modes are exactly the entireγ (1) andγ (2) operators.
Appendix H: Level-spacing statistics of perturbed spin Hamiltonians The Hamiltonians discussed in Sec. IV have many integrals of motion. Using the eigenstates of these integrals of motion, we block diagonalized the Hamiltonians according to their quantum number sectors and analyzed the level-spacing statistics in particular sectors. Many of these Hamiltonians have significant eigenstate degeneracy in these sectors, which makes analysis of the level-spacing statistics unreliable or inconclusive. We observed that particular perturbations of these Hamiltonians possess the same integrals of motions as the unperturbed Hamiltonians, though at the cost of breaking spatial symmetries. These perturbed Hamiltonians, however, have almost no degenerate eigenstates, which allows us to gather more reliable level-spacing statistics. Below, we describe how we computed the level-spacing statistics of these perturbed Hamiltonians while accounting for as many integrals of motion as possible.
On the square lattice, we analyzed the perturbed Hamiltonian of Eq. (43). Like the unperturbed model, the perturbed model commutes with straight-line Z Wilson loops and products of two straight-line X loops, which are listed in Eq. (40).
We diagonalized these Hamiltonians in a basis of states that are eigenstates of the known integrals of motion.
Here we describe how we found these eigenstates. The +1 eigenstates of theẐ L integrals of motion, for example, are simply product state spin configurations |S that satisfyẐ L |S = +|S , i.e., spin configurations that have an even number of down spins on the sites of loop L. The +1 eigenstates |S of theX LXL operators can be constructed from spin configurations |S in the following way: |S = 1 √ 2 (Î +X LXL )|S . Using these two facts, we are able to construct all +1 eigenstates of theẐ L and X LXL integrals of motion (with slight modification, we can build −1 eigenstates as well). In addition to the Z and X loop integrals of motion, there is also a global integral of motionB, slightly modified from the operator in Eq. (42), which is a sum of the terms of Hamiltonian Eq. (43) that lie on B-sublattice squares. We determine the eigenstates ofB by perturbing the Hamiltonian bŷ B. In particular, instead of diagonalizingĤ square + δĤ, we diagonalizeĤ square + δĤ + λB for a large constant λ. TheB-perturbation separates out the eigenstates of B so that they are far away from one another in energy. This allows us to numerically identify eigenstates with the same eigenvalue ofB.
We also attempted to account for integrals of motion that we were not able to identify directly. Such unknown integrals of motion, if left unaccounted for, would give rise to Poisson level-spacing statistics. This would occur when neighboring energy levels E n , E n+1 are in different quantum number sectors of these integrals of motion. Note that here we are considering a family of perturbed Hamiltonians,Ĥ 0 + δĤ, where [Ĥ 0 , δĤ] = 0. If there is a hidden integral of motionÔ that commutes withĤ 0 + δĤ for all , then [Ĥ 0 ,Ô] = [δĤ,Ô] = 0. This means that one can block diagonalize δĤ according to the eigenstates ofÔ and that eigenstates |ψ 1 , |ψ 2 ofÔ with different eigenvalues satisfy ψ 2 |δĤ|ψ 1 = 0. Moreover, energy eigenstates |ψ j ofĤ 0 + δĤ that are non-degenerate are also eigenstates ofÔ. Using these observations, we computed the perturbation overlaps δH ij = ψ i |δĤ|ψ j in the basis of {|ψ j } obtained with ED. We then reordered the states so as to block diagonalize the δH ij matrix. When block diagonalizing δH ij , we considered entries of the matrix smaller than 10 −6 to be zero. The eigenstates within the same block were considered as a "sector" of the hidden integrals of motion. We performed our level-spacing ratio calculations using the energies in these sectors, which are also contained within the known integral of motion quantum number sectors mentioned above.
Accounting for both known and unknown integrals of motion as described above, we computed the levelspacing ratios of 4096 eigenvalues in the +1 sectors of the Z loops and X product loops for the Hamiltonians of Eq. (43) on an 8 × 4 square lattice. We did this for 10 random realizations of the perturbed Hamiltonians for from 0.05 to 2. After accounting for the global integral of motion, the energy eigenstates typically separated into sectors of 64 states. After accounting for the hidden integrals of motion, these 64 states separated further into three sectors containing 28 states, 35 states, and 1 state. We computed the level-spacing ratios between states of neighboring energies in the 28 and 35-dimensional sectors and averaged the results over all such sectors and over random realizations of the Hamiltonians. The resulting average level-spacing ratios for from 0 to 6 are depicted in Fig. A1.
On the kagome lattice, we analyzed the perturbed Hamiltonians of Eqs. (52) and (53). Both of these Hamiltonians commute with Z Wilson loops and products of two X loops. However, the Hamiltonian of Eq. (52) also commutes with the local hexagon integrals of motionẐ , while Eq. (53) is a generic Hamiltonian that does not.
Using the same techniques described above to account for known and unknown integrals of motion, we performed ED on the Hamiltonians of Eq. (52)-(53) for 3 × 2 × 3 = 18 site and 3 × 3 × 3 = 27 site kagome lattices. For the Eq. (52) Hamiltonians on the 18-site lattice, we performed ED in a 256-dimensional basis corresponding to the +1 quantum number sector of each of the Wilson loops and hexagon local integrals of motion. We were not able to identify any "hidden" quantum number sectors in this 256-dimensional space and found that all energies in this sector were unique. For the Eq. (52) Hamiltonians on the 27-site lattice, we performed ED in a 8192-dimensional basis. In this case, we did find "hidden" integral of motion sectors that were 2048-dimensional and contained no degeneracies. For the Eq. (53) Hamiltonians on the 18-site lattice, we performed ED in a 1024-dimensional basis. In this case, the global integral of motion (a modified version of op-eratorD from Eq. (51) that only contains the hexagon interactions of the Hamiltonian) split the space into three sectors, which are 256, 256, and 512-dimensional. Interestingly, while the first two sectors have no degeneracies, the third sector is doubly-degenerate in each of its energy eigenstates. In this case, we ignored the degeneracy when computing the level-spacing statistics. We did not find any hidden integral of motion sectors for these Hamiltonians. For the Hamiltonians of Eq. (52)-(53) at a particular disorder strength , we averaged the level-spacing ratios r of the sectors over many random realizations of the random variables h (100 realizations for the 18-site kagome lattice and 10 realizations for the 27-site lattice). The average level-spacing ratios as a function of are shown in Fig. 9.