Absence of Finite Temperature Phase Transitions in the X-Cube Model and its $\mathbb{Z}_{p}$ Generalization

TL;DR

This paper analyzes the thermal behavior of the X-Cube model and its $ ext{Z}_p$ generalization, demonstrating the absence of finite temperature phase transitions and suggesting their classical analogs lack glassy dynamics.

Contribution

It shows that the X-Cube model and its $ ext{Z}_p$ extension do not exhibit finite temperature phase transitions and are thermally fragile, extending understanding of fracton models.

Findings

No finite temperature phase transitions in these models.

Thermal fluctuations do not enable size-dependent autocorrelations.

Models map to classical Ising or clock chains without glassy dynamics.

Abstract

We investigate thermal properties of the X-Cube model and its $\mathbb{Z}_{p}$ `clock-type' ($p$X-Cube) extension. In the latter, the elementary spin-1/2 operators of the X-Cube model are replaced by elements of the Weyl algebra. We study different boundary condition realizations of these models and analyze their finite temperature dynamics and thermodynamics. We find that (i) no finite temperature phase transitions occur in these systems. In tandem, employing bond-algebraic dualities, we show that for Glauber type solvable baths, (ii) thermal fluctuations might not enable system size dependent time autocorrelations at all positive temperatures (i.e., they are thermally fragile). Qualitatively, our results demonstrate that similar to Kitaev's Toric code model, the X-Cube model (and its $p$-state clock-type descendants) may be mapped to simple classical Ising ($p$-state clock) chains in which neither phase transitions nor anomalously slow glassy dynamics might appear.

Figures

Equations

$$A_c\equiv \prod _{n\in \partial c}{\sigma }_n^x,$$

$$B_v^x\equiv {\sigma }_j^z{\sigma }_n^z{\sigma }_{\ell }^z{\sigma }_m^z,B_v^y\equiv {\sigma }_i^z{\sigma }_n^z{\sigma }_k^z{\sigma }_m^z,B_v^z\equiv {\sigma }_i^z{\sigma }_j^z{\sigma }_k^z{\sigma }_{\ell }^z.$$

$$B_v^x{B}_v^y{B}_v^z=\text{mathds}1.$$

$$B_v^x{B}_v^y=B_v^z,B_v^x{B}_v^z=B_v^y,B_v^y{B}_v^z=B_v^x.$$

$$[A_c,B_v^{\mu }]=0,\forall c,v,\mu .$$

$$H=-a\sum _c{A}_c-b\sum _{\mu ,v}{B}_v^{\mu },$$

$$\mathcal{Z}=\mathrm{Tr}\left[\exp \left(\beta a\sum _c{A}_c+\beta b\sum _{\mu ,v}{B}_v^{\mu }\right)\right],$$

$$\mathcal{Z}=\mathrm{Tr}\left[\prod _c\left(\exp (\beta a{A}_c)\right)\prod _{\mu ,v}\left(\exp (\beta b{B}_v^{\mu })\right)\right].$$

$$\begin{array}{rl}\exp (\beta a{A}_c)& =\text{mathds}1+\beta a{A}_c+\frac{1}{2}(\beta a{)}^2\text{mathds}1+\frac{1}{6}(\beta a{)}^3{A}_c+\dots \\ & =\text{mathds}1\cosh (\beta a)+A_c\sinh (\beta a).\end{array}$$

$$\prod _c\exp (\beta a{A}_c)={\mathsf{C}}_a^{L^3}\prod _c\left[\text{mathds}1+A_c{\mathsf{T}}_a\right],$$

$$\prod _c\left[\text{mathds}1+A_c{\mathsf{T}}_a\right]=\text{mathds}1+\sum _c{A}_c{\mathsf{T}}_a+\sum _{c<d}{A}_c{A}_d{\mathsf{T}}_a^2+\sum _{c<d<e}{A}_c{A}_d{A}_e{\mathsf{T}}_a^3+\dots .$$

$$f(\beta )=-\frac{1}{\beta L^3}\log \mathcal{Z}$$

$$\begin{array}{rl}{\mathcal{Z}}_{\text{Open}}& ={\mathsf{C}}_a^{L^3}\mathrm{Tr}\left[\prod _c\left[\text{mathds}1+A_c{\mathsf{T}}_a\right]\left(\prod _{\mu ,v}\exp (\beta b{B}_v^{\mu })\right)\right]\\ & ={\mathsf{C}}_a^{L^3}\mathrm{Tr}\left[\prod _{\mu ,v}\exp (\beta b{B}_v^{\mu })\right],\end{array}$$

$$\begin{array}{rl}\prod _{\mu ,v}\exp (\beta b{B}_v^{\mu })& ={\mathsf{C}}_b^{3(L-1{)}^3}\prod _{\mu ,v}\left[\text{mathds}1+B_v^{\mu }{\mathsf{T}}_b\right]\\ & ={\mathsf{C}}_b^{3(L-1{)}^3}\left[\text{mathds}1+\sum _{\mu ,v}{B}_v^{\mu }{\mathsf{T}}_b+\dots \right],\end{array}$$

$$\prod _{\mu ,v}[\text{mathds}1+B_v^{\mu }{\mathsf{T}}_b]=\left[\sum _{n=0}^{(L-1{)}^3}\left(\genfrac{}{}{0ex}{}{(L-1{)}^3}{n}\right){\mathsf{T}}_b^{3n}\right]\text{mathds}1+\text{t.t.},$$

$$\mathrm{Tr}\left[\prod _{\mu ,v}\exp (\beta b{B}_v^{\mu })\right]={\mathsf{C}}_b^{3(L-1{)}^3}{\left[1+{\mathsf{T}}_b^3\right]}^{(L-1{)}^3}\mathrm{Tr}[\text{mathds}1].$$

$${\mathcal{Z}}_{\text{Open}}=2^{3L^3+6L^2+3L}{\mathsf{C}}_a^{L^3}{\left[{\mathsf{C}}_b^3+{\mathsf{S}}_b^3\right]}^{(L-1{)}^3}.$$

$${\mathsf{C}}_b^3+{\mathsf{S}}_b^3=\frac{1}{4}(e^{3\beta b}+3e^{-\beta b}).$$

$${\mathcal{Z}}_{\text{Open}}=2^{L^3+12L^2-3L+2}{\mathsf{C}}_a^{L^3}(e^{3\beta b}+3e^{-\beta b}{)}^{(L-1{)}^3}.$$

$$f_{\text{Open}}=-\frac{1}{\beta }\left[\frac{3L^3+6L^2+3L}{L^3}\log 2+\log {\mathsf{C}}_a+\frac{(L-1{)}^3}{L^3}\log ({\mathsf{C}}_b^3+{\mathsf{S}}_b^3)\right].$$

$$\begin{array}{rl}{A}_c\to r_m,& 1\le m\le L^3,\\ B_v^x\to s_1^n,B_v^y\to s_2^n,B_v^z& \to s_1^n{s}_2^n,1\le n\le (L-1{)}^3.\end{array}$$

$${\mathcal{Z}}_{\text{Open}}=2^{\bar{N}}\sum _{\{r_m,s_j^n\}}\prod _{m=1}^{L^3}\prod _{n=1}^{(L-1{)}^3}\exp [\beta a{r}_m+\beta b(s_1^n+s_2^n+s_1^n{s}_2^n)].$$

$${\mathcal{Z}}_{\text{Open}}=2^{\bar{N}}(2{\mathsf{C}}_a{)}^{L^3}(e^{3\beta b}+3e^{-\beta b}{)}^{(L-1{)}^3}.$$

$$B_v^x\to s_1^n{s}_2^n,B_v^y\to s_2^n{s}_3^n,B_v^z\to s_3^n{s}_1^n.$$

$$⟨A_{c_1}\dots A_{c_{\mathcal{M}}}{B}_{v_1}^{{\mu }_1}\dots B_{v_{\mathcal{N}}}^{{\mu }_{\mathcal{N}}}⟩=\frac{\mathrm{Tr}[e^{-\beta H}{A}_{c_1}\dots A_{c_{\mathcal{M}}}{B}_{v_1}^{{\mu }_1}\dots B_{v_{\mathcal{N}}}^{{\mu }_{\mathcal{N}}}]}{\mathrm{Tr}[e^{-\beta H}]},$$

$$⟨\mathbb{AB}⟩=\frac{\mathrm{Tr}\left[{\prod }_c[\text{mathds}1+A_c{\mathsf{T}}_a]{\prod }_{\mu ,v}[\text{mathds}1+B_v^{\mu }{\mathsf{T}}_b]\mathbb{AB}\right]}{\mathrm{Tr}\left[{\prod }_c[\text{mathds}1+A_c{\mathsf{T}}_a]{\prod }_{\mu ,v}[\text{mathds}1+B_v^{\mu }{\mathsf{T}}_b]\right]},$$

$$⟨\mathbb{AB}⟩=\frac{{\mathsf{T}}_a^{\mathcal{M}}\mathrm{Tr}\left[{\prod }_{\mu ,v}[\text{mathds}1+B_v^{\mu }{\mathsf{T}}_b]\mathbb{B}\right]}{\mathrm{Tr}\left[{\prod }_{\mu ,v}[\text{mathds}1+B_v^{\mu }{\mathsf{T}}_b]\right]}.$$

$$⟨\mathbb{AB}⟩={\mathsf{T}}_a^{\mathcal{M}}\frac{{\left({\mathsf{T}}_b+{\mathsf{T}}_b^2\right)}^{\mathcal{N}}\left({\sum }_{n=0}^{(L-1{)}^3-\mathcal{N}}\left(\genfrac{}{}{0ex}{}{(L-1{)}^3-\mathcal{N}}{n}\right){\mathsf{T}}_b^{3n}\right)\mathrm{Tr}[\text{mathds}1]}{\left({\sum }_{n=0}^{(L-1{)}^3}\left(\genfrac{}{}{0ex}{}{(L-1{)}^3}{n}\right){\mathsf{T}}_b^{3n}\right)\mathrm{Tr}[\text{mathds}1]}.$$

$$⟨\mathbb{AB}⟩={\mathsf{T}}_a^{\mathcal{M}}{\left(\frac{{\mathsf{T}}_b+{\mathsf{T}}_b^2}{1+{\mathsf{T}}_b^3}\right)}^{\mathcal{N}}.$$

$$⟨A_{c_1}\dots A_{c_{\mathcal{M}}}{B}_{v_1}^{{\mu }_1}\dots B_{v_{\mathcal{N}}}^{{\mu }_{\mathcal{N}}}⟩=⟨A_{c_1}⟩\dots ⟨A_{c_{\mathcal{M}}}⟩⟨B_{v_1}^{{\mu }_1}⟩\dots ⟨B_{v_{\mathcal{N}}}^{{\mu }_{\mathcal{N}}}⟩$$

Full text

††journal: Annals of Physics

1. 1 Introduction

2. 2 Main Results of This Work

3. 3 Outline

4. 4 General Elements of the X-Cube Model and its Partition Function

5. 5 Open Boundary Conditions

6. 5.1 Partition Function

7. 5.1.1 High Temperature Series Expansion

8. 5.1.2 Bond Algebraic Duality

9. 5.2 Cube-Star Correlation Functions

10. 6 Cylindrical Boundary Conditions

11. 6.1 Partition Function

12. 6.1.1 High Temperature Series Expansion

13. 6.1.2 Bond Algebraic Duality

14. 6.2 Cube-Star Correlation Functions

15. 7 Periodic Boundary Conditions

16. 7.1 Partition Function

17. 7.2 Absence of Phase Transitions

18. 8 Dynamics of the X-Cube Model at Finite Temperature

19. 9 Equilibrium Thermodynamics of the $p$X-Cube Model

20. 9.1 Open Boundary Conditions

21. 9.2 Cylindrical Boundary Conditions

22. 9.3 Large-$p$ Limit

23. 10 Implications of Dualities on the Nature of Fracton Excitations

24. 11 Conclusion and Outlook

25. 12 Acknowledgments

26. A Ground States and Low Energy Excitations: A Review and General Remarks

1 Introduction

Notwithstanding the triumphs of the Landau symmetry breaking paradigm [1, 2], it is not powerful enough to describe numerous physical systems. A detailed pedagogical discussion of these non-Landau symmetry breaking aspects appears in [2]. Let us briefly mention several of these. In gauge theories, including those describing the fundamental interactions, symmetry breaking is prohibited by Elitzur’s theorem [3, 4]. Consequently, in the absence of matter fields, only expectation values associated with closed loops (so-called “Wilson loops”) may acquire finite expectation values and differentiate between the various phases [5, 6, 7]. Another important theory in which symmetry breaking cannot occur is the classical two-dimensional XY model [8, 9]. Herein, the binding and unbinding of topological excitations (vortices) characterizes the different phases of the system.

Along very different lines, the remarkably accurate quantization of the conductance plateau in quantum Hall systems can stem from topological invariance [10]. Motivated by these and related considerations and the prospect of spin-liquids [11, 12], the notion of “topological order” has been introduced, e.g., [13, 14]. This endeavor has been bolstered by the prospect of employing topological matter for quantum computation and by the construction of elegant soluble spin models in which basic notions of topological quantum information come to life [15, 16]. Indeed, these and other early investigations, e.g., [17, 18, 19] discussed or were partially motivated by the quest of achieving “fault-tolerant” topological quantum hardware. The last decade has also witnessed a flurry of experimental findings of materials that exhibit topological effects [20].

More recently, various quantum spin models believed to exhibit “fracton topological order” [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35] were studied. Fracton topological order is a proposed state of matter in which fundamental excitations, termed “fractons”, exhibit rich behaviors. In particular, fractons may exhibit confinement along certain spatial directions and further display hindered dynamics at low energies [23]. Individual fractons may be somewhat immobile. However, collectively, several fractons may more readily move together in a constrained correlated fashion at low energies. As a consequence of these constraints on their motion, fractons have been believed to exhibit slow, glassy dynamics, ideal for finite temperature quantum memories [21, 23, 36, 37]. For a more general current perspective of related issues in localization, symmetry, and topology, see [38].

Prototypical fracton models display symmetries originally known as “$d$-dimensional Gauge Like symmetries”. These symmetries include, what have been later termed, “Generalized Global Symmetries”, “subsystem symmetries”, “$p$-form symmetries”, or “higher symmmetries” [4, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48]. The relation between $d$-dimensional Gauge Like symmetries and topological quantum order has been established in Refs. [39, 40, 41]. A key ingredient is the generalization of Elitzur’s theorem from conventional gauge (i.e., local ($d=0$)) symmetries to $d>0$ Gauge Like symmetries. These are symmetries that are neither local nor global. For a system residing in $D$ spatial dimensions, the latter symmetries act on a region of intermediate spatial dimensionality $d$ such that $0<d<D$. These symmetries often allow for low-energy excitations that are typically common to $d$ spatial dimensions to appear in the full $D$-dimensional system. Just as excitations in low-dimensional systems eradicate long lived structures, as a consequence of these symmetries, certain topological memories (such as Kitaev’s toric code model [15]) may similarly become susceptible to thermal fluctuations - a phenomenon known as “thermal fragility” [39, 40, 41, 49, 50, 51, 52, 53]. This susceptibility is not to be confused with the existence of topological order at finite temperatures [39, 40]. Imprints of thermal fragility also appear in measures such as the entanglement entropy [54]. Different three-dimensional variants of Kitaev’s original (two-dimensional) toric code model [49] can further exhibit finite temperature robust correlations coexisting with thermally fragile properties [49, 54]. Compass models [55, 56] similarly exhibit topological order and finite temperature transitions and were also suggested as candidate systems for topological memories [57]. A broad perspective on quantum memories is given in [58]. The current paper extends these earlier studies of topological order at finite temperatures to fracton topological order.

Various conservation laws for effective charges and constrained mobility are often associated with pure gauge symmetries and/or $d>0$ Gauge like symmetries [39, 40, 41, 59, 60, 61, 62, 63, 64]. Restricted dynamics are also known to arise in tensorial gauge field theoretic formulations of elasticity. In these theories [62, 63, 64], conservation laws that are captured by the tensorial gauge fields allow for motion of elastic defects only along certain directions (in particular, the well known “glide” of dislocations). Herein, dynamics were included by examining elasticity of a medium in space-time (with the kinetic energy corresponding to the energy of elastic deformations along the temporal direction). In this approach, Noether currents may be computed and conservation laws can be reformulated in terms of a gauge invariance [63]. Such a gauge invariance leads to constraints on the motion of dislocations in solids and allows the derivation of the lower dimensional restricted “glide” motion. Illuminating works [31, 65] related these constraints to those pertaining to fractons.

One particularly popular fracton lattice system is the X-Cube model [22], often regarded as a quintessential example of “type-I fracton topological order” [23]. The model is given by an exactly solvable commuting Pauli Hamiltonian, traditionally defined on a cubic lattice, but the model can alternatively be defined on many different lattices [30]. Compared to other fracton models, such as Haah’s code and the Chamon model [36, 37, 41, 66], the X-Cube model is arguably more intuitive in its structure and excitations, likely leading to its widespread popularity.

The primary focus of this paper is to investigate the properties of the X-Cube model at finite temperatures by exactly solving its partition function for both open and cylindrical boundary conditions. Additionally, we will use these results to investigate corresponding correlation functions and the finite temperature dynamics of the model. We will employ two different approaches towards this end: first, we will solve each partition function utilizing a brute-force trace calculation, making liberal use of the binomial theorem; second, we will apply a bond-algebraic mapping, in which each Pauli operator contained in the Hamiltonian is mapped to a classical Ising-like spin variable in such a way as to preserve the algebra of the model. We will find that these two approaches yield exactly identical solutions. We will also exactly compute the partition function of the $p$X-Cube model, the natural $\mathbb{Z}_{p}$ generalization of the ($\mathbb{Z}_{2}$) X-Cube model, under the same boundary conditions.

Bond-algebraic mappings akin to those contained in this paper have frequently been used to investigate the finite temperature properties of quantum spin models [39, 40, 41, 49, 67, 68, 69, 70, 71, 72, 73, 74]. Using this strategy, the properties of many complex quantum spin models may be determined from the familiar properties of classical Ising and Ising-like models. In particular, the application of bond algebras first demonstrated that the finite temperature partition function of a prototypical example of topological order - Kitaev’s toric code model [15] - a model that exhibits $d=1$ dimensional Gauge Like Symmetries, is identical to that of classical one-dimensional Ising chains [39]. Similarly, it was shown in Ref. [41] that the $XXYYZZ$, or Chamon, model maps onto four decoupled Ising chains and also displays the phenomenon of “dimensional reduction”.

More explicitly, a bond-algebraic mapping utilizes a (generally non-local) unitary transformation $X\rightarrow U^{\dagger}XU$ to map each operator in the Hamiltonian to a simple “classical” product of $z$-Pauli spins, from which properties such as the partition function, correlation functions, and more, can be evaluated. In particular, this approach will readily enable us to study the finite temperature properties of the X-Cube model.

We note that the X-Cube model is most commonly assumed to possess periodic boundary conditions, and that although we will be able to exactly solve the partition function for the case of open and cylindrical boundaries, we will be unable to completely solve (in closed form) for the partition function of the model under fully periodic boundary conditions. However, because the corrections to each partition function due to its boundary conditions will appear at high orders of system size, the corresponding corrections to the free energy will be negligible in the thermodynamic limit. Indeed, in section 7, we will prove that in the thermodynamic limit, the free energy of the X-Cube model with periodic boundary conditions is identically the same as that with open or cylindrical boundaries.

2 Main Results of This Work

A central result of the current work is that:

(1) Much as for Kitaev’s toric code model which may be mapped onto classical Ising spin chains [39, 40, 49], the X-Cube model does not exhibit a finite temperature phase transition, because of the phenomenon of dimensional reduction.

Along similar lines, we find that:

(2) There are, at least, some thermal bath realizations which do not lead to super-Arrhenius glassy dynamics of fractons but rather to those of the conventional activated form found in classical Ising spin chains [49].

Both of these behaviors also appear in the $\mathbb{Z}_{p}$ generalizations of the X-Cube model, i.e., the $p$X-Cube model.

3 Outline

The remainder of this paper is organized as follows.

In section 4, we review the X-Cube model, write its formal partition function, and set up the framework for the high temperature series expansion that we will employ. In section 5, we compute its partition function (section 5.1, Eq. (17)) in the presence of open boundary conditions. This is achieved using both the high temperature series expansion and a bond-algebraic technique that yield identical results. We further compute (section 5.2) general finite temperature correlation functions. In section 6, we perform similar calculations for the X-Cube model endowed with cylindrical boundary conditions (Eq. (39)). These partition functions lead to free energy densities that exhibit no non-analyticities at any finite temperature. In section 7, we turn to the X-Cube model with periodic boundary conditions. Here, the partition function does not admit a simple closed form expression. However, as we demonstrate (and as is generally anticipated for bulk thermodynamic properties), the free energy density of the periodic system is identical to that when open or cylindrical boundary conditions are present. Thus, also in the presence of periodic boundary conditions, the system exhibits no finite temperature transitions. In section 8, using simple duality transformations, we arrive at closed form expressions for the autocorrelations. These calculations demonstrate that for a heat bath generated by dualizing a general Glauber heat bath, the X-Cube model does not exhibit long time correlations. That is, here, notwithstanding the topological character of the theory and the constraints for low energy motion, memory of the initial state may be lost after a finite (system size independent) autocorrelation time. We caution that our results do not exclude the possibility of long lived glassy memory when other heat baths are assumed for the system. In section 9, we discuss the $p$-state clock type generalization of the X-Cube model (the original X-Cube is the $p=2$ realization of this more general model). We compute the partition functions of this $p$X-Cube model under open and cylindrical boundary conditions. In the $p\to\infty$ limit (section 9.3), the $\mathbb{Z}_{p}$ symmetry of the discrete $p$X-Cube models becomes a continuous $U(1)$ symmetry. Given our exact results concerning the lack of finite temperature phase transitions and possible glassy dynamics, we step back, and discuss in section 10 qualitative aspects of low energy motion. We conclude in section 11 with general remarks. We discuss various qualitative features of low energy excitations in the Appendix.

Apart from deriving results for the X-Cube models, our analysis highlights the utility of the bond-algebraic duality techniques. Indeed, although we derived the partitions functions for these models using both traditional high temperature series expansions and bond-algebraic duality mappings, the careful reader may readily appreciate the ease with which bond-algebras may enable us to derive results that may be far more cumbersome to arrive at by the more conventional high temperature series expansions.

4 General Elements of the X-Cube Model and its Partition Function

In this Section, we first briefly review the Hamiltonian of the X-Cube model as introduced by Vijay, Haah, and Fu [22, 23]. We will then turn to our main objective of analyzing the system at finite temperatures. Towards that end, we will then formally write down its partition function invoking the well known high temperature series expansion.

X*-Cube Model —* Consider an $L\times L\times L$ cubic lattice, with qubits (or spin-$\frac{1}{2}$’s) located at each edge $n$ of the lattice (see figure 1). The total number of qubits $N$ in this lattice depends on the choice of boundary conditions, as will be discussed in their respective sections. Each qubit is associated with a two-dimensional Hilbert space $\mathcal{H}_{n}=\mathbb{C}^{2}$. The total state space of the system is then given by $\bigotimes_{n=1}^{N}\mathcal{H}_{n}$ with dimension $2^{N}$.

For each elementary cube $c$ of the lattice, we define the operator $A_{c}$ by:

[TABLE]

where $\sigma^{x}_{n}$ is the $x$-Pauli operator acting on $\mathcal{H}_{n}$. $A_{c}$ is therefore a product of twelve $x$-Pauli operators, each associated with a qubit on one of the twelve edges of the simple cube shown in figure 1a.

In addition, for each vertex $v$ of the lattice, label the six surrounding qubits as $i$, $j$, $k$, $\ell$, $m$, and $n$, as in figure 1b. The operators $B^{\mu}_{v}$, $\mu\in\{x,y,z\}$, are then the four link “stars” defined by:

[TABLE]

Each $B^{\mu}_{v}$ is the product of the four $z$-Pauli operators surrounding $v$ forming a plane perpendicular to the direction $\mu$. From (2), it is quickly seen that:

[TABLE]

Alternatively, using $(\sigma^{z}_{n})^{2}=\mathds{1}$, this can be written as:

[TABLE]

First note that each of the operators $A_{c}$ and $B^{\mu}_{v}$ commute. If the vertex $v$ is not a vertex of the simple cube $c$, then $A_{c}$ and $B^{\mu}_{v}$ act on no common qubits and therefore must commute. If $v$ is a vertex of $c$, then $A_{c}$ and $B^{\mu}_{v}$ act on two common qubits. Since $\sigma^{x}_{n}$ and $\sigma^{z}_{n}$ anticommute, $\sigma^{x}_{n}\sigma^{x}_{m}$ and $\sigma^{z}_{n}\sigma^{z}_{m}$ will commute, and therefore:

[TABLE]

It is also trivially verified that $A_{c}$ and $B^{\mu}_{v}$ are Hermitian operators with eigenvalues $\pm 1$, and that $(A_{c})^{2}=(B^{\mu}_{v})^{2}=\mathds{1}$.

The X-Cube model is defined by the stabilizer Hamiltonian [22, 23]:

[TABLE]

where $a>0$ and $b>0$ are constant parameters. The first sum is performed over all $L^{3}$ simple cubes of the lattice, while the second sum is performed over each vertex $v$ and each of the three cardinal directions $\mu$. The particular vertices to be included in the sum will depend on the choice of boundary conditions.

Symmetries — Hamiltonian (6) is host to several symmetries. $H$ is considered invariant under the symmetry transformation $|\psi\rangle\rightarrow U|\psi\rangle$ if $U^{\dagger}HU=H$. For instance, because each $A_{c}$ and $B^{\mu}_{v}$ commutes, $H$ is invariant under the local (gauge) symmetries defined by $U=A_{c}$ or $U=B^{\mu}_{v}$ for any $A_{c}$ or $B^{\mu}_{v}$. In addition, let $\mathcal{U}^{\mu}_{i}=\prod_{n\in\bar{P}^{\mu}_{i}}\sigma^{y}_{n}$ be the product of $\sigma^{y}_{n}$ for each qubit in the $i$th plane perpendicular to the direction $\mu$, indicated by $\bar{P}^{\mu}_{i}$ (see figure 2a). $1\leq i\leq L$ for each periodic dimension $\mu$, and $0\leq i\leq L$ for each open dimension — see Sections 5 and 6 for the precise meaning of open and periodic dimensions. For any individual $n$th qubit in the plane, $\sigma^{x}_{n}\rightarrow\sigma^{y}_{n}\sigma^{x}_{n}\sigma^{y}_{n}=-\sigma^{x}_{n}$ and $\sigma^{z}_{n}\rightarrow\sigma^{y}_{n}\sigma^{z}_{n}\sigma^{y}_{n}=-\sigma^{z}_{n}$, so $\mathcal{U}^{\mu}_{i}$ flips the sign of $\langle\sigma^{x}_{n}\rangle$ and $\langle\sigma^{z}_{n}\rangle$. However, since any $A_{c}$ is composed of either zero or four spins in any $i$th plane, and any $B^{\mu}_{v}$ is composed of zero, two, or four spins in any $i$th plane, $H$ is left invariant under the transformations defined by $\mathcal{U}^{\mu}_{i}$; these are $d=2$-dimensional Gauge like symmetries. Finally, consider the product $\mathcal{V}$ of each $z$-Pauli operator forming a line parallel to the direction $\mu$, lying on links perpendicular to $\mu$ (see figure 2b). For instance, labeling one corner of the lattice as the origin, these may be the $z$-Pauli operators lying at locations $(i,J,K+\frac{1}{2})$ for all $1\leq i\leq L$ and for some particular $J$ and $K$. Because $\mathcal{V}$ is composed of $z$-Pauli operators, it automatically commutes with each $B^{\mu}_{v}$, and because $\mathcal{V}$ shares zero or two common spin sites with each $A_{c}$, $\mathcal{V}$ commutes with each $A_{c}$ as well, and $H$ is therefore invariant under the transformations defined by $\mathcal{V}$; these are $d=1$-dimensional Gauge like symmetries.

Partition function — The X-Cube model partition function is:

[TABLE]

where $\beta\equiv 1/k_{B}T$ is the inverse temperature. Since each operator in (6) commutes, we may rewrite (7) as:

[TABLE]

Throughout the following sections, when performing a high temperature (small $\beta$) series expansion, we will utilize the following properties for simplifying (8): Using $(A_{c})^{2}=\mathds{1}$, we can rewrite each exponential as

[TABLE]

The first product in (8) can then be written as:

[TABLE]

where we have defined $\mathsf{C}_{a}\equiv\cosh(\beta a)$ and $\mathsf{T}_{a}\equiv\tanh(\beta a)$ (and later $\mathsf{S}_{a}\equiv\sinh(\beta a)$) for the sake of brevity. The righthand product in (10) will contain one linear term for every possible combination of $A_{c}$ operators. For clarity, the first few terms of the product are:

[TABLE]

We thus obtain a series in powers of $\mathsf{T}_{a}$ realizing a high temperature series expansion. Using the linearity of the trace, the following feature will be used to discard significant portions of (7) which contribute no trace. Because each Pauli operator is traceless, and each product of Pauli operators is also a Pauli operator, the only terms in (8) contributing a trace will be those proportional to the identity. These terms will depend on the choice of boundary conditions.

Free energy — Once each partition function is found, the corresponding free energy density is given by:

[TABLE]

For our purposes, the free energy density is important for two reasons: first, we will show that the free energy density of the X-Cube model is independent of our choice of boundary conditions in the thermodynamic limit; second, we will show that this thermodynamic free energy density is completely analytic, indicating the absence of finite temperature phase transitions.

5 Open Boundary Conditions

We first assume our $L\times L\times L$ lattice has open boundary conditions. By this, we mean that our lattice is fully non-periodic and has boundaries. This choice of boundary conditions will lead to the simplest solution for (7), but the result is no less meaningful: we will soon see that in the large system size (thermodynamic) limit, any corrections (arising from our choice of boundary conditions) to an extensive physical quantity calculated via derivatives of the partition function will only appear to order $L^{2}$ and higher.

The $L\times L\times L$ lattice is constructed using $L+1$ horizontal square $L\times L$ lattices, with $(L+1)^{2}$ vertical rungs of $L$ edges each connecting the horizontal lattices, as seen in figure 3a. Each square lattice consists of $2L(L+1)$ edges: $L+1$ rungs of $L$ edges horizontally, and $L+1$ rungs of $L$ edges vertically, as seen in figure 3b. Then, the total number of qubits $N$ is given by $3L^{3}+6L^{2}+3L$.

Using open boundary conditions, we must be careful in how we define the second sum in (6): for vertices on the boundary of the lattice, $B^{\mu}_{v}$ will only be properly defined for $\mu$ perpendicular to the boundary. Therefore, we include only the $(L-1)^{3}$ interior vertices in the sum. We emphasize that this choice will not significantly affect the resulting partition function, and will have no effect whatsoever on the free energy in the thermodynamic limit. We will remark on how the partition function trivially changes for different choices after obtaining our solution.

5.1 Partition Function

5.1.1 High Temperature Series Expansion

To evaluate (8), we begin by expanding the first set of exponentials using (9). We then note that, since each $\sigma^{x}_{n}$ and $\sigma^{z}_{n}$ is traceless, and $\sigma^{x}_{n}\sigma^{z}_{n}=-i\sigma^{y}_{n}$ is traceless, the only terms contributing to the trace will be those which are proportional to $\mathds{1}$. However, no terms of (11) containing $A_{c}$ operators can yield the identity – for any connected section of elementary cubes, the operators $\sigma^{x}_{n}$ corresponding to the boundaries of the section will appear in the product only once. Furthermore, no product of non-identity $A_{c}$ and $B^{\mu}_{v}$ operators will yield the identity: because $(\sigma^{x}_{n})^{2}=(\sigma^{z}_{n})^{2}=\mathds{1}$, any non-identity product of $A_{c}$ operators will simply be a product of $\sigma^{x}_{n}$ operators to the first power, each of which cannot be canceled by a $\sigma^{z}_{n}$ of power zero or one. Equation (7) therefore reduces to:

[TABLE]

To further simplify (13), we apply the same procedure, but we must be careful of the constraint (3). We start by expanding the exponentials in the same manner as (9) and (11):

[TABLE]

where $\mathsf{C}_{b}\equiv\cosh(\beta b)$ and $\mathsf{T}_{b}\equiv\tanh(\beta b)$ (and later $\mathsf{S}_{b}\equiv\sinh(\beta b)$). The trace-contributing terms are those proportional to the identity. In particular, each term of (14) will be proportional to the identity if and only if all three $B^{\mu}_{v}$ are included or excluded at once for each vertex. Each term containing all three $B^{\mu}_{v}$ for a given set of $n$ vertices will carry a factor of $\mathsf{T}_{b}^{3n}$, and there are $(L-1)^{3}\choose n$ configurations of $n$ vertices out of $(L-1)$ total vertices. Equation (14) can therefore be written as:

[TABLE]

where “t.t.” stands for “traceless terms”. Using the binomial theorem, and including the factor of $\mathsf{C}_{b}^{3(L-1)^{3}}$, we can evaluate the trace in (13):

[TABLE]

Note that the second product after the equality can be expanded by choosing, for each of $(L-1)^{3}$ factors of $[1+\mathsf{T}_{b}^{3}]$, either a factor of 1 or of $\mathsf{T}_{b}^{3}$. This term altogether can therefore be interpreted to represent the sum of all possible choices for including or excluding all three $B^{\mu}_{v}$ at each of $(L-1)^{3}$ vertices. We will therefore sometimes skip the binomial theorem altogether, when the meaning of such a term is clear.

The trace of the identity is $2^{N}=2^{3L^{3}+6L^{2}+3L}$, the dimension of the total state space. By combining the $\cosh$ and $\tanh$ terms, (13) is finally given by:

[TABLE]

Alternatively, we may use the following to rewrite the latter product in terms of exponentials:

[TABLE]

Rewriting in terms of exponentials, and combining powers of two, (17) is given in the form:

[TABLE]

The partition function found by a high temperature series expansion is a regular function. Thus, for all finite $L$ and $\beta$, the high temperature series expansion that we invoked leads to a convergent answer.

Additionally, we may now consider alternative definitions of the vertex sum in (6). In addition to the aforementioned treatment of considering only the $(L-1)^{3}$ interior vertices, there are two other sensible methods of performing this sum. First, we may suppose that for each boundary vertex, we only include $B^{\mu}_{v}$ operators which are already properly defined – for instance, a vertex on a boundary perpendicular to the $x$ direction only has a properly defined $B^{x}_{v}$ operator, and no vertex on any edge or corner of a boundary has any properly defined $B^{\mu}_{v}$ operators. In this case, (17) is exactly identical to its current form: these newly introduced boundary $B^{\mu}_{v}$ operators are traceless, cannot be fully canceled by any product of operators, and have no such constraint (3). Therefore, any term in (14) containing a boundary $B^{\mu}_{v}$ operator has trace zero.

Alternatively, we could additionally include partial $B^{\mu}_{v}$ operators at the boundary – that is, if not all four qubits necessary to define a given $B^{\mu}_{v}$ are present, we simply include in every $B^{\mu}_{v}$ only the $z$-Pauli operators which are defined. For instance, in the extreme case, the $B^{\mu}_{v}$ operators corresponding to each of the eight corner vertices consist of only two $z$-Pauli operators each. In this case, each boundary vertex has three $B^{\mu}_{v}$ operators satisfying (3). This choice makes (13) particularly difficult to solve, as it introduces $3(L+1)$ additional planar constraints. Such planar constraints will be discussed in following sections, and we will see in section 7 that such a partition function is very difficult to solve. However, we will still find that such a choice leads to exactly the same free energy in the thermodynamic limit.

Given (17), the corresponding free energy density is:

[TABLE]

Note that in the thermodynamic, $L\rightarrow\infty$, limit, our particular choice of definition of the model at the boundaries of the system does not affect the free energy density, indicating the bulk nature of this physical quantity. Also note that (20) is a regular function for all finite $\beta$.

5.1.2 Bond Algebraic Duality

In addition to the brute-force high temperature series expansion, (17) can be obtained using a far simpler approach, the bond-algebraic method. Since each operator in (6) commutes, we may map each operator to the classical spins $r_{m}$ and $s^{n}_{j}$, $j=1,2$, as follows (see figures 4 and 5):

[TABLE]

The mapping (21) preserves the bond algebra of the system – namely, it respects the conditions (4). We therefore expect this classical mapping to preserve the spectrum of (6), with degeneracies that can only differ by a global power of two. By extension, (21) will maintain the form of (7) up to a power of two. By explicitly writing the sum over $\mu$ in (7) and applying the mapping (21), we obtain the partition function:

[TABLE]

This expression is easily summed directly to yield

[TABLE]

The value of ${\bar{N}}$ is then determined to be $12L^{2}-3L+2$ by taking the infinite temperature ($\beta\rightarrow 0$) limit and demanding that $\mathcal{Z}_{\text{Open}}$ is equal to the dimension of the total state space in this limit. Equation (23) is then quickly verified to agree with (17) and therefore its free energy with (20). This factor of $2^{\bar{N}}$ indicates that the two dual models have the same spectra, but degeneracies at each energy level that differ by a global factor of $2^{\bar{N}}$.

We also note that the mapping (21) could also have been accomplished by mapping each $B^{\mu}_{v}$ to the classical spins (see figure 5):

[TABLE]

These mappings and the corresponding partition function calculations show that (6) under open boundary conditions is dual to the classical Hamiltonian of a single open Ising chain of length $L^{3}+1$ with bond variables $r_{m}$ and $(L-1)^{3}$ two-site Ising chains under a magnetic field (or equivalently, $(L-1)^{3}$ periodic three-site Ising chains).

5.2 Cube-Star Correlation Functions

Correlation functions can also be easily determined under open boundary conditions, using either brute-force trace calculations or bond-algebraic mappings. For $\cal{M}$ cubic operators $A_{c_{1}},A_{c_{2}},\ldots A_{c_{\cal{M}}}$ and $\cal{N}$ vertex operators $B^{\mu_{1}}_{v_{1}},B^{\mu_{2}}_{v_{2}},\ldots B^{\mu_{\cal{N}}}_{v_{\cal{N}}}$, we wish to calculate:

[TABLE]

which is the most general nonvanishing correlation function one can compute: any spin product not expressible as a product of $A_{c}$ and $B^{\mu}_{v}$ operators will have a spectrum of $\pm 1$ with equal Boltzmann weights on each eigenvalue, and will therefore have an expectation value of zero.

In this calculation, we will assume without loss of generality that each of the $\cal{M}$ operators $A_{c_{m}}$ and each of the $\cal{N}$ operators $B^{\mu_{n}}_{v_{n}}$ are distinct. Because $(A_{c_{m}})^{2}=(B^{\mu_{n}}_{v_{n}})^{2}=\mathds{1}$, and because all operators commute, any instances of indistinct operators can be discarded. In addition, we may assume that each of the $\cal{N}$ $B^{\mu_{n}}_{v_{n}}$ operators correspond to distinct vertices. If $v_{n}=v_{n^{\prime}}$ for $n\neq n^{\prime}$, then we may use (4) to reduce $B^{\mu_{n}}_{v_{n}}B^{\mu_{n^{\prime}}}_{v_{n^{\prime}}}$ to a single operator. If $v_{n}=v_{n^{\prime}}=v_{n^{\prime\prime}}$ for $n\neq n^{\prime}\neq n^{\prime\prime}$, then we may use (3) to eliminate the three operators altogether. Finally, for the sake of brevity, we define $\mathbb{A}\equiv A_{c_{1}}\ldots A_{c_{\cal{M}}}$ and $\mathbb{B}\equiv B^{\mu_{1}}_{v_{1}}\ldots B^{\mu_{\cal{N}}}_{v_{\cal{N}}}$.

In order to calculate (25), we start by using (10) and (14):

[TABLE]

where we have preemptively canceled the $\mathsf{C}_{a}$ and $\mathsf{C}_{b}$ powers in the numerator and denominator. As before, the only terms inside each trace are those proportional to the identity. Starting with the $A_{c}$ operators, because no nontrivial product of $A_{c}$ operators can yield the identity, the only term of (11) which can yield the identity when multiplied by $\mathbb{A}$ is the term containing $\mathbb{A}$ itself, which carries with it a factor of $\mathsf{T}_{a}^{\cal{M}}$. In the denominator, the only term of (11) proportional to the identity is $\mathds{1}$ itself. Equation (26) therefore reduces to:

[TABLE]

For each $B^{\mu_{n}}_{v_{n}}$, let $\mu_{n}^{\prime}$ and $\mu_{n}^{\prime\prime}$ denote the two cardinal directions of $\{x,y,z\}$ not equal to $\mu_{n}$. Each term in the numerator of (27) with nonzero trace must cancel each $B^{\mu_{n}}_{v_{n}}$, which can only be canceled by $B^{\mu_{n}}_{v_{n}}$ itself or by $B^{\mu_{n}^{\prime}}_{v_{n}}B^{\mu_{n}^{\prime\prime}}_{v_{n}}$. Each term with nonvanishing trace in the numerator of (27) equates to picking either $B^{\mu_{n}}_{v_{n}}$ with a factor of $\mathsf{T}_{b}$ or $B^{\mu_{n}^{\prime}}_{v_{n}}B^{\mu_{n}^{\prime\prime}}_{v_{n}}$ with a factor of $\mathsf{T}_{b}^{2}$ for each $1\leq n\leq\cal{N}$. In addition, just as in calculating (15), we may pick any number of vertices $v$ out of the $(L-1)^{3}-\cal{N}$ vertices not included in $\{v_{1},\ldots v_{\cal{N}}\}$ and include $B^{x}_{v}B^{y}_{v}B^{z}_{v}$. Therefore, (27) is written as:

[TABLE]

Evaluating each summation using the binomial theorem and simplifying yields the solution:

[TABLE]

From (29), we immediately see that each $A_{c}$ and each $B^{\mu}_{v}$ on distinct vertices are entirely uncorrelated. That is, the expectation value (25) is simply given by the product of individual expectation values:

[TABLE]

This is independent of both the size of the system and the relative placements of cubes and vertices. The expectation values of $A_{c}$ operators for two adjacent cubes are just as uncorrelated as those of cubes very far apart, and the same is true of vertex operators $B^{\mu}_{v}$ for distinct vertices. Under open boundary conditions, neither the quantum mechanics nor the thermodynamics of (6) “know” about the geometry of the lattice: each cube and each vertex form their own isolated system.

The same solution is easily obtained using the mapping (21). First, we map each $B^{\mu}_{v}$ to some $s^{n}_{1}$: although we previously mapped $B^{x}_{v}$ in particular to each $s^{n}_{1}$, we may instead choose for each vertex individually how we map bijectively from $\{B^{x}_{v},B^{y}_{v},B^{z}_{v}\}$ to $\{s^{n}_{1},s^{n}_{2},s^{n}_{1}s^{n}_{2}\}$. We obtain the same result with any such mapping, but for ease of calculations, we choose our mapping such that each $B^{\mu}_{v}$ included in (25) is mapped to an $s^{n}_{1}$. In addition, we will without loss of generality enumerate our mapping such that $A_{c_{m}}\rightarrow r_{m}$ for each $1\leq m\leq\cal{M}$ and $B^{\mu_{n}}_{v_{n}}\rightarrow s^{n}_{1}$ for each $1\leq n\leq\cal{N}$. Defining $\mathbb{R}\equiv r_{1}\ldots r_{\cal{M}}$ and $\mathbb{S}\equiv s^{1}_{1}\ldots s^{\cal{N}}_{1}$ for brevity, (25) is then mapped to:

[TABLE]

This expression is easily evaluated by directly summing the first $\cal{M}$ variables $r_{m}$ and the first $\cal{N}$ variables $s^{n}_{1}$ and $s^{n}_{2}$. In the numerator, each factor of $r_{m}e^{\beta ar_{m}}$, when summed over $r_{m}=\pm 1$, yields a factor of $e^{\beta a}-e^{-\beta a}$. Similarly, each factor of $s^{n}_{1}\exp[\beta b(s^{n}_{1}+s^{n}_{2}+s^{n}_{1}s^{n}_{2})]$, when summed over $s^{n}_{j}=\pm 1$, yields a factor of $e^{3\beta b}-e^{-\beta b}$. In the denominator, each factor of $e^{\beta ar_{m}}$ sums to a factor of $e^{\beta a}+e^{-\beta a}$, and each factor of $\exp[\beta b(s^{n}_{1}+s^{n}_{2}+s^{n}_{1}s^{n}_{2})]$ sums to a factor of $e^{3\beta b}+3e^{-\beta b}$. Canceling the remaining sums over $\{r_{m},s^{n}_{j}\}$ for $m>\cal{M}$ and $n>\cal{N}$, which are identical in the numerator and denominator, (31) therefore evaluates to:

[TABLE]

6 Cylindrical Boundary Conditions

Next, we consider a cubic lattice periodic along its $y$ and $z$ axes, but open along its $x$ axis. Under this cylindrical topology, a plane of simple cubes extending in the $y$ and $z$ directions (that is, perpendicular to the $x$ direction) will contain cubes sharing common edges across the periodic boundaries. In this plane of cubes, all edges parallel to the $x$ direction will be shared by four common cubes, while all edges parallel to the $y$ and $z$ directions will be shared by two common cubes — see figure 6a. The $A_{c}$ operators corresponding to these cubes will therefore cancel along a plane $P^{x}_{i}$ of simple cubes perpendicular to the $x$ axis, forming $L$ additional constraints of the form:

[TABLE]

Similarly, any set of $B^{x}_{v}$ operators corresponding to a plane of vertices perpendicular to the $x$ axis will share each of its four qubits with exactly one other $B^{x}_{v}$ due to the periodicity in the $y$ and $z$ directions. We therefore have an additional $L-1$ constraints given by:

[TABLE]

where $\bar{P}^{x}_{i}$ here refers to a plane of vertices perpendicular to the $x$ axis — see figure 6b.

Under this topology, the number $N$ of qubits on the lattice is given by $3L^{3}+2L^{2}$. This is obtained by starting with the $3L^{3}$ qubits in the periodic lattice, and closing the lattice in the $x$ direction with an additional $2L^{2}$ edges (see figure 7). In addition, we again must be careful in defining the second sum in (6). Whereas each of the $L^{2}$ vertices on an $x$-plane interior to the lattice have all three $B^{\mu}_{v}$ operators properly defined, vertices on the two exterior $x$-planes have no properly defined $B^{y}_{v}$ or $B^{z}_{v}$ operators. Therefore, we include only the $L^{2}(L-1)$ interior vertices in the sum. Other sensible definitions, as well as their implications, are the same as those discussed in the end of section 5.1.

6.1 Partition Function

6.1.1 High Temperature Series Expansion

First, we evaluate (7) directly, using (10) and (14). Starting with the $A_{c}$ operators, the terms proportional to the identity in (11) are those satisfying (33) for some number $n$ of $x$-planes with $0\leq n\leq L$. Out of $L$ total $x$-planes, there are $L\choose n$ configurations of $n$ planes, with each plane contributing a factor of $\mathsf{T}_{a}^{L^{2}}$. We may therefore rewrite (11) under cylindrical boundary conditions as:

[TABLE]

Using the binomial theorem and combining with the factor of $\mathsf{C}_{a}^{L^{3}}$, (8) is reduced to:

[TABLE]

The terms contributing to a nonvanishing trace in (36) are those satisfying (3) and (34). However, we must be careful when using (34): due to (4), we may replace $B^{x}_{v}$ in any product forming an $x$-plane with $B^{y}_{v}B^{z}_{v}$ and still yield the identity. These terms will have an additional factor of $\mathsf{T}_{b}$, since two operators are included for the vertex $v$ instead of one. We therefore sum (14) for cylindrical boundary conditions as follows: given $n$ $x$-planes of vertices, with $0\leq n\leq L-1$, there are $L-1\choose n$ configurations of these planes. A priori, each of these planes is composed of $B^{x}_{v}$ operators alone, and contributes a factor of $\mathsf{T}_{b}^{L^{2}}$. However, for each of the $nL^{2}$ vertices on the $n$ planes, we may replace $m$ of these $B^{x}_{v}$ operators with $B^{y}_{v}B^{z}_{v}$ in $nL^{2}\choose m$ different configurations, with each vertex replacement contributing an additional factor of $\mathsf{T}_{b}$. Finally, out of the $L^{2}(L-1)-nL^{2}$ vertices not lying on any $x$-plane, we may include any $\ell$ products of the form (3) in $L^{2}(L-1)-nL^{2}\choose\ell$ different configurations, each contributing a factor of $\mathsf{T}_{b}^{3}$. (14) can therefore be written under cylindrical boundary conditions as:

[TABLE]

where $Q\equiv L^{2}(L-1)-nL^{2}$. Equation (37) is easily summed using three applications of the binomial theorem. First, we evaluate the sums over $m$ and $\ell$ at once, since they share no dependencies. Equation (36) becomes:

[TABLE]

Isolating the $n$ dependencies, applying the binomial theorem once again, and combining $\mathsf{C}_{b}$ and $\mathsf{T}_{b}$ terms, we obtain our solution:

[TABLE]

where the righthand trace is $2^{N}=2^{3L^{3}+2L^{2}}$, the dimension of the total state space. In addition, the $\mathsf{C}_{b}$ and $\mathsf{S}_{b}$ terms may be rewritten as exponentials using (18) and:

[TABLE]

Combining powers of two, the partition function is finally given by:

[TABLE]

Note that this expression appears very similar to that of (17) in the high-temperature ($\beta\rightarrow 0$) limit. In fact, by dividing (41) by (17), we obtain:

[TABLE]

This function is perfectly regular for all $\beta\geq 0$. This illustrates, that similar to the case of open boundary conditions, our high temperature series expansion converges for all $\beta$.

The free energy density corresponding to (39) is:

[TABLE]

which could alternatively be written as:

[TABLE]

Written this way, we see that (43) rapidly converges with (20) as $L$ increases, as the final two terms of (44) quickly approach zero. This is yet another indication that our choice of boundary conditions should not affect the thermodynamics of the model at large system sizes.

6.1.2 Bond Algebraic Duality

Once again, the most economical way to evaluate (7) is via a duality mapping each $A_{c}$ and $B^{\mu}_{v}$ to classical spins preserving the original bond algebra. This is achieved through the following mapping:

[TABLE]

where the upper $i$ index in $r^{i}_{m}$ is constant along each of the $L$ cubic $x$-planes, and the upper $j$ index in $s^{j}_{n}$ and $t^{j}_{n}$ is constant along each of the $L-1$ vertex $x$-planes. Effectively, each $x$-plane of $A_{c}$ operators is mapped to a periodic Ising chain of length $L^{2}$ (see figure 8), while each $x$-plane of $B^{\mu}_{v}$ operators is mapped to what might be thought of as an Ising-gauge chain (see figure 9). Under this mapping, (33) and (34) are mapped to:

[TABLE]

where $L^{2}+1\equiv 1$ in the above products. Similarly, the constraint (3) is mapped to:

[TABLE]

We therefore see that this mapping indeed preserves the bond algebra of the operators in (6).

The partition function is then given by:

[TABLE]

where the first and third terms inside the parenthesis represent the Ising-gauge chain (see figure 9). This expression is most easily summed by first summing the $r^{i}_{m}$ spins, then the $t^{j}_{n}$ spins, and finally the $s^{j}_{n}$ spins. The $r^{i}_{m}$ spins simply form $L$ independent periodic Ising chains each of length $L^{2}$:

[TABLE]

Next, because each $t^{j}_{n}$ is uncoupled, they can be summed directly:

[TABLE]

The remaining sum is evaluated using the transfer matrix method: let $T$ be a $2\times 2$ matrix given by:

[TABLE]

Then (50) is reduced to the form:

[TABLE]

The eigenvalues of $T$ are given by $e^{3\beta b}+3e^{-\beta b}$ and $e^{3\beta b}-e^{-\beta b}$, reducing (52) to:

[TABLE]

Once again, the value of $\tilde{N}$ is determined as $4L^{2}$ by taking the infinite temperature ($\beta\rightarrow 0$) limit and demanding that $\mathcal{Z}_{\text{Cylindrical}}$ is given by $2^{N}$ in this limit, where $N=3L^{3}+2L^{2}$. Combining powers of two, we see that the solution given by the mapping (45) exactly matches that of (41), showing that (6) under cylindrical boundary conditions is dual to the classical Hamiltonian of $L$ periodic Ising chains of length $L^{2}$ and $L-1$ periodic Ising-gauge chains. The factor of $2^{\tilde{N}}$ again indicates that the dual models have the same energy levels, with degeneracies that differ only by a global factor of $2^{\tilde{N}}$.

6.2 Cube-Star Correlation Functions

To compute (25) under this topology, we must be particularly careful of $A_{c}$ operators lying on the same cubic $x$-plane and $B^{\mu}_{v}$ operators lying on the same vertex $x$-plane. Labeling the $L$ cubic $x$-planes 1 through $L$ and the $L-1$ vertex $x$-planes 1 through $L-1$, let ${\cal{M}}_{i}$ for $1\leq i\leq L$ denote the number of $A_{c}$ operators lying in the $i$th plane that are included in the expectation value (25), and let ${\cal{N}}_{j}$ for $1\leq j\leq L-1$ denote the number of $B^{\mu}_{v}$ operators lying in the $j$th plane that are included in (25). Furthermore, let $A_{i;m}$ denote the $m$th included operator in the $i$th cubic plane, and let $B^{\mu_{j;n}}_{j;n}$ denote the $n$th included operator in the $j$th vertex plane. Equation (25) is then given by:

[TABLE]

Once again, (54) is the most general nonvanishing correlation function: any spin product not composed of the $A_{c}$ and $B^{\mu}_{v}$ operators composing (6) will necessarily have an expectation value of zero. In addition, we again assume in our calculation of (54) that each $A_{i;m}$ and $B^{\mu_{j;n}}_{j;n}$ is distinct, and that each $B^{\mu_{j;n}}_{j;n}$ lies on distinct vertices: any product of $A_{c}$ and $B^{\mu}_{v}$ operators can be reduced to this form using $(A_{c})^{2}=(B^{\mu}_{v})^{2}=\mathds{1}$ and (4). Just as in the case of open boundary conditions, we calculate (54) using an expansion of the form (26). Under cylindrical boundary conditions, each cubic planar product $\prod_{m=1}^{{\cal{M}}_{i}}A_{i;m}$ can be canceled by one of two products: either $\prod_{m=1}^{{\cal{M}}_{i}}A_{i;m}$ itself, or the remaining $L^{2}-{\cal{M}}_{i}$ cubic operators in the $i$th plane not included in $\prod_{m=1}^{{\cal{M}}_{i}}A_{i;m}$. Therefore, each term of (11) contributing nonvanishing trace is given by choosing for each plane either the product $\prod_{m=1}^{{\cal{M}}_{i}}A_{i;m}$ with a factor of $\mathsf{T}_{a}^{{\cal{M}}_{i}}$ or the entire plane except $\prod_{m=1}^{{\cal{M}}_{i}}A_{i;m}$ with a factor of $\mathsf{T}_{a}^{L^{2}-{\cal{M}}_{i}}$. The $A_{c}$ operators in the denominator are handled as in (36), reducing (54) to:

[TABLE]

Each product $\prod_{n=1}^{{\cal{N}}_{j}}B^{\mu_{j;n}}_{j;n}$ can be canceled in one of two ways. One possibility is that they are canceled pointwise, where each $B^{\mu}_{v}$ is matched by an identical $B^{\mu}_{v}$ or by (4). A priori, this incurs a factor of $\mathsf{T}_{b}^{{\cal{N}}_{j}}$, with an additional factor of $\mathsf{T}_{b}^{n}$ for each choice of $n$ operators canceled using (4). We may also include $m$ factors of $\mathsf{T}_{b}^{3}$ for $0\leq m\leq L^{2}-{\cal{N}}_{j}$ by applying (3) to $m$ of the $L^{2}-{\cal{N}}_{j}$ vertices not included in $\prod_{n=1}^{{\cal{N}}_{j}}B^{\mu_{j;n}}_{j;n}$. The other method of canceling $\prod_{n=1}^{{\cal{N}}_{j}}B^{\mu_{j;n}}_{j;n}$ is to include the remaining factors necessary to satisfy (34). To do this, we must account for the particular directions $\mu_{j;n}$ involved: let ${\cal{X}}_{j}$ be the number of $B^{x}_{v}$ included in the $j$th $x$-plane. For each such operator, we may either simply include each such $B^{x}_{j;n}$ in the product, or we may include (3) with a factor of $\mathsf{T}_{b}^{3}$. For each additional ${\cal{N}}_{j}-{\cal{X}}_{j}$ operators included in $\prod_{n=1}^{{\cal{N}}_{j}}B^{\mu_{j;n}}_{j;n}$, we may either include $B^{\mu_{j;n}^{\prime}}_{j;n}$ for $\mu_{j;n}^{\prime}\neq x$, or we may include $B^{\mu_{j;n}}_{j;n}B^{x}_{j;n}$, each converting $B^{\mu_{j;n}}_{j;n}$ to the required $B^{x}_{j;n}$ to satisfy (34). For the remaining $L^{2}-{\cal{N}}_{j}$ vertices in the $j$th plane, we may either include $B^{x}_{v}$ or $B^{y}_{v}B^{z}_{v}$ to satisfy (34), contributing factors of $\mathsf{T}_{b}$ and $\mathsf{T}_{b}^{2}$ respectively. The numerator trace of (55) is therefore given by:

[TABLE]

Using the binomial theorem and combining with the results of (39), (55) can be simplified to:

[TABLE]

First, note that cubes and stars are once again uncorrelated: $\langle A_{c}B^{\mu}_{v}\rangle=\langle A_{c}\rangle\langle B^{\mu}_{v}\rangle$ for each $c,\mu,v$. This is a general feature of the X-Cube model, independent of boundary conditions, as the partition function will always factor into the product of a partition function for the $A_{c}$ operators and a partition function for the $B^{\mu}_{v}$ operators. Additionally, note that each cubic $x$-plane and each vertex $x$-plane are entirely uncorrelated: for any cubes $c$ and $d$ on distinct cubic $x$-planes, $\langle A_{c}A_{d}\rangle=\langle A_{c}\rangle\langle A_{d}\rangle$, and the same holds for two $B^{\mu}_{v}$ operators on different vertex $x$-planes. Each $x$-plane of cubes or vertices forms its own isolated system, as suggested by the duality (45). Finally, note that cubes and stars within the same $x$-plane are very weakly correlated in large system sizes: (57) is trivially verified to agree with (29) in the thermodynamic limit.

7 Periodic Boundary Conditions

Care must be exercised in evaluating (7) when using fully periodic boundary conditions. Whereas the constraints (35) were entirely independent of each other, the corresponding constraints under boundary conditions periodic in all three cardinal directions are not so simple.

In particular, let $P^{\mu}_{i}$ denote a set of simple cubes $c$ forming an $L\times L$ plane perpendicular to the direction $\mu$. Then, we have the $3L$ constraints:

[TABLE]

Note that, unlike the constraints (33), these constraints do not all function independently from one another: for instance, two perpendicular planes of simple cubes share a common line of $L$ cubes at their intersection, and three perpendicular planes additionally share one common cube at their shared intersection. We will find as a consequence that properly counting the possible products in (11) with nonzero trace becomes far less trivial than in the previous two cases, and that the binomial theorem alone will prove insufficient for obtaining a closed-form partition function under this topology.

Similarly, let $\bar{P}^{\mu}_{i}$ denote any set of vertices $v$ forming a plane perpendicular to the direction $\mu$. We then also have the $3L$ constraints:

[TABLE]

These constraints will pose different difficulties than the cubic constraints (58): while two perpendicular planes of $B^{\mu}_{v}$ operators will not share any common operators, the additional constraints (3) will lead to many more possible trace-containing factors, further complicating the trace calculation. These difficulties, along with those from the constraints (58), will render our previous strategies insufficient in finding a closed form solution to (7).

Finally, we note that fully periodic boundary conditions leads to a very natural choice of the second sum in (6): we simply sum over all $L^{3}$ vertices in the lattice, as the elimination of hard boundaries allow us to properly define each $B^{\mu}_{v}$ for all vertices. Additionally, the counting of spins is trivial: each of the $L^{3}$ vertices is simply associated with the three edges located directly in the positive $x$, $y$, and $z$ directions from the vertex, giving the total number of qubits $N$ as simply $3L^{3}$.

7.1 Partition Function

We might attempt to simplify (8) using (10) as follows: the terms in (11) proportional to the identity are those for which (58) can be used to eliminate products of $A_{c}$ operators – that is, each product corresponding to some set of planes of cubes will be proportional to the identity. A priori, if a given linear term in (11) contains the products of $n$ $x$-planes, $m$ $y$-planes, and $\ell$ $z$-planes, it will carry a factor of $\mathsf{T}_{a}^{(n+m+\ell)L^{2}}$ corresponding to the $L^{2}$ cubes in each of the $n+m+\ell$ planes. However, we must be cautious when dealing with perpendicular planes: because each linear term in (11) may contain at most one of each $A_{c}$ operator, and two perpendicular planes share a common intersection of $L$ cubes, two perpendicular planes cannot each individually contribute $L^{2}$ factors of $\mathsf{T}_{a}$. Instead, using $A_{c}^{2}=\mathds{1}$, we include no operators corresponding to the linear intersection of two given planes. Because each of the $nm+n\ell+m\ell$ intersections of planes are double-counted initially, we remove $2L$ factors of $\mathsf{T}_{a}$ at each intersection. Finally, at each of the $nm\ell$ intersections of three perpendicular planes, we must add back in the $A_{c}$ operator and its corresponding factor of $\mathsf{T}_{a}$: using $A_{c}=A_{c}^{3}$, only one $A_{c}$ is necessary to satisfy all three constraints (58) for three perpendicular planes. Since these intersections are initially counted three times and then removed six times, we must add four factors of $\mathsf{T}_{a}$ back to the product for these $nm\ell$ planes. Therefore, each product corresponding to $n$ $x$-planes, $m$ $y$-planes, and $\ell$ $z$-planes carries a factor of:

[TABLE]

Unfortunately, more care is needed than simply summing each of $n$, $m$, and $\ell$ from [math] to $L$ with the appropriate binomial coefficients: because the constraints (58) are not all independent of each other, a given trace-containing product of $A_{c}$ operators can be interpreted as one of several possible planar configurations. As a simple example, the product of all $L^{3}$ cubes can be interpreted as $L$ $x$-planes, or $L$ $y$-planes, or $L$ $z$-planes, or the intersection of all $3L$ planes. Indeed, while this particular example is trivial to see, there exist many more nontrivial examples of a given operator product having multiple interpretations in terms of the constraints (58). As a result, the simple counting argument employed in Sections 5.1 and 6.1 will lead to dramatic overcounting of the terms proportional to the identity.

These same sort of difficulties are not present immediately in counting the trace-containing terms in (14), but the resulting computation is no less difficult. We might attempt to proceed as follows: a given product corresponding to $n$ $x$-planes, $m$ $y$-planes, and $\ell$ $z$-planes of vertices will a priori contribute a factor of $\mathsf{T}_{b}^{(n+m+\ell)L^{2}}$. Then, each vertex in the lattice lies on exactly zero, one, two, or three planes. Each vertex lying on zero planes could contribute either zero or three additional factors of $\mathsf{T}_{b}$ by choosing to include factors of the form (3). Each vertex lying on exactly one plane may contribute an additional factor of $\mathsf{T}_{b}$ by substituting $B^{\mu}_{v}$ contained in (59) with $B^{\mu^{\prime}}_{v}B^{\mu^{\prime\prime}}_{v}$ using (4). Each vertex lying on exactly two perpendicular planes may contribute one less factor of $\mathsf{T}_{b}$ using (4) in the same manner. Finally, each vertex lying on three perpendicular planes may contribute three fewer factors of $\mathsf{T}_{b}$ using (3). The problem of counting the trace-containing terms of (14) is then simply reduced to properly counting the number of vertices lying on each of zero, one, two, or three planes given $n$ $x$-planes, $m$ $y$-planes, and $\ell$ $z$-planes. Unfortunately, the constraints (3) introduce the same overcounting problems as before: as a trivial example, the product of all $3L^{3}$ $B^{\mu}_{v}$ operators can be thought of as either the product of all $3L$ planes and zero vertices or of zero planes and of all $L^{3}$ vertices.

7.2 Absence of Phase Transitions

In place of a closed-form partition function, we will instead give a more general argument against the presence of finite temperature phase transitions in the model. Starting from (11), we gather all terms proportional to the identity and denote their prefactors by the function $\mathcal{T}_{a}$:

[TABLE]

$\mathcal{T}_{a}$ is a finite-order power series with nonnegative integer coefficients. Each coefficient $a_{n}$ represents the number of unique products of $n$ $A_{c}$ operators yielding the identity. We do the same for (14), denoting the prefactors of the trace-containing terms by $\mathcal{T}_{b}$:

[TABLE]

As usual, the traceless terms may not combine to yield any terms containing a trace. Using the functions $\mathcal{T}_{a}$ and $\mathcal{T}_{b}$, the fully periodic partition function is given by:

[TABLE]

In order to investigate the possibility of a phase transition, we wish to understand how $\mathcal{T}_{a}$ and $\mathcal{T}_{b}$ behave in the thermodynamic ($L\rightarrow\infty$) limit. We start with $\mathcal{T}_{a}$ by noting that each $a_{n}$ must be bounded above by $2^{3L}$. Because each trace-containing term of (11) is reduced to the identity using the constraints (58), each trace-containing term may be found using the $2^{3L}$ possible products of these $3L$ constraints. Although we know that each of these $2^{3L}$ products will not yield a unique trace-containing term, there can certainly be no greater than $2^{3L}$ trace-containing terms in (11), and therefore no greater than $2^{3L}$ trace-containing terms for each power of $\mathsf{T}_{a}^{n}$. Additionally, we note that $a_{n}=0$ for $1\leq n<L^{2}$, as there are no nontrivial trace-containing terms outside the order $L^{2}$ terms arising from (58). We therefore find:

[TABLE]

We may now safely take the limit of $L\rightarrow\infty$: as $\mathsf{T}_{a}<1$ for all nonzero temperatures, the above sum converges to a finite value while $\mathsf{T}_{a}^{L^{2}}2^{3L}$ goes to zero. We therefore find that $\mathcal{T}_{a}=1$ in the thermodynamic limit. This result is exactly as expected: it is simply a restatement of the idea that the partition function of (6) should be asymptotically the same in the thermodynamic limit for any choice of boundary conditions.

Unfortunately, the same strategy does not follow as easily for $\mathcal{T}_{b}$. Because we expect the free energy to be asymptotically the same as (20), we expect $\mathcal{T}_{b}$ to go asymptotically as:

[TABLE]

Note the exponent of $L^{3}$ in place of $(L-1)^{3}$, arising from our definition of the second sum in (6) under periodic boundary conditions. Of course, this change has no effect on the free energy density in the thermodynamic limit. From this, we do not expect the coefficients $b_{n}$ in $\mathcal{T}_{b}$ to be bounded so easily as in (64).

Instead, we examine the differences in the free energy density to that obtained by presupposing a large $L$ $\mathcal{T}_{b}$ limit of the form (65). That is, we wish to examine:

[TABLE]

In order to bound $\mathcal{T}_{b}$, we perform the calculation described in section 7.1 without regard for possible overcounting. That is, we count every possible combination of planar and vertex products proportional to the identity without worrying about potentially multi-counting a particular configuration in order to obtain an upper bound. This is done as follows: given a product of $n$ $x$-planes, $m$ $y$-planes, and $\ell$ $z$-planes, there exist $nm\ell$ vertices lying on all three planes, $(nm+n\ell+m\ell)L-3nm\ell$ vertices lying on exactly two planes, $(n+m+\ell)L^{2}-2(nm+n\ell+m\ell)L+3nm\ell$ vertices lying on exactly one plane, and $L^{3}-(n+m+\ell)L^{2}+(nm+n\ell+m\ell)L-nm\ell$ vertices lying on no planes. Each vertex lying on zero or three planes may contribute either a factor of $1$ or $\mathsf{T}_{b}^{3}$, and each vertex lying on exactly one or two planes may contribute either a factor of $\mathsf{T}_{b}$ or $\mathsf{T}_{b}^{2}$. We therefore find that $\mathcal{T}_{b}$ is bounded above by:

[TABLE]

The term in parenthesis is less than or equal to one for all temperatures, and its exponent is non-negative. Additionally, each binomial is strictly less than $2^{L}$. We therefore have:

[TABLE]

We also have that $\mathcal{T}_{b}$ is strictly larger than $[1+\mathsf{T}_{b}^{3}]^{L^{3}}$, as $\mathcal{T}_{b}$ contains all zero plane trace-containing terms. Dividing by $[1+\mathsf{T}_{b}^{3}]^{L^{3}}$ and raising to the power of $1/L^{3}$, we find:

[TABLE]

In the thermodynamic limit, we see that the difference in free energies (66) goes to zero. Then, in that limit, the free energy density is simply given by:

[TABLE]

As expected, this is the same thermodynamic limit as in (20) and (44), proving that the order $L^{2}$ constraints (58) and (59) do not matter in the thermodynamic limit, and that (6) under periodic boundary conditions cannot display a finite temperature phase transition.

8 Dynamics of the X-Cube Model at Finite Temperature

In addition to equilibrium correlation functions, we also wish to estimate the time dependent autocorrelation function of each $A_{c}$ and $B^{\mu}_{v}$. The autocorrelation function of a generic operator $X$ as a function of time $t$ is given by:

[TABLE]

where $\rho$ is the density matrix of the system, and $U$ is a unitary time-evolution operator. A priori, since both $A_{c}$ and $B^{\mu}_{v}$ commute with $H$, neither will evolve in time under the Heisenberg picture. To allow for thermalization, we therefore imagine perturbing (6) into a thermal Hamiltonian with an infinitesimal applied field given by:

[TABLE]

where $\lambda>0$ and $\gamma>0$ determine field strengths, and the sum over $n$ includes all edges.

Equation (71) is most easily estimated using the duality map defined by (21): explicitly, it defines an isometry realized by a unitary transformation $\mathcal{U}$, mapping each operator $X$ corresponding to the X-Cube model to the operator $X_{\sf d}$ in the dual representation via:

[TABLE]

While the dual representations given by (21) and (24) are written in terms of classical variables, the same bond algebra is achieved if each classical variable is thought of as a $\sigma^{z}$ operator. This representation of the duality allows us to estimate (71) by calculating the corresponding autocorrelations of simple Ising models using Glauber dynamics [75]:

[TABLE]

where ${\rho}_{\sf d}$ and ${U}_{\sf d}$ are the corresponding Ising density matrix and time-evolution operator. We may therefore find the time evolution of each $A_{c}$ and $B^{\mu}_{v}$ from the time evolution of their classical Ising duals. Although the following results are derived explicitly assuming open boundary conditions, corrections to these results due to boundary conditions such as (58) or (59) are of order $L^{2}$. We therefore expect that these results will accurately describe the bulk material of the system to very high order in the thermodynamic limit of very large $L$.

The dynamics of Ising-type models are investigated thoroughly in [75], with important results cited here. The most important assumption of Glauber dynamics is that of detailed balance: given a system of classical spins $\{s_{i}\}$ and a corresponding Hamiltonian $H$, the transition rate $w_{i}(s)$ of each $i$th spin $s_{i}$ in the state $s$ is related to the probability $P(s_{i}=s)$ of finding the spin in that state by:

[TABLE]

That is, the ratio of the rates at which each spin transitions out of and into the state $s$, and correspondingly into and out of state $-s$, is equal to the ratio of probabilities of finding that spin in the state $s$ or $-s$ in the first place.

Each probability in (75) is given by the corresponding Boltzmann factor. Utilizing the mappings (21) and (24), we write the classical dual Hamiltonian as:

[TABLE]

where $H^{A}_{m}$ and $H^{B}_{n}$ are simply given by:

[TABLE]

From (75), we see that the Glauber dynamics of each $r_{m}$ are governed solely by the effective Hamiltonian $H^{A}_{m}$, and those of each $s^{n}_{j}$ are governed solely by the effective Hamiltonian $H^{B}_{n}$, as each Boltzmann factor contribution from uncorrelated spins divides out on the right-hand side. This is simply a restatement of the result that (6) under open boundary conditions is dual to $L^{3}$ single spins in a magnetic field and $(L-1)^{3}$ three-site Ising chains.

First, to find the autocorrelation function of each $A_{c}$, the expectation value of each $r_{m}$ is given by that of a single spin in a constant magnetic field. This problem is investigated in [75], with the following result:

[TABLE]

where $\alpha>0$ is a parameter corresponding to the transition probability between states. A priori, this parameter could be constant or may have a nontrival temperature dependence. Note that $\langle r_{m}(t)\rangle$ at long times corresponds with the result of (29). We may find (71) using the law of total expectation:

[TABLE]

where $\langle\cdot|\cdot\rangle$ here denotes conditional expectation, and the equilibrium probability of finding $r_{m}$ at $\pm 1$ is proportional to the corresponding Boltzmann factor. We therefore find our $A_{c}$ autocorrelation function:

[TABLE]

In the standard case of a constant, temperature independent, $\alpha$, the autocorrelation function (80) exhibits trivial high and low temperature limits. Specifically, for all $\beta$, the autocorrelation function (80) has an asymptotic long time decay set by $1/\alpha$. By contrast, the spin-spin correlations in conventional, infinite size, Ising chains and other systems become increasingly longer (in both spatial and temporal separation) as the temperature is lowered. In finite size systems (and especially for effective single site systems such as those associated with the decoupled terms $H^{A}_{m}$ on a lattice with open boundary conditions), by construction, long range spatial correlations are impossible and long time autocorrelations are more readily destroyed by thermal fluctuations. Effectively, in the equations of [75] determining the autocorrelations, the absence of coupling to spins that would appear in the infinite chain Ising model and are absent in the trivial finite size system can be emulated by setting $\beta=0$. The infinite chain case (typically examined by Glauber dynamics) leads to coupled differential equations whose solution exhibits a divergent correlation length (and time) as the temperature veers to zero. In our case, because the “chain” is finite (an effective single site in our calculation for the $A_{c}$ autocorrelations), the correlation length is trivially bounded by the effective length of the system (a single spin $r_{m}$) so that the infinite set of recursive coupled equations that appear in the standard textbook Markov chain analysis for Ising chains is truncated and no divergent correlations appear (in space and thus also weaker correlations in time). In finite size systems, autocorrelations are weaker than those on infinite size lattices. Given that even infinite length Ising chains do not exhibit divergent memory times at any non-zero temperature, it is no surprise, then, that in (80) we find that the correlations do not persist on arbitrarily long time scales. Even if the basic elementary “clock move” time $1/\alpha$ of Glauber dynamics is adjusted to be of an activated Arrhenius form of an exponential in the inverse temperature, (80) will not display long time correlations. Exactly at zero temperature, the righthand side of (80) is identical unity; at all other positive temperatures, no long memory times appear.

We find qualitatively similar behaviors for the $B^{\mu}_{v}$ operators. That is, similar to the $A_{c}$ autocorrelations, at finite temperatures, no long memory times appear in the autocorrelation function of $B^{\mu}_{v}$ operators. To explicitly compute the Glauber dynamics of the $B^{\mu}_{v}$ terms, we investigate the dynamics of the classical $s^{n}_{j}$ spins. The thermalization of each $s^{n}_{j}$ is not investigated explicitly in [75], but we may reproduce Glauber’s arguments easily by applying his techniques to the effective Hamiltonian $H^{B}_{n}$. Because each $B^{\mu}_{v}$ is mapped to some $s^{n}_{i}s^{n}_{j}$ for $i\neq j$, we seek expectation values of the form $\langle s_{i}(t)s_{j}(t)\rangle$, which are given by the functions $r_{i,j}$ in [75]. For the particular case of an ordinary Ising chain, each $r_{i,j}$ satisfies the differential equation in equation (31) of Glauber’s paper:

[TABLE]

where $\mathsf{T}_{2b}=\tanh(2\beta b)$. This differential equation holds a particularly simple form for the case of a periodic three-spin system: using $r_{i,i}(t)=\langle[s_{i}(t)]^{2}\rangle=1$ and $r_{i,j}=r_{j,i}$, we have for $r_{1,2}$:

[TABLE]

Similar differential equations are found for $r_{1,3}$ and $r_{2,3}$ by exchanging $r_{1,2}$ with $r_{1,3}$ and $r_{2,3}$ respectively. These three equations yield a system of three linear inhomogeneous differential equations, most easily solved in matrix form:

[TABLE]

We solve this system by diagonalizing the matrix. Its eigenvalues and eigenvectors are given by:

[TABLE]

In this eigenbasis, (83) is given by:

[TABLE]

where $X(t),Y(t),$ and $Z(t)$ are given by:

[TABLE]

Equation (85) consists of three uncoupled ordinary differential equations, which are easily solved:

[TABLE]

This yields the following solutions for the expectation values:

[TABLE]

To determine the constants $X_{0}$, $Y_{0}$, and $Z_{0}$, we demand that each expectation value is equal to its initial value at $t=0$:

[TABLE]

Solving this system, we find the solution for $r_{1,2}$ is given by:

[TABLE]

Similar expressions are given for $r_{2,3}$ and $r_{1,3}$ by respectively exchanging $r_{1,2}(0)$ and $r_{2,3}(0)$ or $r_{1,2}(0)$ and $r_{1,3}(0)$. Note that we may rewrite the final term of (90) as follows:

[TABLE]

In this way, we see that (90) also coincides at long times with the result of (29).

Since this expression doesn’t assume knowledge of the initial conditions $r_{i,j}(0)$, it can be used to calculate conditional expectation values as well. To finally obtain each autocorrelation function, we again employ the law of total expectation:

[TABLE]

Each probability is given by the respective Boltzmann constants, as determined by the effective Hamiltonian $H^{B}_{n}$. In this Hamiltonian, $H^{B}_{n}$ attains a value of $-3b$ for the two configurations corresponding to $s_{1}=s_{2}=s_{3}=\pm 1$, and attains a value of $b$ for all six remaining configurations. The former two states and two of the six latter states correspond to $s_{1}(0)s_{2}(0)=+1$, while the remaining four correspond to $s_{1}(0)s_{2}(0)=-1$. Therefore the two probabilities are:

[TABLE]

Since $r_{1,2}(t)=\langle s_{1}(t)s_{2}(t)\rangle$, $r_{1,2}(0)$ is determined immediately by the conditionals $s_{1}(0)s_{2}(0)=\pm 1$. Additionally, note that $r_{2,3}(0)$ and $r_{1,3}(0)$ will take on different values under these two conditionals: in particular, if $s_{1}(0)s_{2}(0)=+1$, then $H^{B}_{n}=-3b$ in the two configurations for which $s_{3}(0)=s_{2}(0)$, and $H^{B}_{n}=b$ for the two configurations in which $s_{3}(0)\neq s_{2}(0)$. On the other hand, $H^{B}_{n}=b$ for all four configurations of spins in which $s_{1}(0)s_{2}(0)=-1$, giving $s_{2}(0)s_{3}(0)$ an expectation value of zero. Letting $r_{i,j}^{\pm}$ denote the conditional expectation of $s_{i}s_{j}$ given $s_{1}s_{2}=\pm 1$, we arrive at the result:

[TABLE]

Utilizing (91), we therefore evaluate (92) as:

[TABLE]

Perusing (95), we indeed see that, similar to (80), for Glauber dynamics on a lattice with open boundary conditions, the autocorrelation function for the $B^{\mu}_{v}$ operators indeed does not exhibit divergent autocorrelation times. Even if $1/\alpha$ is made to be exponential in the inverse temperature (as in activated finite temperature dynamics), the correlations of (95) decay with a finite lifetime at all non-zero temperatures.

We next analyze a $p$-state generalization of the X-Cube model (the “$p$X-Cube model”).

9 Equilibrium Thermodynamics of the $p$X-Cube Model

We may easily generalize the preceding discussion of the ordinary X-Cube model to an X-Cube model built upon $\mathbb{Z}_{p}$ clock and shift operators, see, e.g., [35, 76]. That is, rather than considering qubits at each edge of an $L\times L\times L$ lattice, we consider $p$-qudits with associated $p$-dimensional Hilbert spaces $\mathcal{H}_{n}$, forming a total state space $\bigotimes_{n=1}^{N}\mathcal{H}_{n}$. In place of the ordinary Pauli operators $\sigma^{x}_{n}$ and $\sigma^{z}_{n}$ acting on each $n$th Hilbert space, following the quantum clock operators of [71] (see also the subsequent work of [77]), one may generalize the qubits of the X-Cube model to the $p$-clock operators $X_{n}$ and $Z_{n}$ (which are the $V$ and $U$ operators of Reference [71], respectively). These are traceless operators with eigenvalues $\omega^{m}$ for $\omega=e^{2\pi i/p}$ and $0\leq m\leq p-1$. $X_{n}$ and $Z_{n}$ are further defined by their relation:

[TABLE]

Unlike the ordinary Pauli operators, $X_{n}$ and $Z_{n}$ are not Hermitian for $p>2$. However, it is easy to verify that they are unitary: in a given $\mathcal{H}_{n}$, we may diagonalize either $X_{n}$ or $Z_{n}$ (but not both simultaneously) as $\operatorname{diag}(1,\omega,\ldots\omega^{p-1})$, in which case $X_{n}^{\dagger}$ or $Z_{n}^{\dagger}$ is given by $\operatorname{diag}(1,\overline{\omega},\ldots\overline{\omega}^{p-1})$, where the bar denotes complex conjugation. We therefore quickly see that $X_{n}^{\dagger}X_{n}=Z_{n}^{\dagger}Z_{n}=\mathds{1}$. Using this result and (96), we may also derive the relations:

[TABLE]

From above, and additionally Hermitian conjugating (96), we find:

[TABLE]

Additionally, it can easily be verified that $X_{n}Z_{n}$ is traceless by taking the trace of both sides of (96), and both $X_{n}^{m}$ and $Z_{n}^{m}$ can be verified to be traceless for $1\leq m\leq p-1$ by computing in their respective diagonal bases.

To generalize our $A_{c}$ and $B^{\mu}_{v}$ operators to arbitrary dimension $p$, we construct the operators $\mathcal{A}_{c}$ and $\mathcal{B}^{\mu}_{v}$ as depicted in figure 10. Note that each $\mathcal{A}_{c}$ and $\mathcal{B}^{\mu}_{v}$ operator commutes: for any given cube $c$ containing a vertex $v$, $\mathcal{A}_{c}$ and any of $\mathcal{B}^{\mu}_{v}$ share two common qudits. On one qudit, $\mathcal{A}_{c}$ and $\mathcal{B}^{\mu}_{v}$ utilize either $X_{n}$ and $Z_{n}$, or $X_{n}^{\dagger}$ and $Z_{n}^{\dagger}$. On the other qudit, $\mathcal{A}_{c}$ and $\mathcal{B}^{\mu}_{v}$ utilize either $X_{n}$ and $Z_{n}^{\dagger}$, or $X_{n}^{\dagger}$ and $Z_{n}$. When $\mathcal{A}_{c}$ and $\mathcal{B}^{\mu}_{v}$ are interchanged, the former will yield a factor of $\omega$ while the latter will yield a factor of $\overline{\omega}$, yielding an overall factor of 1.

Given these $\mathcal{A}_{c}$ and $\mathcal{B}^{\mu}_{v}$ operators, we define the corresponding $p$-clock X-Cube, or $p$X-Cube, Hamiltonian as:

[TABLE]

As in (6), the first sum is over each elementary cube $c$, while the second sum is over vertices $v$ and cardinal directions $\mu\in\{x,y,z\}$. Note that $\mathcal{A}_{c}^{\dagger}$ and $(\mathcal{B}^{\mu}_{v})^{\dagger}$ are included in (99) to ensure its Hermiticity. Since each operator in (99) commutes, we may write the corresponding partition function as:

[TABLE]

In the ordinary Pauli case of $p=2$, we evaluated this sum by noting that $e^{\beta\sigma}=\mathds{1}\cosh(\beta)+\sigma\sinh(\beta)$ for some operator $\sigma$ with $\sigma^{2}=\mathds{1}$. Here, the case is not so simple: let $\Sigma$ be an operator with $\Sigma^{p}=\mathds{1}$. Then, we have:

[TABLE]

where we have defined the functions $Q^{p}_{m}$ via power series and in closed form:

[TABLE]

$0\leq m\leq p-1$. These functions are natural generalizations of $\cosh$ and $\sinh$: for $p=2$, $Q^{2}_{0}$ is simply $\cosh$ and $Q^{2}_{1}$ is $\sinh$.

Using these results in (100), we have:

[TABLE]

We may simplify this expression further: since $\mathcal{A}_{c}^{p}=\mathcal{A}_{c}\mathcal{A}_{c}^{p-1}=\mathds{1}$, we have that $\mathcal{A}_{c}^{\dagger}=\mathcal{A}_{c}^{p-1}$, and we can rewrite the first two summation products (and the latter two similarly) as:

[TABLE]

Because this sum includes negative powers of $\mathcal{A}_{c}$, it contains redundancies: $\mathcal{A}_{c}^{-n}$ is the same as $\mathcal{A}_{c}^{p-n}$. We therefore wish to rewrite this sum in terms of only positive powers of $\mathcal{A}_{c}$. Toward this end, we first split the sum into three components:

[TABLE]

Then, we rewrite the first sum as a sum over positive powers of $\mathcal{A}_{c}$:

[TABLE]

Finally, we recombine powers of $\mathcal{A}_{c}$:

[TABLE]

In order to proceed, we once again introduce new functions to manage the algebra. Let $R^{p}_{m}$ be defined by:

[TABLE]

We may then rewrite (102) as:

[TABLE]

In defining the functions $R^{p}_{m}$, we have somewhat obscured the meaning of our calculation, but we have made the computation significantly easier by reducing the number of constraints needed to consider. As written in (108), each product term for a given $c$ or $(\mu,v)$ contains a sum over all powers zero through $p-1$ of the respective operator $\mathcal{A}_{c}$ or $\mathcal{B}^{\mu}_{v}$. The total product expansion will therefore contain exactly one linear term for every possible combination of operators $\mathcal{A}_{c}$ and $\mathcal{B}^{\mu}_{v}$ with each operator raised to a power zero through $p-1$. All terms not proportional to the identity will be traceless. Had we attempted to expand directly from (102), this product would also contain negative powers of $\mathcal{A}_{c}$ and $\mathcal{B}^{\mu}_{v}$, and the number of possible product combinations proportional to the identity would be significantly larger. From here, the calculation will depend on the choice of boundary conditions.

9.1 Open Boundary Conditions

Under open boundary conditions, no product of $\mathcal{A}_{c}$ operators can yield the identity. As can be seen in figure 10a, each $\mathcal{A}_{c}$ can only be canceled by adjacent $\mathcal{A}_{c}$ operators. But for product of $\mathcal{A}_{c}$ operators forming a connected section of elementary cubes, the $X_{n}$ operators lying at its boundary will appear in the product only once. We may therefore immediately reduce (108) to the form:

[TABLE]

As in the $p=2$ case, the remaining sum is performed over only the $(L-1)^{3}$ interior vertices of the system. Additionally, we see from figure 10 that the constraint (3) carries over to the general case. However, because each $\mathcal{B}^{\mu}_{v}$ may now be raised to nontrivial powers, we have the more general condition:

[TABLE]

Just as in the $p=2$ case, the only products of $\mathcal{B}^{\mu}_{v}$ operators proportional to the identity are those satisfying (110). Therefore, we obtain all terms proportional to the identity by choosing, for each of $(L-1)^{3}$ vertices, the power to which (110) is given. Each product of the form (110) carries a factor of $[R^{p}_{m}(\beta)]^{3}$. Therefore, we may expand the remaining product as:

[TABLE]

As before, the number $N$ of $p$-qudits in our system is given by $3L^{3}+6L^{2}+3L$. Because $\operatorname{Tr}[\mathds{1}]$ is given by $p^{N}$, (109) is finally given by:

[TABLE]

As before, we remark that the thermodynamic properties of the model captured by this partition function will accurately describe the bulk material in the thermodynamic limit of very large $L$, regardless of our choice of boundary conditions: any corrections due to constraints such as (58) and (59) will appear only at orders $L^{2}$ and higher.

If $p=2$, we expect this expression to yield the prior solution (17). Indeed, $Q^{2}_{0}(x)$ is simply $\cosh(x)$ and $Q^{2}_{1}(x)$ is $\sinh(x)$. $R^{2}_{0}$ and $R^{2}_{1}$ are therefore given by:

[TABLE]

By substituting into (112), we obtain:

[TABLE]

This matches exactly with our prior solution (17): for $p=2$, we have $\mathcal{A}_{c}^{\dagger}=\mathcal{A}_{c}$ and $(\mathcal{B}^{\mu}_{v})^{\dagger}=\mathcal{B}^{\mu}_{v}$, so (99) is identical to (6) with $a\rightarrow 2a$ and $b\rightarrow 2b$.

9.2 Cylindrical Boundary Conditions

Just as in the $p=2$ case, we gain two additional constraints by imposing cylindrical boundary conditions: letting our system become periodic in the $y$ and $z$ directions, we find:

[TABLE]

where $P^{x}_{i}$ refers to a plane of elementary cubes perpendicular to the $x$-direction in the first equality, and $\bar{P}^{x}_{i}$ refers to a plane of vertices perpendicular to the $x$-direction in the second equality. Once again, the only products of $\mathcal{A}_{c}$ and $\mathcal{B}^{\mu}_{v}$ operators proportional to the identity are those satisfying (110) and/or (115). Therefore, the products of $\mathcal{A}_{c}$ operators in (108) proportional to the identity are found by choosing, for each of $L$ $x$-planes, the multiplicity $m$ of the entire plane included in the product:

[TABLE]

Handling the $\mathcal{B}^{\mu}_{v}$ operators is slightly more complex: first, for any $x$-plane of vertices not satisfying (115), we must choose for each vertex to satisfy (110) for $0\leq m\leq p-1$. Second, for any vertex on an $x$-plane satisfying (115), we may choose to replace $(\mathcal{B}^{x}_{v})^{m}$ with $(\mathcal{B}^{y}_{v}\mathcal{B}^{z}_{v})^{p-m}$. We therefore construct each $\mathcal{B}^{\mu}_{v}$ product proportional to the identity as follows: first, for each of $L-1$ $x$-planes, we choose the degree $m$ with $0\leq m\leq p-1$ with which a given $x$-plane satisfies (115). Then, if $m=0$ for a given plane, we choose for each vertex on that plane to include a factor of $[R^{p}_{n}(\beta b)]^{3}$ with $0\leq n\leq p-1$ in order to satisfy (110). Next, if $m\geq 1$, we pick for each vertex on the plane to include either a factor of $R^{p}_{m}(\beta b)[R^{p}_{0}(\beta b)]^{2}$ corresponding to the use of $(\mathcal{B}^{x}_{v})^{m}$ for that vertex or a factor of $R^{p}_{0}(\beta b)[R^{p}_{p-m}(\beta b)]^{2}$ corresponding to the use of $(\mathcal{B}^{y}_{v}\mathcal{B}^{z}_{v})^{p-m}$ for the vertex. This gives the following product expansion:

[TABLE]

Finally, since the number $N$ of $p$-qudits in our system under cylindrical boundary conditions is $3L^{3}+2L^{2}$, the final partition function is given by:

[TABLE]

Once again, we check the $p=2$ case to ensure it matches (39). Using (113), we have:

[TABLE]

This is exactly the solution we previously derived, after rescaling $a\rightarrow 2a$ and $b\rightarrow 2b$.

9.3 Large-$p$ Limit

Another limit of particular interest is when $p$ becomes very large. In this case, the $\mathbb{Z}_{p}$ theory becomes approximated by a $U(1)$ theory (see, e.g., [71] for a discussion of clock models). Let $Q^{\infty}_{m}$ and $R^{\infty}_{m}$ denote the limits as $p$ goes to infinity of $Q^{p}_{m}$ and $R^{p}_{m}$, respectively. From its power series representation in (101), it can easily be seen that $Q^{\infty}_{m}$ is simply given by the $m$th term of the exponential power series:

[TABLE]

Then, each $R^{\infty}_{m}$ can quickly be verified to be given by modified Bessel functions by comparing Taylor series:

[TABLE]

That $R^{\infty}_{0}$ becomes the 0th modified Bessel function may not be surprising – since the eigenvalues of $\mathcal{A}_{c}$ and $\mathcal{A}_{c}^{\dagger}$ are given respectively by $\omega^{m}$ and $\overline{\omega}^{m}$, we may evaluate the partition function of a single elementary cube via:

[TABLE]

In the large $p$ limit, we may approximate this sum by an integral, in which the above sum is interpreted as a left-handed Riemann sum [71]:

[TABLE]

This result appears in particular in (112), in which the first nontrivial factor is simply the product of $L^{3}$ such Bessel functions in the large $p$ limit.

10 Implications of Dualities on the Nature of Fracton Excitations

It has been well appreciated, for some time by now, that the low energy excitations of the X-Cube model may propagate in a correlated manner in order to avoid further energy penalties. This intriguing feature raised the possibility of glassy dynamics and associated “protection” of quantum information that may be coded in the low energy states of this model [21, 23, 36, 37]. Given the results of our duality mappings, we now revisit these notions and point to a simple consequence of our dualities. As we demonstrated in section 8, our dualities imply that with (dual) thermal baths, no excessively long time autocorrelations may persist at positive temperatures. That is, the finite temperature autocorrelations display, at low energies (or, equivalently, low temperatures), the hallmarks of conventional activated dynamics. Our computed results do not feature any indications of exotic behaviors or particularly slow constrained dynamics. We caution anew that our results, invoking dualities, relate to (generally non-local) duals of Glauber thermal baths. Thus, the Glauber dynamics that we derived in section 8 might differ from those for other baths. However, if the baths are relatively featureless (as typical thermal baths are) then duality transformations might not be expected to alter the system dynamics. Additionally, it is sometimes possible to explicitly construct dualities for which local heat baths in the original model are mapped to local heat baths in the dual model. For a demonstration using the 2D toric code, see the supplementary material of Ref. [78].

The low energy dynamics are rooted in the character of the corresponding excitations. With that in mind, we wish to stress a simple conceptual point. Our duality mappings establish that the spectra of the X-Cube model and those of Ising chains are identical. Stated more precisely, the dualities (21) and (45) imply that for open and cylindrical boundary conditions, the X-Cube model has a spectrum which is none other than that of Ising chains, with degeneracies that differ only by a global power of two (as is also confirmed by our high temperature series results of (19) and (39)). All exact dualities are unitary maps that preserve the spectrum [69, 70]. Equivalently, the equivalence of the partition function,

[TABLE]

of two dual models implies that the spectra of the dual models is the same [2, 69, 70]. Here, $g(E_{n})$ denotes the degeneracy of each energy $E_{n}$. Thus, for both cylindrical and open boundary conditions, the spectrum of the X-Cube model including the degeneracy (modulo a global power of two) of all its low energy excitations is precisely the same as that of Ising chains. (From our results in section 7, the same holds true only in an asymptotic sense for the X-Cube model with periodic boundary conditions.)

Thus, if the arguments concerning immobility of low energy excitations do not involve the boundaries, one might expect that since the excitations map in a one-to-one manner between the dual models, the energetics of defects in the X-Cube model in the system bulk (including arguments favoring low energy dynamics of one type or another) will have exact counterparts for the classical Ising chains that are dual to the X-Cube model. However, because defects in standard classical Ising chains (i.e., domain walls) do not feature unusual dynamics, the same may be expected for their exact X-Cube duals. Dualities are, generally, non-local unitary transformations. Thus, a priori, one might anticipate that a sequence of states in which the energy is progressively altered or remains the same as defects locally move in a given system may involve, in its dual counterpart, a very different sequence of (non-local) moves for the corresponding defects in the dual system. With this in mind, we remark that Lieb-Robinson (LR) bound [79] type arguments (that fundamentally restrict the propagation of correlations) are suggestive of local defect motion in two systems that are dual to each other, so long as these systems exhibit local interactions and local operators may be used to define/measure the defects [80]. Thus, any local dynamics of defects in a given system (with this locality also required by the LR bounds) that change the energy in some way mandate the corresponding existence of local dynamics of any locally discernible defects in the dual model. By the unitary character of the duality transformation, the dynamics of the defects in the dual model alter the energy by exactly the same amount. In the Appendix, we very qualitatively discuss in some more detail several aspects of low energy excitations of the X-Cube model. Unlike the calculations in our work thus far, the arguments in that Appendix and in the current Section are by no means rigorous, and are only suggestive.

11 Conclusion and Outlook

In this paper, we investigated a prototypical fractonic model, the X-Cube model and its $p$-state generalizations, at finite temperature. We computed the partition function of the models in closed form for open and cylindrical boundary conditions, and we showed that these partition functions agree in the thermodynamic limit and agree with that for periodic boundary conditions. These calculations provide compelling evidence that the X-Cube models have a straightforward thermodynamic limit, insensitive to boundary conditions as one usually assumes in elementary statistical mechanics. Moreover, we find the absence of finite temperature phase transitions and thermal fragility in these models [39, 40, 41, 49, 50, 51, 52, 53]. From a dynamical point of view, the elementary excitations of the X-Cube model display highly constrained mobility. Thus, one may hope that the approach to equilibrium might be extraordinarily slow (“glassy behavior”). To explore this possibility, we leveraged a duality transformation to set up a simple Glauber model of equilibration for the X-Cube model, and found that conventional activated dynamics may appear, as opposed to glassy dynamics. It is possible that the Glauber dynamics that we study are not generic, but they suffice to demonstrate as a proof of principle that glassy dynamics are not mandatory for fractonic matter.

Why is fractonic matter susceptible to thermal fragility? At the most basic level, the problem is that all the efforts in designing fracton models go into engineering the energy barrier for low-lying excitations (as in, for example, this type of excitation can only be created in quartets), but the issue at finite temperature is, of course, the free energy. To design matter at finite temperature, one must keep track of both the energy and the entropy, or, microscopically, the energy levels and the density of states. Now, more concretely, thermal fragility in numerous models stems from the same specific source - that of dimensional reduction. Entropic effects prohibit stable finite temperature order in conventional low dimensional systems. In a similar manner, thermal fluctuations eradicate stable orders in models dual to these low dimensional systems. This underscores the importance of entropic effects in this family of models of topological quantum matter. We are currently further investigating the notion of effective dimensionality and encoding from the point of view of bond algebras; we reserve further remarks for a future publication.

Looking forward, we believe that our results highlight the need for an extended set of designing principles towards topological quantum memories.

12 Acknowledgments

This research was largely supported by NSF Grant 1411229 (CMMT). We also gratefully acknowledge NSF PHY-1607611 for work at the Aspen Center for Physics and grant NSF PHY-1748958 for work at the Kavli Institute for Theoretical Physics (KITP).

Appendix A

Ground States and Low Energy Excitations:

A Review and General Remarks

In this Appendix, we review and further discuss the ground states of the X-Cube model and their low energy excitations. The aim of this Appendix is to ground the general considerations of section 10. Unlike most other sections of this paper, the following discussion is largely qualitative.

The X-Cube model Hamiltonian of (6) is a sum of commuting terms (so-called “stabilizers”). Any ground state $|\psi_{0}\rangle$ of the X-Cube model satisfies the “frustration free” condition,

[TABLE]

From (6), $|\psi_{0}\rangle$ clearly has the lowest possible energy of the system. A ground state is given by

[TABLE]

where ${\cal{N}}_{0}$ is a normalization factor and $|0\rangle$ is a simultaneous ($+1$) $\sigma^{z}$ eigenstate of each link. The state $|0\rangle$ is a trivial eigenstate of all $B^{\mu}_{v}$ operators with eigenvalue $+1$ and the commuting projectors $\frac{1}{2}(\mathds{1}+A_{c})$ ensure that $|\psi_{0}\rangle$ is an eigenstate of $A_{c}$ (with eigenvalue $+1$) for all cubes $c$. Thus, Eq. (125) is satisfied. The X-Cube model exhibits an exponential (in system length $L$) degeneracy. Notice that such an exponentially large degeneracy may appear in classical models that do not display topological order [56, 81].

We may construct particular excited states by applying $\sigma_{n}^{z}$ to $|\psi_{0}\rangle$. The operator $\sigma_{n}^{z}$ commutes with all but four $A_{c}$ operators. Thus, $\sigma_{n}^{z}|\psi_{0}\rangle$ remains an eigenstate of these operators (with eigenvalue $+1$). However, $\sigma^{z}_{n}$ will anticommute with the four $A_{c}$ operators containing $\sigma^{x}_{n}$. Thus,

[TABLE]

It follows that $\sigma^{z}_{n}$ flips the eigenvalues of $A_{c}$ for each cube connected to the $n$th qubit. This can be expanded by instead considering a “membrane” of $\sigma^{z}_{n}$ operators as in figure 11, flipping the eigenvalues of the four cubes at the corners of the membrane. The energy cost of creating these four excitations is $8a$, from flipping the eigenvalue of four $aA_{c}$ terms in (6) from $-a$ to $+a$.

Three of these excitations can be moved “off to infinity” by increasing the size of the membrane operator at no additional energy cost. A single localized excitation cannot move freely by application of a local operator. However, such an excitation can move while creating two additional excitations as in figure 12, at an energy cost of $4a$. Note that a pair of excitations can move freely in two dimensions, as shown in figure 13, but not all three. In addition, a quartet of excitations is fully mobile in all three dimensions via applications of appropriate $\sigma^{z}$ operators.

As argued in section 10, we can construct spectra for the open and cylindrical X-Cube models identical to those of their classical duals using nonlocal membrane operators.

First, consider the case of open boundary conditions. The mapping (21) identifies each cubic operator $A_{c}$ with a bond variable $r_{m}$ of an $L^{3}+1$ site open Ising chain. The duality suggests that the spectra of $A_{c}$ operators in the open X-Cube model ought to be identical to that of the open Ising chain, in which each $n$th energy level can be achieved by choosing any arbitrary arrangement of $n$ “bad bonds” corresponding to $r_{m}=-1$ (see figure 14). In particular, the duality implies that the excitations of the $A_{c}$ operators in the open X-Cube model are given not by four-fold excitations at the corners of a membrane operator, but by any arbitrary arrangement of excitations corresponding to $A_{c}|\psi\rangle=-|\psi\rangle$.

Indeed, we can explicitly construct each of these excited states using nonlocal membrane operators: to place a lone excitation in any particular cubic location, start with the ground state $|\psi_{0}\rangle$, and apply the membrane operator $\cal{M}$ described in figure 11 with one corner at the desired locaton. Then, expand $\cal{M}$ to move the three extraneous excitations to the boundaries, as in figure 15. With $\cal{M}$ extending to the boundaries at all but one corner, only one excitation remains. By overlaying multiple such operators, we may place any number of excitations in any possible arrangement: for instance, we can create four excitations in a non-rectangular arrangement (figure 16), or we can move a single excitation at no energy cost (figure 17).

The case of cylindrical boundary conditions is not as simple. The mapping (45) identifies each plane of cubic operators with an $L^{2}$ site periodic Ising chain, with each given $A_{c}$ mapped to a single bond $r^{i}_{m}r^{i}_{m+1}$. Excitations in the bonds of a periodic Ising chain come in multiples of two: any domain wall must be accompanied by another domain wall to return to the original spin direction. The duality implies that the excitations of the cylindrical X-Cube model are also given by arbitrary arrangements of excitations, so long as each $x$-plane has excitations in multiples of two. We may construct such excited states using two nonlocal membrane operators: first, to construct two cubic excitations in the same $x$-plane and lying along a line, simply apply to the ground state a membrane operator $\cal{M}$ perpendicular to the $y$ or $z$ direction, with two corners of the membrane lying in the desired plane and the other two corners extending to the open $x$ boundary. Then, to create two excitations in any location within the same $x$-plane, simply multiply two such perpendicular membrane operators sharing a common corner, as in figure 18. This procedure may be continued to construct arbitrary multiple-of-two arrangements of excitations within each $x$-plane of the lattice.

We have previously shown that the free energy density of the X-Cube model is independent of any choice in boundary conditions in the thermodynamic limit. However, this does not explicitly mean that the spectrum is independent of our choice in boundary conditions: indeed, just as the open and periodic one-dimensional Ising chain have different spectra, the open and cylindrical X-Cube models also do not have the same spectra. For this reason, we do not suggest that the periodic X-Cube spectra is the same as that of the open or cylindrical models — in fact, under fully periodic boundary conditions, there exists no operator (local or nonlocal) which can create a single cubic excitation in a ground-state wavefunction [22]. That said, the spectra of the open and cylindrical systems do suggest that we must be careful when discussing the nature of excitations in the X-Cube and other similar models. In particular, while the above depiction of fractons in the model discusses creation and mobility of fractons using local operators, the results of this paper suggest that it may not be sufficient to consider local operators alone. Once nonlocal operators are considered, it’s perfectly clear that the open boundary X-Cube model has no special thermodynamics or constrained mobility: all cubic excitations are completely decoupled from one another, and any combination of cubic excitations may be realized from any other via some combination of nonlocal membrane operators. The same is true of the cylindrical X-Cube model: while the additional constraints (33) introduce some correlation among $A_{c}$ expectation values within a given $x$-plane, these correlations vanish in the thermodynamic limit, in the same way as the bond variables of the periodic Ising chains to which the cylindrical X-Cube model is dual.

While we cannot provide a closed form partition function, spectrum, or free energy of the periodic X-Cube model at finite $L$, the results for open and cylindrical boundary conditions at least indicate that we must be careful with how we discuss excited states and their dynamics under periodic boundary conditions as well. The lack of finite temperature phase transitions in all models, as well as the seeming irrelevance of order $L^{2}$ constraints such as (33), (34), (58), and (59) in the thermodynamic limit, suggest that the dynamics and equilibrium thermodynamics of the periodic X-Cube model may be simpler than originally thought. An earlier analysis [21] found that “when a zero temperature type I fracton model is placed in contact with a finite temperature heat bath, the approach to equilibrium is a logarithmic function of time over an exponentially wide window of time scales.” By contrast, the results of section 8 for the open X-Cube model with Glauber dynamics suggest, instead, an exponential decay at all times.

What’s more, if the physical motivation for studying the X-Cube model is as a stabilizer code with a robust quantum memory, then it’s likely that a physical realization of the X-Cube model has boundary conditions closer to the open boundary conditions discussed in section 5 than the fully periodic boundary conditions discussed in section 7. In this case, our results suggest that such a physical system may be subject to the same thermal fragility as the Kitaev model [49].