Protection of parity-time symmetry in topological many-body systems: non-Hermitian toric code and fracton models

In the study of $\mathcal{P}\mathcal{T}$-symmetric quantum systems with non-Hermitian perturbations, one of the most important questions is whether eigenvalues stay real or whether $\mathcal{P}\mathcal{T}$-symmetry is spontaneously broken when eigenvalues meet. A particularly interesting set of eigenstates is provided by the degenerate ground-state subspace of systems with topological order. In this paper, we present simple criteria that guarantee the protection of $\mathcal{P}\mathcal{T}$-symmetry and, thus, the reality of the eigenvalues in topological many-body systems. We formulate these criteria in both geometric and algebraic form, and demonstrate them using the toric code and several different fracton models as examples. Our analysis reveals that $\mathcal{P}\mathcal{T}$-symmetry is robust against a remarkably large class of non-Hermitian perturbations in these models; this is particularly striking in the case of fracton models due to the exponentially large number of degenerate states.

Among these models, a particularly important and commonly studied class of non-Hermitian Hamiltonians is provided by PT -symmetric Hamiltonians which are invariant under a combination of parity and time-reversal.Despite being non-Hermitian, these Hamiltonians can exhibit real spectra [1-3, 89, 90].Intuitively, this may be attributed to a balance of gain and loss between the system and its environment.Mathematically, the protection is related to the fact that PT symmetry implies that eigenvalues come in complex-conjugate pairs such that isolated real eigenvalues cannot become complex immediately.When they "meet" with another eigenvalue, they can either stay on the real axis or form complex-conjugate partners; when the latter happens, PT is said to be bro-ken.Therefore, the analysis of PT -symmetry breaking is particularly subtle in systems with (approximate) degeneracies.
For symmetry-imposed degeneracies, the reality of the eigenvalues can be simply protected by the symmetry itself and the fact that eigenvalues must come in complex conjugate partners.A priori, this is different for degeneracies related to intrinsic topological order [91,92]: for instance, the toric code model [93] has four ground states on a torus, that are guaranteed to be (exponentially) close in energy, even if all unitary symmetries are broken; similar statements apply to other spin-liquid phases.An even more dramatic ground-state degeneracy (GSD), that scales exponentially with linear system size, is realized in fracton models-novel quantum states of matter that are characterized by excitations with restricted mobility [94,95].Similar to spin-liquids, the GSD of fracton phases is topological in the sense that the different ground states are locally indistinguishable.One might be tempted to conclude that turning on a non-Hermitian, PT -symmetric perturbation in such systems will immediately lead to complex ground-state energies.Contrary to these expectations, we demonstrate in this paper that the reality of the ground-state eigenvalues in these phases can be surprisingly robust against a large class of such perturbations, even if all unitary symmetries are broken and in the presence of exponentially many degenerate states.
More specifically, we study under which conditions the eigenvalues of a given (almost degenerate) subspace of a Hermitian quantum system will stay real upon adiabatically turning on a non-Hermitian perturbation such that the total Hamiltonian commutes with a generalized PT symmetry.Here, "adiabatically" refers to keeping the gap to all other states finite and "generalized PT " indicates that P does not have to be spatial inversion, but might be any unitary operator.We first discuss a general mathematical condition for the eigenvalues to stay real and, hence, PT symmetry to be protected.We then demonstrate that this condition has strong implications for the protection of PT symmetry in the ground-state manifold of systems with topological GSDs, taking the toric code [93], the X-cube model [96], the checkerboard models [96,97], Haah's 17 CSS cubic codes [98], and the large class of quantum fractal liquids of Ref. [99] as examples.It is found that PT symmetry will be preserved on systems with even linear system sizes, L j , (in some Haah codes, divisibility by 4 is required) for a large class of perturbations, while it is generically fragile in systems with odd L j .
We emphasize that understanding the preservation or breaking of PT symmetry is not only one of the central theoretical questions of PT -symmetric quantum mechanics, but also of practical relevance for experimental realizations and potential applications of effectively non-Hermitian systems.We hope that our framework for predicting the stability of the reality of eigenvalues and the presence or absence of exceptional points will provide greater control over the effects of non-Hermitian perturbations, which is, e.g., important for the observation of power-law oscillations [18,[100][101][102] and the potential applications as topological lasers [103][104][105] and sensing devices [106,107].
The remainder of the paper is organized as follows.In Sec.I, we define the type of non-Hermitian Hamiltonians we are interested in, PT symmetry, and the more general condition of pseudo-Hermiticity.We also discuss the general, mathematical condition for colliding eigenvalues to stay real.It is first applied to the toric code, in Sec.II, to the X-cube, checkerboard models, and Haah's codes in Sec.III, and finally to the fractal liquids of Ref. [99] in Sec.IV.A summary of our findings is provided in Sec.V.

I. PSEUDO-HERMITIAN PERTURBATIONS
We start with a general explanation of the class of non-Hermitian perturbations under consideration.To this end, let us first assume that our non-Hermitian Hamiltonian H admits a complete biorthonormal eigenbasis {|ψ , |φ } [45], which means that (1) This is equivalent to the statement that H is diagonalizable, which is a very natural assumption for a generic (non-Hermitian) Hamiltonian of a physical system.Note, however, that it can be violated, most importantly at exceptional points [43,44], which we will discuss separately below.
In the study of non-Hermitian perturbations to quantum systems, it is common to further assume that these Hamiltonians are PT -symmetric [1-3, 89, 90], i.e., [H, PT ] = 0, where P can be abstractly defined as any unitary operator that squares to 1, and T is complex conjugation in a certain basis.Doing so imposes additional restrictions on the spectrum of H. Eigenvalues must come in complex conjugate pairs, as H(PT ) |ψ n = E * n (PT ) |ψ n .Importantly, this means that if one starts with a Hermitian, PT -symmetric Hamiltonian and applies a PT -symmetric non-Hermitian perturbation, isolated eigenvalues cannot become complex on their ownthey must merge with another eigenvalue on the real axis before becoming complex.This feature leads to the reality of energy spectra generally being robust to sufficiently small PT -symmetric perturbations, although degenerate subspaces are not necessarily protected from becoming complex.When PT |ψ n ∝ |ψ n , PT symmetry is said to be "unbroken" and the associated eigenvalues are real.Once eigenvalues meet and become complex, PT symmetry is "broken" and |ψ n is not an eigenstate of PT any more.
In this work, however, we do not restrict ourselves to PT symmetry, and instead impose a closely related but more general condition of pseudo-Hermiticity [108][109][110].A Hamiltonian H is pseudo-Hermitian if there exists a linear operator η, which we will refer to as the metric operator, such that In this paper, we take H to consist of a Hermitian component, H 0 , and a non-Hermitian perturbation, V , with magnitude that we control with ∈ R: For the Hermitian part H 0 , Eq. (2a) implies [η, H 0 ] = 0, i.e., η is a symmetry of the unperturbed Hamiltonian.We also take η to be unitary, so that Eq. (2a) is equivalent to ηHη † = H † .The purpose of this work is to derive and discuss general conditions under which (certain physically relevant parts of) the spectrum of H in Eq. (2b) can stay real upon adiabatically turning on .This condition of pseudo-Hermiticity (2a) is manifestly identical to PT symmetry with P ≡ η −1 provided H is symmetric, H = H T .In fact, it was shown [111] that any PT symmetric, finite-dimensional Hamiltonian is also pseudo-Hermitian.For this reason and since it does not involve any anti-linear operators and, thus, does not require a choice of basis, we focus on pseudo-Hermiticity in this work.Moreover, pseudo-Hermiticity gives a more systematic way of constructing non-Hermitian perturbations V to Hermitian models: one can immediately obtain all the possible choices of η as it has to be a symmetry of the unperturbed, Hermitian part, H 0 , of the model, which then specifies the suitable non-Hermitian perturbations.

A. Protection of reality of energies
If H is pseudo-Hermitian, complex eigenvalues also must come in conjugate pairs, since the combination of Eqs. ( 1) and (2a) implies Hη −1 |φ n = E * n η −1 |φ n .As is the case with PT -symmetric perturbations, this means that if a non-Hermitian perturbation is applied, the reality of isolated eigenvalues is stable to small pseudo-Hermitian perturbations.If a group of eigenvalues are degenerate (or almost degenerate) under H 0 -as is common in models involving symmetries or topological superselection sectors-they are generally not stable to pseudo-Hermitian perturbations.In these cases, we identify two main mechanisms by which these degenerate eigenvalues can stay real under pseudo-Hermitian perturbations: (I) The first method of ensuring degenerate eigenvalues stay real is simply to preserve the degeneracy under pseudo-Hermitian perturbations.
Pseudo-Hermiticity implies that if degenerate eigenstates are going to become complex, they must acquire imaginary parts with opposite signs.If one forces the (in general complex) eigenvalues to remain degenerate, this can never be satisfied for a non-zero imaginary component, unless the eigenvalues meet with another set of symmetry-unrelated eigenvalues.The latter, however, requires a sufficiently large value of , as symmetry-unrelated states are generically not degenerate for = 0.The symmetries enforcing the degeneracy can be unitary symmetries, fermionic timereversal symmetry [56], or even bosonic time-reversal symmetries unique to pseudo-Hermitian systems [112].
(II) The second mechanism is more subtle and our main focus in this work.If a pseudo-Hermitian term breaks all symmetries protecting the degeneracy, the eigenvalue splitting will generally be nonzero.This splitting can be either real or imaginary.However, one can show that if all the eigenstates of H 0 of the degenerate (or almost degenerate) subspace of interest have the same eigenvalue under η, then this splitting will always be real.This mathematical fact can be readily understood within the framework of G-Hamiltonian systems developed by Krein, Gel'fand and Lidskii [113,114] in the 1950s for the case of Hermitian η.In Appendix A, we provide a simple and physically insightful proof to all orders of perturbation theory that works for η being Hermitian or unitary.Furthermore, our analysis shows that, if the eigenvalues of η are identical, the projections of the associated eigenstates to the (almost) degenerate subspace of H 0 will be orthogonal to first order in and to zeroth order in the limit of a large gap to the rest of the spectrum; it also follows that, as long as the energetic separation of the subspace of interest from the rest of the spectrum is sufficiently large, they will be approximately orthogonal in the entire Hilbert space, even though the Hamiltonian is not Hermitian any more.This is very different when the eigenvalues of η are not the same.In that case, there can be exceptional points [43,44], where the Hamiltonian is defective, eigenstates coalesce and become identical, irrespective of how large the gap to the other states of the system is.
Intuitively, this is related to the fact that the Hamiltonian restricted to the degenerate subspace is Hermitian: denoting the degenerate eigenfunctions by |ψ i and writ- The complete argument of Appendix A involves constructing a full effective Hamiltonian for the degenerate subspace, and showing that it is Hermitian via similar reasoning.

B. Remarks on condition for reality
Some comments should be made regarding the possibility of multiple degeneracies and multiple metric operators.In the case of a two-fold degeneracy, i.e., two eigenvalues being identical, there are only two possibilities for the eigenvalues of η-either they both have the same eigenvalue under η, or they are different.In the former case, the splitting is always real.In the latter case, the energy splitting can be real or complex depending on the magnitudes of the matrix elements in the effective Hamiltonian.When there are more than two degenerate eigenstates, the full criteria becomes more complicated, as some eigenvalues may become complex while others stay real.For a concrete system, it should always be possible to determine the nature of the splitting through perturbation theory, using the methods described in Appendix A. However, we note that it is always the case that if all the unperturbed eigenstates have the same eigenvalue under η, their energies will stay real.
One can also consider a case where there is a two-fold degeneracy, but multiple possible choices of metric operators.If there are two metric operators, η 1 and η 2 , such that both eigenstates have the same eigenvalue under η 1 and different eigenvalues under η 2 , the degeneracy can be protected as a consequence of the mechanism (I) above: S = η −1 1 η 2 is, by construction, a symmetry of H and if S |ψ 1 = |ψ 2 , the eigenvalues will remain identical for = 0.However, the pseudo-Hermiticity of H with respect to η 2 can be broken without causing the eigenvalues to become complex.
In this paper, we focus on the protection mechanism (II) for the reality of the eigenvalues, i.e., on cases where pseudo-Hermitian perturbations break all the relevant symmetries, eigenvalue degeneracies are not preserved and hence the interplay between the metric operator and the unperturbed eigenstates are important in deducing whether the energies stay real.The general procedure for utilizing this phenomenon goes as follows.First, specify a subspace of interest, whose energies are separated from the rest of the spectrum.Next, identify the unitary symmetries under which the subspace has a definite eigenvalue under.These symmetries will yield a class of non-Hermitian perturbations-namely, those that are pseudo-Hermitian with the symmetry as a metric operator-for which the degenerate eigenvalues will stay real.
This notion of stability is useful in quantum systems when the subspace under consideration is well-separated The toric code is defined on a cubic lattice, with Pauli spins on each edge.The Hamiltonian is a sum of stabilizers, consisting of either a product of X operators on a plaquette (red), or a product of Z operators adjacent to a vertex (blue).We here show the case Lx = 3 and Ly = 2. (Bottom) On an even-by-even lattice (depicted here for a 2×2 lattice), every site can be covered by a combination of nonoverlapping plaquette operators (red).The four sites that are seemingly not covered by plaquette operators are redundant due to periodic boundary conditions.This covering is also possible with vertex operators.The coverings of larger lattices can be accomplished by sewing together copies of this 2 × 2 covering-of course, this only works for even-by-even lattices.
from the rest of the spectrum.In the remainder of the paper we will be concerned with gapped many-body systems with several degenerate ground states and discuss under which conditions the ground-state energies can remain real, provided the perturbations do not close the gap between the ground and excited states.

II. NON-HERMITIAN TORIC CODES
We begin with a study of non-Hermitian perturbations to the two-dimensional toric code [93], focusing on the reality of the ground-state subspace.Non-Hermitian generalizations of the toric code [81,83] or closely related models [82] have recently been studied; these works, however, have a different focus and a systematic understanding of the stability of PT symmetry or, more generally, of the reality of the spectrum in the ground-state subspace remains unexplored.
The toric code is defined on a square lattice, with Pauli spins on every edge, see Fig. 1.We denote the number of sites along the x and y directions by L x and L y and focus on periodic boundary conditions.The Hamiltonian is where we have introduced the vertex operators, A c , which cover the four spins adjacent to a vertex c, and the plaquette operators, B p , which cover the four spins on a plaquette p, Unless stated otherwise, we will use α = β = 1.In accordance with quantum code terminology, we refer to A c and B p collectively as "stabilizers."Each term in Eq. ( 4) commutes with the rest of the Hamiltonian, so the ground states can be obtained by minimizing the energy of each operator independently.Any state |ψ in the ground-state subspace satisfies A c |ψ = B p |ψ = |ψ .If defined on a torus, one can define loops of Z or X operators that wind around either of the two cycles of the torus.These logical string operators, that cannot be deformed to the identity by applications of stabilizers, imply a fourfold degenerate ground state, with string operators acting irreducibly within that subspace.

A. Pseudo-Hermitian perturbations
We are interested in pseudo-Hermitian perturbations to Eq. ( 4) and how they affect the degenerate ground states.To this end, let us first focus on three possible choices of η, where the product involves all sites of the system, and postpone the discussion of other options to Sec.II D below.One can easily check that H T C , η = 0.In contrast with many other features of the toric code, which only depend on the topology of the manifold, the eigenvalues of the ground states under η in Eq. ( 5) are highly sensitive to the system size.On an even-by-even lattice, the entire ground-state subspace has the same eigenvalue under η.This can most easily be seen by the fact that η can be written as a product of plaquette and vertex operators, and must give eigenvalue +1 in the ground-state subspace as a result, see Fig. 1.This cannot be accomplished on any other lattice size for the following reason: a straight line drawn along the x (y) direction going through the centers of the plaquettes will intersect L x (L y ) sites.If we attempt to cover the full lattice with plaquette operators, the placement of an additional operator will always change the number of covered sites on the line by an even amount.The same holds true for vertex operators and lines drawn through the vertices.Therefore, if either L x or L y is odd, the full lattice can never be assembled solely from stabilizers.The fact that η cannot be written as a product of stabilizers is sufficient to show that not all ground states can have the same eigenvalue under η.To see this, suppose that all ground states have the same eigenvalue under η.If this holds, then we can add η to the group of stabilizers of the toric code without modifying the GSD.If η is independent from the rest of the stabilizers, we arrive at a contradiction, since increasing the number of independent stabilizers lowers the GSD.
The observation that all the ground states have the same sign under η can also be seen by noting that η commutes with all the logical string operators on an evenby-even lattice, which take the system between different ground states.On an odd-by-even or an odd-by-odd lattice, η anti-commutes with at least one of the logical string operators, which in both cases lead to two ground states having eigenvalue +1 and the other two having eigenvalue −1.
What sort of perturbations, V , can we add to our Hamiltonian for which ηV η † = V † ?Writing V = iO, this requires which reduces to {η, O} = 0 for Hermitian O. Taking η to be the product of Y operators for concreteness, this means that O can be a sum, O = t g t O t , g t ∈ R, over terms O t which are products of Pauli matrices, only constrained to contain an odd number of X i and Z i .This includes a large class of perturbations such as random, planar fields, V = i i (g i1 X i + g i3 Z i ), g i1 , g i3 ∈ R, and highly non-local terms, such as i i<j<k g ijk X i X j X k , or Eq. ( 6) separately, there is no relation between the prefactors of the different terms required and we can think of them as random, non-Hermitian disorder, that in general breaks all symmetries of the system (other than PT ).
In combination with our results of Sec.I, this implies that on an even-by-even lattice, the ground-state subspace of the toric code remains real under the large class of pseudo-Hermitian perturbations that satisfy Eq. ( 6) with η given by Eq. ( 5).As the eigenvalues must stay real for small perturbations, they never exhibit any square root singularities [44] and exceptional points are avoided.This is verified by exact diagonalization (ED) of the toric code spectrum in Fig. 2(a,b), where it can be seen that the ground-state energies can only become complex when meeting with the excited states.As such, the PT symmetry of the ground-state manifold is protected by the gap to the excited states.
In contrast, on a lattice that is not even-by-even, the ground states generically become complex immediately upon applying the same non-Hermitian perturbations.This sensitivity of the ground state to the system size FIG. 2. Spectrum of the toric with non-Hermitian random field perturbation, i i giXi, where gi was initialized randomly according to a Gaussian distribution with mean and variance 1.In (a,b) and (c,d) we take the bare toric-code Hamiltonian ( 4) and the perturbed one, Eq. ( 7) with Gaussian distributed hi (mean 0 and standard deviation 0.4), as starting point, respectively.In (a,c), the real part of the energy is shown with red and gray referring to real-valued ground and excited energy levels, whereas eigenvalues with a complex part (broken PT symmetry) are indicated in blue.The corresponding imaginary parts can be found in (b,d) with red indicating the ground states, defined as those four states with the lowest Re(Ei).can be thought of as representative of the highly entangled nature of the toric code ground states.Even if one was to consider arbitrarily large system sizes, the toric code ground states are still able to "detect" whether the system size is even or odd.A similar interpretation of this phenomenon is that even for local perturbations, the order in perturbation theory in which the ground state energy splitting will occur necessarily involves a non-local operator which winds around the torus and, as such, can be sensitive to (the parity of) the system size.

B. Starting with perturbed toric code
One might wonder whether the remarkable protection of PT and reality of the ground-state energies is just a consequence of the highly fine-tuned and exactly solvable toric code Hamiltonian (4) or a more general property of the underlying topologically ordered phase.To investigate this, let us take as our base Hamiltonian the toric code with some small Hermitian perturbation, for example a field along the Z-direction with in general spatially varying amplitude, The perturbation in Eq. ( 7) forces us to choose η = i Z i in Eq. ( 5), as it is the only one that commutes with the Hermitian Hamiltonian.Note that, of course, a completely random Hermitian field will break all symmetries and no η is possible; we are, however, not interested in this case as the Hamiltonian would break PT explicitly and the question of whether it is broken spontaneously would become ill defined.
With the additional perturbation in Eq. ( 7), we no longer have exactly degenerate ground states for = 0, but a finite energy splitting that is exponentially suppressed by the system size.The four low-energy states will still all be even under η for an even-by-even lattice, since our perturbation respects the η symmetry.Consequently, the ground-state energies will stay real, even when they "meet" each other at finite , as long as the gap to excited states stays finite.The protection of PT symmetry is, thus, a more general property of the underlying phase with topological order.We also demonstrate this with a concrete example in Fig. 2(c,d).

C. Exceptional points
A surprising observation is that, while these Hermitian perturbations do not change whether the ground states become complex, they do change the nature of how they become complex.Non-Hermitian Hamiltonians can exhibit exceptional points [43,44].Here, the eigenvalues coalesce, the matrix becomes defective, i.e., also the eigenvectors become degenerate, and the eigenvalues exhibit a square-root singularity in the tuning parameter, in our case , in the sense that the difference of eigenvalues scales with √ 0 − .This has crucial consequences, e.g., for the Green's function that exhibits a pole of second order in addition to the conventional first-order pole [44].For pseudo-Hermitian or PT -symmetric Hamiltonians, exceptional points typically arise at the moment when two eigenvalues meet on the real line and become complex.If we start with an unperturbed toric code on an even-by-odd lattice and apply a non-Hermitian, pseudo-Hermitian perturbation V , such as an imaginary transverse field, the degenerate ground states can immediately become complex.However, this degeneracy is not an exceptional point, since the degeneracy occurs in the Hermitian limit and must admit a complete basis of eigenvectors.In contrast, if one first applies a Hermitian perturbation, such as in Eq. ( 7), and then V , we have verified by ED on a 2×3 lattice that the ground states will form an exceptional point when they meet each other on the real line to become complex, and the corresponding eigenstates become identical.We emphasize that this is true for arbitrarily small Hermitian perturbations.This is in stark contrast to systems with even L x , L y ; here eigenvalues must stay real for small perturbations, they never exhibit any square root singularities, and exceptional points are avoided.
To illustrate this subtle behavior of perturbed systems with odd system sizes, let us take a two-level system as an effective description of two ground states with opposite eigenvalue of η meeting to become complex.Denoting Pauli matrices acting in this subspace by σ x,y,z , we have η = σ z and the most general pseudo-Hermitian Hamiltonian has the form with the real-valued parameters E 0 , ∆, α, and ; the latter parameterizes the strength of anti-Hermitian perturbations as before.Note that the model is PT symmetric only if 2α/π ∈ Z.The right eigenvalues and eigenvectors of h in Eq. ( 8) are given by The eigenvalues meet when = ±∆ and become complex for | | > |∆|.When ∆ = 0, however, this is not an exceptional point as ψ ± → (±1, sign( )e iα ) T / √ 2, forming an orthonormal basis, and ∆E = E + − E − → 2| | scaling linearly with , for ∆ → 0. For ∆ = 0, instead, we get ψ + → ψ − when → ±∆, showing that the matrix becomes defective, and the difference of eigenvalues scales as ∆E ∼ 2 √ 2 0 √ 0 − , for near 0 = ±∆.It is also readily verified that the overlap, φ ± |ψ ± , with the corresponding left eigenvector is non-zero except for the exceptional points = ±∆ = 0, where it vanishes; this "self-orthogonality" rules out the construction of a biorthogonal basis as in Eq. (1).In summary, we should think of the special case of vanishing splitting, ∆ = 0 or of the unperturbed toric code, as a fine-tuned limit where two lines of exceptional points, = ±∆, meet and give rise to a non-defective Hamiltonian, as required by Hermiticity.
We finally point out that this behavior is also visible on an even-by-even lattice when taking into account the excited states: as can be seen in Fig. 2(b,d), the imaginary part of the excited states that become complex at infinitesimal scales linearly in , whereas the PT symmetry breaking at finite exhibits the aforementioned square-root singularity.

D. Other metric operators
So far, we have focused on the three different choices of η in Eq. ( 5), but there are in principle many more possibilities for the bare toric code model (4), as it possesses many other symmetries.Here, we argue that our choices of η are unique provided we assume our anti-Hermitian perturbations can be disordered and are not required to have a specific spatial structure.
As a starting point, one might use spatial symmetries-lattice translations T x,y , four-fold rotation C 4 , and inversion I and combinations thereof.
For instance, η = I with is clearly a symmetry, [H T C , I] = 0, and it is easy to see that all ground states have the same eigenvalue under it for any system size (the same holds for T x,y but not for C 4 ).However, it is not a natural choice for a generic system with spatially varying Hermitian or non-Hermitian perturbations, such as those discussed above.For example, for an imaginary field, , it would require g iµ = −g −iµ and, hence, fine-tuning between spatially distant sites.Not even a site-independent complex field is possible.
Having established that choosing an η which relates spatially distant sites requires fine tuning, we focus on η that commute with all stabilizers separately.This requirement can alternatively be thought of as a restriction to symmetries that are preserved in the presence of spatial disorder in the couplings of the bare toric code, i.e., α → α c > 0, β → β p > 0 in Eq. ( 4).This leads to two distinct classes of possible η, schematically given by where "logical strings" stands for strings of X i or Z i operators along a non-contractible loop of the torus connecting the different ground states [115].
Clearly, the ground states will have the same eigenvalues under η in Eq. (10a) and, thus, stay real.For even system sizes, η in Eq. ( 5) are of this form and, as we have seen above, indeed admit a large class of non-Hermitian perturbations.
This is different for η of the form of Eq. (10b): the ground states will have different eigenvalues under η and PT symmetry is in general fragile.However, since η in Eq. ( 5) can be written in the form (10a), it is clear that Eq. (10b) cannot be spatially homogeneous on an even-by-even lattice, but must be distinct on a noncontractible loop around the torus; the same must hold for the associated non-Hermitian perturbation, which requires, again, significant spatial fine-tuning.Let us illustrate this latter point using the concrete example of η = i X i j∈P Z j , where P is a non-contractible closed path through the centers of the plaquettes.In that case, an imaginary field, V = i i along non-contractible loops, suitable metric operators are of the form of Eq. (10a) for even-by-even lattices.As the ground states will always have eigenvalue +1 under any such η, the reality of their eigenvalues and, thus, PT symmetry are protected.

E. Arbitrary system sizes
So far, we have focused our attention on even-byeven lattices since the homogeneous metric operators in Eq. ( 5) can be written as a product of stabilizers, while this is not possible on even-by-odd or odd-by-odd lattices; nevertheless, if one naively applies the covering shown in Fig. 1 on these lattices, one can obtain a modified metric operator η, defined as the product of Pauli operators, X i , Y i , or Z i , on all sites except for a single line (in the evenby-odd case) or two lines (in the odd-by-odd case) that wind around the odd lengths of the torus.In other words, η in Eq. ( 5) is necessarily of the form of Eq. (10b) on a lattice with at least one of L x , L y odd.Based on our previous discussion, this implies that the reality of the ground-state eigenvalues and PT symmetry are generically fragile on even-by-odd and odd-by-odd lattices.
We finally mention, for completeness, one less general but potentially useful immediate consequence.As follows from using η as metric operator, any pseudo-Hermitian perturbation V with anti-Hermitian part, (V − V † )/2, that has support only in a subregion of the system that is contractible around the odd lengths of the torus, will leave the ground-state eigenvalues real.

III. NON-HERMITIAN FRACTON MODELS
Our analysis of the toric code carries over to many well-known fracton models in three dimensions.Fracton models [94][95][96][97][98][99][116][117][118][119][120] constitute a unique phase of matter, characterized by excitations with restricted mobility, either by being immobile or only mobile in certain directions.These systems are typically gapped and have GSDs exponential in linear system size.In this section, we analyze various models with fracton ordernamely, the X-cube model, checkerboard model, and Haah's codes-and show that, like the toric code, the full ground-state subspaces are stable against a large class of non-Hermitian perturbations provided the linear system sizes along all directions are even.Unless stated otherwise, we take η to be defined in the same way as in Eq. ( 5), i.e., as a product of X, Y , or Z operators over all qubits in the system; as motivated in Sec.II above in the context of the toric code, these η provide the largest class of allowed non-Hermitian perturbations by virtue of being spatially homogeneous.

A. X-cube model
The X-cube model [96] is defined on a cubic lattice, with qubits living on the edges of the lattice.It has a Hamiltonian composed of mutually commuting terms (11) where A c = j∈∂c X j is the product of X operators on the 12 edges of the cube labelled by c, and B i v is a vertex operator, composed of four Z operators at vertex v in the plane perpendicular to the i'th direction.On an even-byeven-by-even lattice, our η operators in Eq. ( 5) can be assembled from these terms, thereby showing that the entire ground-state subspace has eigenvalue +1 under η, see Fig. 3.An identical argument as in the toric code case implies that η cannot be assembled from stabilizers on a lattice with any odd length.In combination with the fact that it commutes with all stabilizers A c , B i v separately, it must be of the form of Eq. (10b) for odd system lengths, with "logical strings" here referring to the logical stringlike operators of the X-cube model.
By our analysis of the toric code, this immediately implies that the X-cube ground states on an even-byeven-by-even lattice stay real under the non-Hermitian perturbations permitted by η, which includes the application of imaginary transverse fields, non-local terms like the ones considered for the toric code, and many others.One can check that all other features of non-Hermitian toric code perturbations, such as their additional stability against real perturbations and the ability to add contractible perturbations on lattices with odd system sizes, also hold.However, these features are more striking for fracton models: instead of a four-dimensional code subspace being protected against these perturbations, fracton models have a GSD that grows exponentially with system size; for the X-cube model on a three-dimensional torus, the GSD obeys The reality of the code subspace in the presence of pseudo-Hermitian perturbations holds for the Xcube model defined on general three-dimensional manifolds [119], provided the full space can be covered by plaquette or star operators.This sensitivity to system size may be surprising, since the X-cube model exhibits foliated fracton order [119].This means that the length of any of the sides of the Xcube model can always be extended by attaching layers of toric code and applying a series of local unitary transformations.In Appendix B, we present a detailed study of how the metric operators η behave under foliations.The end result is that, while the ground states can be extended by this foliation procedure, the foliation acts non-trivially on η, meaning that the interplay between η and the X-cube ground states can change depending on the system size.FIG. 3.For the X-cube model defined on an even-by-evenby-even lattice (shown here for 2 × 2 × 2), the full lattice can be covered by non-overlapping plaquette operators (red).The lattice can also be covered by vertex operators.Although this is more difficult to visualize, it can be generated by covering each 2D layer in a definite plane by a toric code covering.The remaining sites on the edges connecting these 2D layers can then be covered by chains of vertex operators.

B. Checkerboard model
The checkerboard model [96] is another example of a system with fracton excitations.This model has spins defined on the vertices of a three-dimensional cubic lattice, as opposed to the edges.By separating the cubes of the lattice with alternating labels A and B, each forming a three-dimensional checkerboard lattice, and denoting the cubic operators i∈∂c Z i and i∈∂c X i as Z c and X c respectively, the checkerboard model is given by the Hamiltonian The geometry of the checkerboard model requires it to be defined on an even-by-even-by-even lattice if periodic boundary conditions are imposed, since otherwise one cannot uniformly partition the cubes into A and B labels.On even-by-even-by-even lattices, the entire lattice can be covered by non-overlapping stabilizers, and therefore the ground-state subspace is even under any η in Eq. ( 5).Again, this implies that the ground-state energies of the checkerboard model always remain real under non-Hermitian perturbations that are pseudo-Hermitian under η.Although only small system sizes are accessible via ED, we have checked these predictions numerically for the checkerboard model on a 2 × 2 × 2 lattice.A Majorana version of the checkerboard model has also been studied [97], which simply replaces the Pauli spins with Majorana fermions γ i , i.e., the model has one Ma-jorana fermion per site i of the cubic lattice.By defining i∈c γ i = γ c , the Hamiltonian of the Majorana checkerboard model is Because the entire lattice can be covered with γ c∈A , all ground states of the system are even under the operator which can be interpreted as the total fermion parity, η ∝ α (c † α c α − 1/2), when combining pairs of Majorana fermions into auxiliary complex fermions c α .Therefore, the ground states remain real under perturbations of the form i O, where each term in O contains an odd number of Majorana operators, i.e., changes the total occupation of auxiliary complex fermions by an odd amount.

C. Haah's codes
Finally, we consider Haah's 17 CSS cubic codes [98], all of which are defined on a cubic lattice with two qubits per site i.Each cube has two stabilizers: one is built up of tensor products of Z and 1 operators on each site i, such as Z i ⊗ 1 i or Z i ⊗ Z i ; the other one involves tensor products of X and identity operators, e.g., X i ⊗ 1 i .The exact form of the stabilizers differs from code to code, but all have a sub-extensive GSD.We defer a more detailed discussion of these codes to Appendix C-our conclusion is that, with the choice of η analogous to Eq. ( 5), the behavior of the code subspace under pseudo-Hermitian perturbations is sensitive not only to whether the system lengths are even or odd, but also whether the system lengths are divisible by 4.Moreover, since not all cubic codes are symmetric under rotations, this behavior is dependent on which directions are even or odd, and which are divisible by 4. This admits eight different classes of codes, based on the relation between their code subspace stability under pseudo-Hermitian perturbations and their system sizes.These classes range from cubic code 7, whose code subspace stays real on all system sizes other than odd-by-odd-by-odd, and cubic code 17, where the code subspace only stays real if L x , L y are divisible by 4. We refer to Appendix C for a complete characterization of this behavior.

IV. NON-HERMITIAN QUANTUM FRACTAL LIQUIDS
In this section, we will generalize the previous analysis to also include another class of fracton models dubbed "quantum fractal liquids" [99] and reformulate the criterion of stability against non-Hermitian perturbations using a polynomial representation of Pauli operators.In this way, we will recover the criterion of stability of the toric code in an algebraic way and show that the reality of eigenvalues of the exponentially large number of ground states of quantum fractal liquids is protected against a wide range of non-Hermitian terms in the Hamiltonian.

A. Polynomial representation of operators
To set up the notation, we will briefly introduce the polynomial representation of operators, a commonly used technique [121].To this end, consider a polynomial of three variables, x, y, z, f = j,k, ∈Z c jkl x j y k z , c jk = 0, 1, over F 2 , meaning that all coefficients are to be understood modulo 2. This allows to define a corresponding Pauli operator whose components lie on the vertices of a cubic lattice in three dimensions in the following way Here Z jk (X jk ) is the Z (X) operator acting at vertex (j, k, ).For example, a stabilizer of the checkerboard model, given by the product of Pauli matrices on the eight vertices of a cube, corresponds to the polynomial f = 1 + x + y + z + xy + yz + xz + xyz.On a finite lattice, periodic boundary conditions are specified by imposing x Lx = y Ly = z Lz = 1.We denote the dual of f , obtained by taking x → x −1 , and likewise for y and z, by f .Certain relations can be expressed more concisely with this polynomial representation.Translating an operator Z(f ) one lattice site along the x-direction is simply given by Z(xf ), and likewise for translations in the y and zdirection.Additionally, the polynomials defined over F 2 naturally encode the commutation relations of the Pauli operators.To see this, consider the commutation polynomial, defined as f ḡ for two polynomials f and g.Writing f ḡ as Quantum fractal liquids are defined on a cubic lattice with two spins on every vertex.The form of their stabilizers is given by [99] Z(α, β), X( β, ᾱ), and translations thereof, where the two arguments of Z and X denote operators on the two distinct spins per site.Different choices of polynomials f and g define different models.Clearly, all stabilizers commute, as follows from the associated commutation polynomial, α β + β ᾱ = 2αβ = 0.
For codes defined by stabilizers of this type, the logical operators take the form for integer i = 0, 1, . . ., L x − 1, where we define It is straightforward to verify that the operators in Eq. ( 16) commute with the stabilizers and constitute logical operators if There are various ways to satisfy Eq. ( 18): the "trivial" solution, that works for any set, L x , L y , L z , of system sizes, is f = g = 1.This corresponds to layers of toric code in the (ŷ,ẑ) plane, upon noting that the bond variables of the toric code, see Fig. 1, can be seen as two qubits per vertex.In this case, (Z,X) i and r (Z,X) i in Eq. ( 16) become Z-, X-type string operators in the ith layer along the ŷ and ẑ direction, respectively.Another way of satisfying Eq. ( 18) that works for arbitrary isotropic system sizes, . However, the largest class of possible polynomials f , g and, thus, possible models is allowed in the isotropic case with L = 2 n L , since Eq. ( 18) will hold as long as f (1) = g(1) = 1 [99].Here, we refer to the latter set of models as "quantum fractal liquids," which have been shown to exhibit exponential scaling of the GSD, obeying log 2 GSD(2L) = 2 log 2 GSD(L) [99].Note, however, that the absence of string-like logical operators and mobile quasiparticles further requires that f and g are not algebraically related, i.e., that there are no integers n 1 and n 2 such that f n1 = g n2 (neglecting periodic boundary conditions).An example of a model free of string-like logical operators is provided by f = 1 + x + x 2 and g = 1 + x + x 3 .

B. Pseudo-Hermitian perturbations
As before, we are interested in adding pseudo-Hermitian perturbations to this class of models that will leave the ground-state subspace real.We take η to be defined analogous to Eq. ( 14) or, in polynomial representation, Any η in Eq. ( 19) will commute with all stabilizers (15).This readily follows from the associated commutation polynomial upon noting that h = h is invariant under multiplication by any monomial, physically related to the translation invariance of η, and that the number of monomials in both f and g must be odd.The latter is a consequence of Eq. ( 18) and of the observation that the parity of the number of terms of a polynomial f over F 2 is the same as that of any of its powers, f n with n > 0.
Based on our discussion of Sec.I, we want to analyze under which conditions the ground-state subspace is even under these operators to guarantee that their eigenvalues stay real.Previously, we had verified this by attempting to assemble η via the stabilizers of the model.In the set of models introduced above, the polynomial representation makes it easier to instead verify whether η commutes with all the logical operators (16), which in turn implies that all ground states have the same eigenvalue of η [and that η is of the form of Eq. (10a) rather than Eq. ( 10b)].
The condition for η to commute with all the logical string operators, is given by for any of the three possible choices in Eq. ( 19).This simple expression arises from the fact that h = h and that the logical operators come in exactly the form of operators relevant to the commutation polynomial, so one can verify that η commutes with all the string operators with one equation.It would technically suffice for hg and hf to be only a function of y and z, since one is only concerned with the commutations of operators like Z(h) and X(x i f ), but not those with relative shift along the ŷ or ẑ directions.However, recalling that h is invariant under multiplication by any monomial, there is no way for g or f to conspire to cancel out only the terms independent of y and z in hg and hf without simply giving 0. Another important consequence of h being invariant under the multiplication by any monomial is that Eq. ( 20) is satisfied if and only if f and g contain an even number of monomials.As argued above, Eq. ( 18) implies that f , g and, therefore, also f n , g n must contain an odd number of terms.Taken together, Eq. ( 20) is obeyed and, thus, the reality of the eigenvalues of the ground states is protected against pseudo-Hermitian perturbation with η given in Eq. ( 19) if L y and L z are even.Note that the x-direction is distinguished from the other two directions in this criterion, a reflection of the fact that the stabilizers given by Eq. ( 15) also distinguishes the x-direction.
Let us illustrate this for the different special cases of f and g noted above.Taking f = g = 1 corresponds to L x uncoupled layers of toric code and the above statement implies that the toric code is protected if and only if the number of sites in each in-plane direction is even, reproducing the result of Sec.II.Our current formalism, however, captures many more cases.For instance, we immediately conclude that any model with f = x n f , g = x ng and L x = L y = L z = L is protected only for even L.As the two polynomials are algebraically related, this two-parameter family of models is characterized by string-like logical operators and has excitations mobile along the direction n g ŷ − n f ẑ [99].Finally, as already noted above, quantum fractal liquids with arbitrary f and g, only constrained by f (1) = g(1) = 1, are in general defined on lattices with an even number of sites and, as such, are always protected against pseudo-Hermitian perturbations with metric operator in Eq. ( 19).

V. SUMMARY AND CONCLUSIONS
In this work, we studied the behavior of the eigenvalues of quantum many-body Hamiltonians of the form of Eq. ( 2), i.e., starting from a Hermitian system, H 0 , we turn on a non-Hermitian perturbation, V , and demand that the entire Hamiltonian be pseudo-Hermitian.Using pseudo-Hermiticity rather than PT symmetry is related to the fact that the former is more general than the latter [111]; we note, however, that all of the explicit examples considered here are both PT symmetric and pseudo-Hermitian.We analyzed whether the energies, E i , of a given subspace of interest of H 0 will remain real as long as the gap to other states of the system is finite (PT symmetry protected) or whether they can move into the complex plane without closing the gap (PT symmetry fragile).While symmetries can enforce degeneracies (E i = E i ) and protect eigenvalues from becoming complex in conjunction with pseudo-Hermiticity (E i = E * i ), we discussed that this is also possible in the absence of symmetries: if the eigenvalues of the metric operator η are the same for all states in the subspace of interest, E i are guaranteed to stay real and PT symmetry is protected.
We demonstrated that this criterion can be readily applied to various paradigmatic many-body models with crucial implications.As a first example, we took the toric code model (4) as unperturbed Hamiltonian, H 0 .On a torus, it exhibits four degenerate ground states and one would generically expect them to become complex when turning on V .However, we have shown that η of the form given in Eq. ( 5) allows for a large class of non-Hermitian perturbations; these are shown to leave the ground-state energies real on an even-by-even lattice, even if all symmetries are broken.They can only become complex and PT can only be broken in the ground-state subspace, when the gap to the excited states closes.In fact, we have argued that any sufficiently generic non-Hermitian perturbation (see Sec. II D) in a system with both linear system sizes even (at least one of them odd) will only allow for η of the form of Eq. (10a) [of the form of Eq. (10b)] and the ground-state eigenvalues are protected (not protected) from becoming complex.This sensitivity to system size reflects the highly entangled nature of the toric-code ground states.
We came to the same conclusions for the ground-state manifolds of the X-cube (11), the spin (12) and Majorana (13) checkerboard models, and for the fractal liquids of Ref. [99].In these cases, the stability of PT symmetry is even more surprising due to the enormous GSD that grows exponentially with system size.For Haah's 17 codes, the stabilizers have a slightly more complicated form and the minimal requirement for stability differs from code to code, although we observe several groups of codes which all obey the same requirements.This classification of Haah's codes based on stability of PT symmetry approximately follows previous classifications based on entanglement renormalization [122].
On a more general level, our work illustrates that PT symmetry and the reality of energies can be protected in the degenerate ground-state manifold of correlated many-body systems with different forms of topological order-even in the absence of any symmetries and although exceptional points are generically expected to be ambundant [79].By virtue of being exact and simple, our framework can be readily applied to a large class of systems and provides a systematic method for constructing pseudo-Hermitian perturbations that ensures the reality of the resulting eigenvalues.This is not only relevant for experimental studies [18,[100][101][102] and potential applications [103][104][105][106][107], but might also help deepen our theoretical understanding of non-Hermitian systems hosting exotic phases of matter, e.g., by providing novel ways of classifying spin-liquid or fracton phases according to their sensitivity to such perturbations.
We are interested in the behavior of the eigenvalues of a subset of (orthonormal) eigenstates, {|ψ i , i = 1, . . ., n}, of H 0 , which can be arbitrarily close or identical in energy but are well separated from all other eigenvalues.We refer to the space spanned by {|ψ i , i = 1, . . ., n} as the almost degenerate subspace.
To analyze how their eigenvalues, E i ( ), i = 1, 2, . . ., n, evolve when turning on V in Eq. (A1), we define the projectors P and Q, to the almost degenerate subspace and its complement.We use that the exact eigenstates, |Ψ i ( ) , obeying with the effective Hamiltonian As follows from Eq. (A2), the eigenvalues E i ( ), i = 1, 2, . . ., n, can be obtained by diagonalizing the effective Hamiltonian H eff in the almost degenerate subspace.Of course, this is not straightforward to do as the effective Hamiltonian itself depends on these eigenvalues; however, the effective-Hamiltonian formulation is a good starting point to develop a perturbative expansion.

Expansion in
Since we view the non-Hermitian part V as a perturbation to H 0 in our discussion in the main text, it is very natural to expand in .Note that P H 0 Q = 0, so the second term in the effective Hamiltonian (A3) is O 2 , H eff (E) = P H 0 P P V P + 2 P V QG (E)QV P. (A5) Let us now assume that we can obtain E i ( ) via perturbative expansion in .To keep the notation compact, let us define the operator T N which performs a Taylor expansion on a function or operator up to and including order N , i.e., T N [f (x)] := First, note that Q commutes with η which implies that G (E) and, thus, also V QG (E)QV are pseudo- Hermitian if E ∈ R. Since this holds for a continuous set of values of , this property holds for each order in the Taylor expansion separately.As such, it also applies to iterative diagonalization of Eq. (A7) will always yield real eigenvalues.Taken together we have shown that E i ( ), i = 1, 2, . . ., n, stay real to any order in .
If the eigenstates are exactly degenerate for = 0, the leading non-zero contribution to the energy splitting will determine whether the eigenvalues of H stay real or become complex.In most cases, the first order corrections, given by diagonalizing P H P , break the degeneracy.Since P H P is clearly Hermitian, our result is simple if the first order energy splitting is non-zero.In fact, a mathematical proof to first order in perturbation theory has been provided in Ref. [124].However, topological degeneracies are often broken only at higher orders in perturbation theory, so a more general result is required.
If, however, the degeneracy is already broken for = 0, our current perturbative approach cannot be used to understand whether the eigenvalues stay real or not: by construction, we assume that E i ( ) is an analytic function of and therefore will never be able to reproduce the dependence of real eigenvalues meeting and moving into the complex plane.For this reason, we next present an alternative approach.

Expansion in energy separation
The problem noted above that arises when the eigenstates of H 0 are not exactly degenerate can be reconciled by performing an expansion in the energy gap between the almost degenerate subspace and the rest of spectrum.More formally, we generalize the effective Schrödinger equation (A2) by introduction of a dimensionless parameter λ, where We assume that we can expand E i (λ) in a power series of λ, but treat its -dependence exactly, and show that it stays real to all orders in λ.Since λ multiplies G in Eq. (A9), this expansion is controlled by the gap to the other states of the spectrum being large (compared to V ).The argument proceeds similar to the one above: the zeroth order contribution, T 0 [E i (λ)] = E i (0), is obtained from diagonalization of Eq. (A6) and as such real.
One can compute With the same arguments as above, this implies that T N [E i (λ)] will stay real for any N .Of course, the perturbative approach is expected to break down when the gap between the almost degenerate subspace and another part of the spectrum with different eigenvalue under η closes since G will develop a pole.

Approximate orthogonality
Above, we have argued that the effective Hamiltonians in Eqs.(A3) and (A9) will be Hermitian if the eigenvalues of η are identical in the almost degenerate subspace.This does not only have consequences for the reality of the eigenvalues, but also for their mutual orthogonality.
To first order in and zeroth order in λ, i.e., to leading order in the limit of a large gap to the excited states, the effective Hamiltonian is also independent of E. Therefore, the projections P |Ψ i ( ) , i = 1, 2, . . .n, are obtained as eigenstates of the same Hermitian Hamiltonian and, as such, orthogonal.Naturally, this does not mean that |Ψ i ( ) are orthogonal in the full Hilbert space; however, the differences between the full and the projected states, |Ψ i ( ) −P |Ψ i ( ) = Q |Ψ i ( ) , are also suppressed in the limit of large energetic separation to the rest of the spectrum since [123] as stated in the main text.

Appendix B: Interplay between X-Cube foliation and metric operators
In the main paper, we noted that the ground states of the X-cube model all have the same eigenvalue under our choice of metric operator η in Eq. ( 5), provided all lengths are even.This is because η can be assembled by a collection of stabilizers.While η cannot be assembled by stabilizers on a system with odd lengths, it is known that the X-cube model exhibits foliated fracton order [119], which implies that an L × L × L X-cube model ground state can be enlarged to a ground state of an L × L × L + 1 model by the attachment of an L × L toric code ground state and the application of local unitary operators.If L is even, then the original X-cube ground states and the toric code ground states will all have the same eigenvalue under η.Because of this, one may suspect that the resulting L × L × L + 1 ground states may also have the same eigenvalue under the appropriately enlarged η.However, as we will show, the process of attaching the two states and applying local unitary operators yields an L × L × L + 1 state that is not an eigenstate of the enlarged η.
We first describe the process of adding an extra layer to the X-cube model, illustrated in Fig. 4. We begin with an L × L × L X-cube ground state, |ψ X , an L × L toric code ground state, |ψ T C , and a collection of L 2 additional qubits initialized in the |0 state, Z |0 = |0 .The statement of foliated fracton order is that an L × L × L + 1 X-cube ground state, |ψ X , can be written as where S is a series of local unitary transformations, which in our case is given by a collection of CNOT gates [119].
This foliation allows us to deduce the behavior of η in Eq. ( 5) applied to |ψ X based on the action of S † ηS on the three constituent states, assuming L is even.This behavior is dependent on the form of η.We first begin with an analysis of η Z ≡ i Z i .Carrying out the corresponding CNOT gate transformations, we see in Fig. 5 that the action of S † η Z S on the original X-cube ground state is not simply the product of all Z i operators-some sites are missing in a way that cannot simply be compensated by a product of stabilizers; this means that |ψ will generally not be an eigenstate of S † η Z S. Carrying this through with η X = i X i and η Y = i Y i yields a similar result.In accordance with the analysis of the main text, we conclude that not every ground state of an even-by-even-by-odd X-cube model will be an eigenstate of η, as the foliation process complicates the behavior of the metric operator.
One can take this L × L × L + 1 model and attach additional toric code layers in either of the two remaining directions, and an identical analysis implies that evenby-odd-by-odd and odd-by-odd-by-odd ground states will not all have the same eigenvalue under η.Of course, one can add another toric code layer to give an L × L × L + 2 model, in which case the metric operator does decompose nicely into the metric operators on the constituent ground states.
FIG. 4. The size of an X-cube model ground state can be increased along one axis by adding a layer of toric code (blue) and an additional set of qubits initialized in the |0 state (red).A series of CNOT gates are applied to this tensor product of states to yield an enlarged X-cube model ground state.The application of the CNOT gates is shown above, with arrows pointing from the control to the target qubit.The CNOT gates are applied in two steps-the gates in the top diagram are used first, then the gates in the bottom diagram are applied.FIG. 5.The behavior of the metric operator η on an L × L × L + 1 X-cube model ground state can be calculated by an effective operator, S † ηS, acting on the exfoliated parts of the X-cube model.Here, we show the action of S † ηZ S (top) and S † ηX S (bottom), where the effective operator is the product of Z (X) operators on all the dark sites.The tensor product of the exfoliated parts of the X-cube ground state is not generally an eigenstate of the corresponding effective operators.The decomposition of ηY is identical to that of ηX .polynomials f and g.As stated in the main text, these codes admit a large set of possible pseudo-Hermitian perturbations that leave the code subspace real: in analogy to Eq. ( 5), a very natural set of choices for the metric operator η is given by Eq. ( 14).
Since all the stabilizers in Haah's cubic codes are mu-tually commuting, all ground states have the same eigenvalue under η, provided η can be assembled by stabilizers.
For the toric code, the X-cube and checkerboard model discussed in the main text, it is straightforward both to find the combination of stabilizers that yield η on a lattice with an even number of sites in all directions, and to show that η cannot be made of stabilizers on any other lattice.For Haah's codes, the more complex form of the stabilizers makes the analysis more demanding, but possible using the polynomial representation of stabilizers [121].
Using the same conventions as in Sec.IV A, the metric operators in Eq. ( 14) can be written as (C1) We will first consider η = Z(h, h).For stabilizers Z(f, g), a choice of covering (i.e., a product of stabilizers at different points) can be specified by a covering polynomial k, with the covering given by Z(kf, kg).For example, if k = 1 + x, then the covering Z(kf, kg) would consist of the product of two stabilizers-one at the origin, and one at (x, y, z) = (1, 0, 0).Therefore, the question of whether η can be assembled from stabilizers is equivalent to the question of whether h = kf = kg for some polynomial k.Mathematically, this factorization takes place in the quotient ring P/I, where P is the ring of polynomials of three variables with coefficients over F 2 , and I is the ideal generated by x Lx + 1, y Ly + 1, and z Lz + 1.This quotienting procedure imposes the periodic boundary conditions of the model.
We calculate this factorization with the computer algebra system SageMath.Generically, this factorization procedure will yield two different coverings, h = k f f = k g g.To determine whether these two coverings are compatible, we calculate whether k f + k g can be separated into two polynomials d f + d g , where d f ∈ (I : f ) and d g ∈ (I : g), where (I : f ) is the colon ideal, (I : f ) = {p ∈ P : pf ∈ I}.This is equivalent to checking whether k f + k g belongs to the ideal generated by (I : f ) ∪ (I : g).If such a separation exists, then k f + d f = k g + d g ≡ k, and h = kf = kg in P/I.This covering may not be unique, as k + d f g also works as a covering, where d f g ∈ (I : f ) ∩ (I : g); however, for the purposes of understanding the behavior of non-Hermitian perturbations, we are only interested in the existence of such a covering.We note that this procedure should always be able to find a covering k if it exists, so if a decomposition k f + k g = d f + d g does not exist, it should imply the non-existence of a covering.
Once we have obtained the covering k for Z(h, h), we immediately know that X(h, h) in Eq. (C1) can be assembled from X-stabilizers with the covering k, since X( kḡ, k f ) = X( h, h) = X(h, h).
This calculation is done in SageMath (see supplementary files) for system sizes L x × L y × L z for 1 ≤  I.
L x , L y , L z ≤ 19.Although the existence/non-existence of a covering follows no clear pattern for very small system sizes, we see regular behavior emerge once the system size is larger than 3 × 3 × 3. Specifically, the existence/non-existence of a covering for a certain cubic code is only dependent on whether each length is even or odd and, if it is even, whether it is divisible by 4. This admits 3 3 = 27 different possible classes of system sizes-however, we find that some classes are equivalent in terms of which cubic codes have coverings on them.A full table of this behavior is shown in Table II.We note several trends.On an odd-by-odd-by-odd lattice, none of the 17 cubic codes have code subspaces that stay real under pseudo-Hermitian perturbations.If only a portion of the system lengths are odd, the reality of the code subspace depends on which dimensions have odd lengths, and whether the remaining lengths are divisible by 4. In contrast, if L x and L y are divisible by 4 and L z is even, all the code subspaces stay real under pseudo-Hermitian perturbations.Overall, cubic code 17 is the most unstable to pseudo-Hermitian perturbations, in that its code subspaces will become complex for almost all system sizes.In contrast, cubic code 7 has the most stable code subspace.There are some groups of codes with the same sensitivity to system sizes.If we consider codes with the same behavior up to a lattice rotation, these groups are (11,12,14,15), (5,8,10,16), and (2,3,6,9).It is interesting to note that, with the exception of cubic code 16, all codes within a group transform the same under entanglement renormalization [122].
While our results are purely numerical, an analytic verification of these trends for all system sizes is likely possible if one was to manually follow the factorization processes carried out in SageMath and show that their conclusions are only sensitive to the system sizes' evenness/oddness and whether they are divisible by 4. We do not attempt this, as there are 459 separate cases that must be checked (27 II.The reality of the code subspace of Haah's cubic codes under pseudo-Hermitian perturbations are highly sensitive to the system size.The reality of the subspace depends on whether each dimension length is odd (o), even and divisible by 4 (E), or even and not divisible by 4 (e).Shown are all 17 of Haah's cubic codes and the dependence of the code subspace stability on the system size.Codes with identical dependencies have been grouped together, and some have been redefined by a spatial rotation.The trends listed have been confirmed numerically to hold from system sizes 3×3×3 to 19 × 19 × 19.

FIG. 6 .
FIG. 6.The stabilizers of Haah's cubic codes correspond to cube operators, with generically different operators on each vertex as labelled by A-D and A -D .The operators at each vertex for the 17 different cubic codes are given in TableI.
possible system sizes for the 17 codes), and instead analyze the numerical results which show clear trends up to 19 × 19 × 19 lattices.A B C D A B C D 1 ZI ZZ IZ ZI IZ II ZI IZ 2 IZ ZZ ZI ZI ZI ZZ IZ ZI 3 IZ ZZ ZZ ZI ZZ II IZ IZ 4 IZ ZZ ZI ZI IZ II IZ ZI 5 ZI ZZ II ZZ ZI II IZ IZ 6 ZI II ZI ZZ IZ ZZ II IZ 7 ZI ZZ ZI IZ IZ II II ZZ 8 ZI ZI IZ ZZ IZ II IZ ZI 9 ZI IZ ZZ ZZ IZ ZZ II IZ 10 ZI IZ ZI ZZ IZ ZZ ZI ZI 11 ZI ZZ II IZ ZI II IZ ZZ 12 ZI IZ ZZ ZZ ZI II II IZ 13 ZI ZZ IZ ZI IZ II II ZZ 14 ZI IZ ZZ ZZ IZ II ZZ IZ 15 ZI IZ II ZZ IZ ZZ II ZI 16 ZI ZI II IZ IZ ZZ II ZZ 17 ZI ZZ IZ ZI IZ ZI ZI ZZ TABLEI.The Z stabilizers for Haah's 17 CSS cubic codes, defined on the eight vertices of a cube, with vertices labeled according to Fig.6.The X stabilizers are obtained by exchanging A ↔ A , and likewise for the other vertices, and by exchanging the two Pauli spins on each site.
e × e × o e × o × o E × o × o o × e × E o × E × e o × e × e o × o × E o × o × e o × E × o o × e × o o × o × o TABLE