Competing topological orders in three dimensions: X-cube versus toric code

We study the competition between two different topological orders in three dimensions by considering the X-cube model and the three-dimensional toric code. The corresponding Hamiltonian can be decomposed into two commuting parts, one of which displays a self-dual spectrum. To determine the phase diagram, we compute the high-order series expansions of the ground-state energy in all limiting cases. Apart from the topological order related to the toric code and the fractonic order related to the X-cube model, we found two new phases which are adiabatically connected to classical limits with nontrivial sub-extensive degeneracies. All phase transitions are found to be first order.


Introduction
Quantum systems with topological order are an important research field due to their intriguing physical properties as well as their potential relevance for quantum technological applications. In two dimensions, these systems are essentially characterized by long-range entanglement and exotic excitations called anyons [1,2], which have quantum exchange statistics distinct from bosons and fermions (see Ref. [3] for a review). These anyonic particles are at the heart of topological quantum computing [4,5] and have been observed unambiguously in quantum Hall systems only very recently [6,7]. Other examples of physical systems where topologicallyordered states play an important role are frustrated quantum magnets and synthetic matter in quantum-optical platforms [8][9][10][11][12][13][14][15].
During these last years, topological order in three dimensions (3D) gained a lot of interest. Some properties such as a topology-dependent ground-state degeneracy are very similar to two-dimensional topological order. Furthermore, in three-dimensional topologically-ordered systems point-like anyonic excitations are excluded but nontrivial statistics can be found for extended objects such as membranes. However, in 3D, one must distinguish between two main categories of topologically-ordered long-range entangled ground states [16,17] depending on whether their degeneracy is finite [18][19][20][21] or sub-extensive with the system size [22][23][24][25][26][27][28][29][30][31][32][33][34][35]. The topological order for systems with sub-extensive ground-state degeneracy is called fracton order. One defining characteristic of fracton phases is that their elementary excitations have a restricted mobility under the action of local operators so that they are considered as attractive candidates for 3D quantum memories [26].
Paradigmatic examples of these two categories of topological order are the 3D toric code (TC) model [19,20], which is a direct extension of the celebrated model introduced by Kitaev in two dimensions [4] for fault-tolerant quantum computation and has a finite ground-state degeneracy, and the X-cube model (XC) proposed by Vijay, Haah, and Fu [28], which has a sub-extensive ground-state degeneracy.
In the absence of any local order parameter, the study of transitions between topological quantum phases of matter is a challenging problem. In two dimensions, the concept of anyon condensation [36] and topological symmetry breaking provides a general framework [37] to understand some of these transitions [38][39][40]. Other salient examples are the Kitaev's honeycomb model which exhibits a transition between an achiral topological phase and a chiral topological phase [41], the string-net model [42] that allows to investigate the competition between different topological achiral phases obeying the same fusion rules [43][44][45], or multilayer systems [35,46]. To our knowledge, similar studies are still missing in 3D.
The goal of the present work is to investigate the competition between two different types of topological orders by considering the interplay between the TC and the XC phases on the cubic lattice. As shown below, the corresponding Hamiltonian can be split into two commuting parts that are analyzed separately. Interestingly, one of these parts has a self-dual spectrum. We determine the phase diagram in the full four-dimensional parameter space using this self-duality as well as high-order series expansions of the ground-state energy in all limiting cases. Apart from the TC and XC phases, the phase diagram displays two additional phases, dubbed X-and Z-phases, that are connected to classical limits with sub-extensive degeneracies. All phase transitions are found to be first order.
The paper is organized as follows: in Sec. 2 we introduce the model and we recall the main properties of the limiting cases. Then, we show that one can recast the full Hamiltonian in two sets of commuting operators allowing for a simpler analysis of the phase diagram which is discussed in Sec. 3. We conclude our findings and give some perspectives in Sec. 4.

Model
We consider the interplay between the TC and the XC. Microscopic degrees of freedom are spins-1/2 located on the links of the cubic lattice (see Fig. 1). The corresponding Hamiltonian is given by with where and represent elementary faces and cubes of the cubic lattice, whereas and label the two different vertex operators of the TC and XC, respectively (see Fig. 1 for illustration). Without loss of generality, we assume non-negative couplings J , J , J , and J for the rest of the paper. Denoting by σ α i the usual Pauli matrices with α = x, y, z acting on the link i, the operators of the TC are defined as where the products run over the six spins (four spins) of ( ). The XC operators are where the products run over the four spins (twelve spins) of ( ). The four operators (4)-(7) have eigenvalues ±1. We stress that all pairs of operators commute unless they have different spin flavors and share an odd number of spins. Hence, the Hamiltonian (1) is not exactly solvable for arbitrary couplings but there are some limiting cases where H can be solved analytically.

Limiting cases
These limits are connected to the four different phases present in the phase diagram discussed in the next section. TC phase: For J = J = 0, the system reduces to the pure TC which is exactly solvable [19][20][21] since [X , Z ] = 0 for all and . The TC has a finite ground-state degeneracy which only depends on the genus of the 3D surface [16,18], e.g., the degeneracy is 2 3 on a 3-torus. These ground states have a finite topological entropy and can be seen as 3D generalizations of the loop gas ground state of the conventional 2D toric code [19][20][21]. Furthermore, gapped elementary excitations correspond to point and spatially extended particles displaying a semionic mutual statistics. In the following, we call the phase adiabatically connected to the limit where J = J = 0, the TC-phase.
XC phase: For J = J = 0, the system corresponds to the exactly solvable XC [28]. Its ground-state degeneracy depends not only on the topology but also on the geometry of the system [30,47]. For instance, on a 3-torus with linear extensions L 1 , L 2 , and L 3 , the ground-state degeneracy equals 2 2(L 1 +L 2 +L 3 )−3 [17,30]. The ground states of the XC have a topological entanglement entropy which scales sub-extensively with the linear system size [48,49] and may again be seen as a generalized loop gas [50]. The XC is known to display type-I fracton topological order [28], i.e., its elementary excitations are either immobile or have a dispersion with dimensional reduction upon the action of local operators. The immobile fracton excitation of the XC corresponds to a single cube excitation with eigenvalue −1 of one Z . The other elementary excitations of the XC are one-dimensional particles related to pairs of X operators with eigenvalue −1 at the same vertex. We call the phase adiabatically connected to the limit where J = J = 0, the XC-phase.
Classical limits: In the limit where J = J = 0 (J = J = 0), the only operators in the Hamiltonian are products of σ x i (σ z i ). Eigenstates of H are thus trivial product states and the ground states are all states with eigenvalues +1 for these operators. For periodic boundary conditions (3-torus with linear extension L), one finds a non-trivial ground-state degeneracy 2 3L 2 for J = J = 0 and 2 3(L 2 −L−1) for J = J = 0. For J = J = 0, the different ground states are distinguished by a set of non-local commuting operators defined as product of σ z i on straight lines. For J = J = 0, ground states are distinguished by a set of nonlocal commuting operators defined on non-contractible tubes and planar membranes (details about these operators are given in Appendix A). As a direct consequence, the ground-state degeneracy of these phases is expected to be robust with respect to small perturbations in the couplings J , J or J , J . In the following, phases connected to the two classical limits will be called X-and Z-phase for obvious reasons.

Exact decomposition and self-duality
The essential ingredient to determine the ground-state phase diagram in the full parameter space relies on another decomposition of the Hamiltonian.
and, as can be easily checked, [H A , H B ] = 0. Interestingly, the spectrum of H A is exactly selfdual (up to degeneracies) so that the phase diagram of H A must be symmetric with respect to the self-dual point J = J . To prove this self-duality, let us introduce the following pseudospin-1/2 operators τ z v = X defined on the vertices of the original cubic lattice Λ. Then, it is easy to see that Z acts like v∈ τ x v so that H A becomes which describes a transverse-field spin-1/2 model with eight-spin interactions on a cubic lattice.
Similarly, if one introduces τ z v = Z defined on vertices of the cubic lattice Λ spanned by the center of each elementary cube of Λ, the operator X acts like v∈ τ Since Λ and Λ are both cubic lattices, the spectrum of H τ A and the one from H τ A are obtained from the other by exchanging J ↔ J . As a consequence, the spectrum of H A is invariant under this exchange and, hence, self-dual. Of course, the mapping described above does not preserve the degeneracies of the spectrum. Thus, the self-duality of H A is only exact, up to degeneracies. For a very similar discussion in two dimensions, see Refs. [51,52].
In contrast, H B is not self-dual as can be directly seen in the series expansions of the ground-state energy given in Appendix B.

Phase diagram
The decomposition of H into two commuting parts ([H A , H B ] = 0) allows one to build the full phase diagram from the ones of H A and H B , separately. We emphasize that this decoupling implies that the phase diagram only depends on the two ratios J /J and J /J driving the transition of H A and H B , respectively. Similarly, assuming the existence of a unique transition point in the phase diagram of H B , we determined its position by extrapolating the crossing point of high-order series expansions for the ground-state energy in both limiting cases J J and J J (see Appendix B). We found a transition point at J /J 1.012 ≡ η B (see Fig. 2). We stress that η A and η B are associated to first-order transitions. Hence, assuming a unique transition point for H A and H B , we obtain the complete phase diagram of H shown in Fig. 3 which contains four distinct phases separated by first-order transition lines. For small J /J and J /J the system is in the TC-phase. For large J /J and J /J , one gets the fractonic XC-phase. In the limit J /J 1 and J /J 1, one finds a phase essentially driven by the two operators consisting of σ x i , i.e., the X-phase. Similarly, when J /J 1 and J /J 1, the phase is mainly determined by the two operators consisting of σ z i , namely, the Z-phase.
Interestingly, we emphasize that a direct transition between the XC-and the TC-phase requires a fine-tuning of the parameters.

Conclusions
In this work we investigated the competition between the two most paradigmatic representatives of 3D topological order. An exact decomposition of the system allows for a quantitative determination of the ground-state phase diagram in the full parameter space. Apart from the TC-and XC-phase, the phase diagram contains two additional phases, the X-and Z-phases, that are connected to limiting cases where operators with either σ x i or σ z i Pauli matrices dominate. In the purely classical limits, we find non-trivial sub-extensive ground-state degeneracies which are robust perturbatively. However, a better understanding of the quantum nature of the X-and Z-phases would be valuable but it is beyond the scope of the present work.
In the derivation of the full phase diagram we assumed that H A and H B display a single phase transition which is found to be first order. Unfortunately, the existence of intermediate phases can not be ruled out by our approach. Although we consider unlikely the existence of intermediate phases, an unbiased numerical investigation would be valuable.
KPS acknowledges financial support by the German Science Foundation (DFG) through the grant SCHM 2511/11-1.

A Ground-state degeneracy of X-and Z-phases
In this appendix we compute the ground-state degeneracy (GSD) of the X-and Z-phase in the classical limit on a 3-torus.
In order to access the GSD of the X-phase and the Z-phase in the classical limit on a 3-torus of dimensions L 1 × L 2 × L 3 we describe the operators whose eigenvalues distinguish between the different ground-states. Importantly, these operators are non local. As shown below, we found that log 2 GSD obeys an area law, contrasting with the linear behaviour of the XC-phase [17,30] and the constant value of the TC-phase [19][20][21].

A.1 X-phase
In the X-phase the relevant operators are products of σ z i operators acting on straight lines and forming non-contractible loops (see Fig. 4 for illustration). Since all these operators are independent and mutually commute, there are loops on the 3-torus. We checked these expressions numerically, on (small) finite systems.

A.2 Z-phase
In the Z-phase there are two types of operators whose eigenvalues distinguish between the different ground-states, planes of σ x i which are precisely halfway in between the lattice planes and non-contractible tubes of σ x i which correspond exactly to the action of X -operators along a straight line. Two plane-operators correspond exactly to the product of the tube operators inside a layer, so not all tube operators inside a layer are independent. These non-local tube and plane operators are illustrated in Fig. 5. In total, we have i L i "plane" operators and i<j (L i − 1)(L j − 1) "tube" operators. Accordingly, one gets

B Series expansions
In this Appendix, we give high-order series for the ground-state energy per vertex and briefly comment on the methods to derive them.
We computed several high-order series expansions using the Löwdin method [53]. This method has been applied successfully in related contexts like the robustness of effective cluster states in measurement-based quantum computation [54] as well as topological string-net phases [55]. Furthermore, the application of the Löwdin method to perturbed topological models is well described in Ref. [56] and we therefore focus on the central aspects for the current problem. Here, we perform the series expansion separately for H A and H B and calculate the ground-state energy for all perturbative limits of H A and H B . The series expansions of the ground-state energy of the original Hamiltonian (1) are then simply given by the sum of the ground-state energies of H A and H B in the appropriate limits. We stress that the obtained series for the ground-state energy are valid for the entire ground-state manifold in all considered limits. This is a direct consequence of the fact that ground states of the same manifold are only connected by non-local operators so that the degeneracy remains intact at any finite order of perturbation theory. Accordingly, the calculation of the ground-state energy up to any finite perturbation order can be performed on any state in the ground-state manifold without loss of generality.
The calculation is most efficiently done via a full graph decomposition using a linked-cluster expansion. Since H A and H B only contain multi-spin interactions linking multiple degrees of freedom, a natural formulation of the linked-cluster expansion is therefore performed in terms of hypergraphs [57]. A hypergraph is a generalization of a graph where edges can link more than two vertices. Technically, we generate all linked subclusters up to a given size [58,59] and sort them into isomorphism classes of hypergraphs using their König representation [60,61]. During this procedure non-contributing subclusters are discarded as early as possible using heuristics adapted from Refs. [62,63]. This allows us to determine the ground-state energy per vertex for H A and H B in the four perturbative limits up to order 10. The corresponding series for the full problem H A + H B can then be straightforwardly extracted.
For the ground-state energy of H A per vertex up to order ten, we find Because of the exact self-duality the series in the opposite limit J J is easily obtained by exchanging J and J in the above expression.
For the ground-state energy of H B per vertex up to order ten we obtain