Glassy quantum dynamics in translation invariant fracton models

TL;DR

This paper reveals that translation invariant fracton models exhibit glassy quantum dynamics, including suppressed charge mobility and logarithmic relaxation, resembling disordered systems but without quenched disorder.

Contribution

It demonstrates that fracton models inherently display glassy dynamics and subdiffusive behavior, expanding understanding of quantum relaxation in disorder-free systems.

Findings

Charge mobility is exponentially suppressed at low temperatures.

Relaxation to equilibrium follows a logarithmic time dependence.

Type II fracton models show super-exponential divergence of relaxation times.

Abstract

We investigate relaxation in the recently discovered "fracton" models and discover that these models naturally host glassy quantum dynamics in the absence of quenched disorder. We begin with a discussion of "type I" fracton models, in the taxonomy of Vijay, Haah, and Fu. We demonstrate that in these systems, the mobility of charges is suppressed exponentially in the inverse temperature. We further demonstrate 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. Generalizing to the more complex "type II" fracton models, we find that the charges exhibit subdiffusion upto a relaxation time that diverges at low temperatures as a super-exponential function of inverse temperature. This behaviour is reminiscent of "nearly localized" disordered systems, but occurs with a translation invariant three-dimensional Hamiltonian. We also conjecture that fracton models with conserved charge may support a phase which is a thermal metal but a charge insulator.

Equations

$$H_{XC}=-\sum _c{A}_c-\sum _{v,k}{B}_c^{(k)},$$

$$H=-J\sum _c{A}_c-\sum _{v,k}{B}_c^{(k)}+\Lambda \sum _i{\sigma }_z+\lambda \sum _i{\sigma }_x,$$

$$H=-\sum _c{A}_c-\sum _{v,k}{B}_c^{(k)}+\Lambda \sum _i{\sigma }_z,$$

$$n_f\sim e^{-\frac{W}{2T}},n_b\sim e^{-\frac{W}{T}},n_c\sim e^{-\frac{2W}{T}},$$

$$\Gamma \sim n_f{e}^{-n}\sim n_f{e}^{-W/\tilde{\Lambda }}.$$

$$\Gamma \sim n_f{n}_b\sim n_f{e}^{-W/T},$$

$$\Gamma \sim n_f\max \left(e^{-W/T},e^{-W/\tilde{\Lambda }}\right)\sim n_f{e}^{-W/T}.$$

$$W{e}^{-W/T}>T,$$

$$\Gamma \sim n_f{e}^{-W/\tilde{\Lambda }},$$

$$E_c\sim 2W{n}_c$$

$$E_b\sim W{n}_b$$

$$E_f\sim \frac{W}{2}{n}_f$$

$$\frac{d{E}_b}{dt}=-{\Lambda }^2\left(3n_b^2+n_b{n}_f-n_b{n}_f^2-n_f^3-2n_f^4\right)=-\frac{d{E}_f}{dt}.$$

$$\frac{d{T}_b}{dt}$$

$$\frac{d{T}_f}{dt}$$

$$T_f(t)=-\frac{W/2}{\log \left(\frac{6{\Lambda }^2}{W}t+b\right)-\frac{2W}{T_b^{(0)}}},$$

$$T_b(t)=\frac{W}{\log \left(\frac{3{\Lambda }^2}{W}t+e^{W/T_b^{(0)}}\right)},$$

$$0\le t\lesssim t^{\ast }=\frac{W}{6{\Lambda }^2}\exp \left(\frac{W}{T_b^{(0)}}\right).$$

$$T_b(t)\sim \frac{W^2}{{\Lambda }^{2}t}{e}^{W/T_b^{(0)}},t>t^{\ast }.$$

$$0\le t\lesssim \frac{cW}{{\Lambda }^2}\exp \left(\frac{W}{T_b^{(0)}}\right),$$

$$\delta T(t)\sim \exp \left(-\frac{{\Lambda }^2}{W}{e}^{-\frac{W}{T_b^{(0)}}}t\right),$$

$$\frac{d{T}_f}{dt}=\frac{4{\Lambda }^2{T}_f^2}{W^2}\left(3\frac{n_b^2}{n_f}+n_b-n_b{n}_f-n_f^2-2n_f^3\right),$$

$$T_f(t)=-\frac{W/2}{\log \left(\frac{6{\Lambda }^2}{W}t+b\right)-\frac{2W}{T}},$$

$$0\le t\lesssim t^{\ast }=\frac{W}{6{\Lambda }^2}\exp \left(\frac{W}{T}\right),$$

$$T_f(t)\approx -\frac{W/2}{\frac{2{\Lambda }^2}{W}{e}^{-W/T}t+\log \left(3e^{-W/T}\right)},t>t^{\ast }.$$

$$\delta T\sim \exp \left(-\frac{4{\Lambda }^{2}t}{W}{e}^{-W/T}\right),$$

$$H=-J\sum _c{A}_c-\sum _{v,k}{B}_c^{(k)}+\Lambda \sum _i{\sigma }_z+\lambda \sum _i{\sigma }_x$$

$$H=-J\sum _c{G}_c^X-\sum _c{G}_c^Z$$

$$t={\left(R\right)}^{\frac{cW}{\tilde{\Lambda }}};\tilde{\Lambda }=\min (T,\Lambda ).$$

$$t_{equilibrate}\sim \exp \left(+c^{\prime }\frac{W^2}{T^2}\right),$$

$$\frac{dy}{dt}=a{y}^2{e}^{\pm \frac{1}{y}}.$$

$${\int }_{y(0)}^{y(t)}\frac{dz}{z^2}{e}^{\mp \frac{1}{z}}=at.$$

Full text

I Introduction and Motivation

Recent years have witnessed an explosion of interest in many body localization (MBL) Anderson (1958); Gornyi et al. (2005); Basko et al. (2006); Žnidarič et al. (2008); Pal and Huse (2010); Imbrie (2016); Nandkishore and Huse (2015), whereby isolated quantum systems with quenched disorder can exhibit ergodicity breaking and fail to thermalize even at infinite times. Interest in the phenomenon is in part intrinsic (e.g. many body localized systems can support entirely new types of quantum order Huse et al. (2013); Vosk and Altman (2014); Pekker et al. (2014); von Keyserlingk et al. (2016)), in part practical (e.g. MBL systems can serve as ideal quantum memories Nandkishore and Huse (2015)), and in part due to the tantalizing connections to other fields. For example, MBL can be viewed as a type of “ideal quantum glass” which does not reach equilibrium at infinite times—can insights from MBL inform the study of classical structural glass? (This particular question has been the focus of intensive work Kagan and Maksimov (1984); Grover and Fisher (2014); Schiulaz et al. (2015); Yao et al. (2016); Papic et al. (2015); De Roeck and Huveneers (2014a, b); Hickey et al. (2016) without any definitive conclusions). Recent experimental progress Kondov et al. (2011, 2015); Schreiber et al. (2015); Choi et al. (2016); Bordia et al. (2016); Álvarez et al. (2015); Choi et al. (2016) has further intensified interest in the field. However, theoretical insights have been hard won, due to the fundamental challenges of describing a non-ergodic, non-equilibrium, strongly correlated, highly disordered phase. Indeed most theoretical progress has required either a new idea (such as the notion of emergent integrability Serbyn et al. (2013); Huse et al. (2014); Ros et al. (2015); Imbrie (2016); Geraedts et al. (2016)), a new technique (such as dynamical versions of real space renormalization Vosk and Altman (2013))—or a new class of models which provide access to a new and formerly unforeseen phenomenology (e.g. Vasseur et al. (2015)).

The study of topological phases in three spatial dimensions has brought to light a new class of models Chamon (2005); Bravyi et al. (2011); Haah (2011); Yoshida (2013); Vijay et al. (2015, 2016) which have remarkable properties. These exactly solvable three-dimensional (3D) lattice models have ground states that exhibit a sub-extensive topological degeneracy on the 3D torus and possess point-like excitations, dubbed “fractons” Vijay et al. (2015) that cannot move without creating additional excitations. Such systems are related by duality to spin models with symmetries along lower-dimensional subsystems Vijay et al. (2016); Williamson (2016) and were recently classified into Type I and Type II fracton phases. In the Type I phases, fractons are created by the application of a membrane operator and pairs of fractons form composite topological excitations that can move along lower-dimensional subsystems. In the Type II phases, fractons are created by the application of a fractal operator and all topological excitations are strictly immobile. The fracton models that have been introduced to date have all involved discrete symmetries, although there does not in principle appear to be any obstruction to constructing continuous fracton models—indeed a stimulating series of works from Pretko Pretko (2016a, b) appears to reproduce much of the fracton phenomenology within a continuum field theory with $U(1)$ symmetry. Layer constructions of these phases have also recently been advanced Vijay (2017); Ma et al. (2017).

In this work, we combine emerging ideas from research into (fracton) topological phases and MBL by studying the dynamical behaviour of fracton models, revealing an intimate and provocative connection between the two fields. We begin by discussing type I fracton models at finite energy density and demonstrate that in these models, the mobility of charges is suppressed exponentially in the inverse temperature. When a type I fracton model prepared in its ground space (i.e. at zero energy density) is placed in contact with a finite temperature heat bath, we show that the equilibration exhibits $\log t$ behaviour over an exponentially wide window of time scales—a classic signature of glassy dynamics (see e.g. Amir et al. (2011) and references contained therein). We emphasize that this glassy quantum dynamics occurs in a three-dimensional, translation invariant Hamiltonian. We then turn to type II fracton models, and demonstrate that charges exhibit subdiffusion upto a relaxation time that diverges at low temperatures as a super-exponential function of temperature, reminiscent of “near MBL” disordered systems Gopalakrishnan and Nandkishore (2014); Gornyi et al. (2016a). Finally, we conjecture that fracton models with conserved $U(1)$ charge could realize exotic three-dimensional phases that are thermal metals but charge insulators.

Our work has striking implications for MBL, for the study of fracton phases, and for possible technological applications of both. For MBL, our work illuminates new connections to glasses, introducing a new class of models that exhibit glassy quantum dynamics with translation invariant Hamiltonians. It may also inform investigations of localization and glassy dynamics in higher dimensions Nandkishore (2014); De Roeck and Huveneers (2016); Gornyi et al. (2016a). For the field of fracton phases, our work reveals that not only do these models have an unusual ground state structure, they also support rich quantum dynamics, thus opening a new line of investigation for three-dimensional topological phases. Practically speaking, our work also uncovers a new route to information storage, as well as identifying a potential class of three-dimensional phases that are thermal metals but charge insulators. These last could have applications e.g. in high density electronics, where the problem of heat dissipation currently constrains possibilities for miniaturization.

II Fracton Models

Fracton topological phases are a new class of three-dimensional phases of matter that display features that go beyond those familiar from gauge theory. These phases can be obtained as the quantum duals of three-dimensional systems with symmetries along lower dimensional sub-systems, specifically along planes and fractals. A unified framework, based on a generalized lattice gauge theory, for fracton topological order was recently proposed in Vijay et al. (2016). Given the novelty of these phases, in this section we provide a self-contained exposition of fracton systems, focusing primarily on specific examples to elucidate the features most relevant to the dynamics.

Fracton phases arise in exactly solvable lattice models in three spatial dimensions and exhibit a sub-extensive topological ground state degeneracy on the 3D torus. The distinguishing feature of these systems is the presence of point-like fractional excitations—fractons—that are fundamentally immobile i.e., they cannot move without creating additional topological excitations. In contrast with anyons in two-dimensional topologically ordered systems, where anyons are created at the ends of a Wilson line and are thus allowed to move by application of a local line-like operator, there exists no local line-like operator that creates a pair of fractons. Instead, fractons are created at the ends of membrane or fractal operators, leaving a single fracton immobile. A classification scheme for fracton topological order was recently proposed in Vijay et al. (2016), where these systems were divided into type I and type II phases.

Type I fracton phases, such as the X-Cube model discussed below, host fracton excitations at the ends of membrane operators. While single fractons are immobile, bound-states of fractons form composite topological excitations that are free to move along lower-dimensional subsystems such as a line or a plane. There may also exist additional quasi-particles that are confined to move only along lower-dimensional subsystems. In type II phases, such as Haah’s code Haah (2011), fractons are created by the application of fractal operators and all topological excitations are strictly immobile. The latter feature leads to a fundamental difference between the dynamics of type I and II fracton phases, which we now consider separately in the following sections.

III Type I Fracton Models

The physics of type I fracton topological order is best illustrated through the example of the X-Cube model Vijay et al. (2016) that displays the essential features of these phases. The X-Cube model is an exactly solvable lattice model defined on a cubic lattice with Ising spins living on each link. The Hamiltonian is

[TABLE]

where the first term is the sum over all cubes of a twelve-spin $\sigma_{x}$ interaction and the second term is the sum over all vertices of planar four-spin $\sigma_{z}$ interactions as depicted in Fig. 1. In contrast with two-dimensional topologically ordered states, such as the Toric Code, that have a finite and constant topological ground state degeneracy on the two-torus, the ground state of the pure X-Cube model on the three-torus of linear dimension $L$ has a sub-extensive topological ground state degeneracy $D$, where $\log_{2}D=6L-3$.

A fracton is created by flipping the eigenvalue of the cubic interaction term. However, there exists no local operator that can create a single pair of fractons. Indeed, applying a $\sigma_{z}$ operator to a link flips the eigenvalues of the four cubes sharing that link. Acting on the ground state by $\sigma_{z}$ along a membrane operator $\hat{\mathcal{M}}$ creates four fractons at the corners of the membrane, as shown in Fig. 2i. A single fracton, denoted $e^{(0)}$ (where the superscript denotes that it is a dimension-0 excitation), is thus fundamentally immobile, as moving it would create additional fractons. This is the fundamental “superselection” rule Kim and Haah (2016) that will lead to glassy dynamics. Pairs of fractons are however free to move by repeated application of local membrane operators. A straight Wilson line of $\sigma_{z}$ operators creates a pair of fractons at each end—each pair is a composite excitation that can move in two dimensions, and which we refer to as a dimension-2 (dim-2) excitation $e^{(2)}$, as shown in Fig. 2ii. In the X-Cube model, there exist additional dimension-1 excitations ($m^{(1)}$) created at the ends of a Wilson line of $\sigma_{x}$ operators, that are mobile along one-dimensional sub-manifolds. Henceforth we will refer to the fully mobile four-fracton composites, created by single $\sigma_{z}$ operators, as the topologically neutral sector, and the lower dimensional excitations (fractons and $e^{(2)}$’s) as the topologically charged sector.

In the following sections, we will focus specifically on the X-Cube model in the presence of transverse fields,

[TABLE]

where $i$ goes over all links in the cubic lattice. Since the pure X-Cube model Eq. (1) is a sum of commuting projectors, the relative coefficient $J$ simply sets the energy scale between the $e$ and $m$ excitations when $\Lambda,\lambda=0$. In the presence of the transverse fields, we expect that the fracton phase will survive up to some finite $\Lambda/J$ and $\lambda/J$ since this is a gapped phase of matter that is stable to local perturbations Bravyi et al. (2010); Bravyi and Hastings (2011). In the limit of large transverse fields however, the fracton topological order will be destroyed, but the precise nature of the transition between the fracton phase and the trivial paramagnetic phase has yet to be understood Vijay et al. (2016). Since we are interested in the dynamics within the fracton phase, allowing only for weak local perturbations, we set $J=1$, noting that our analysis holds as long $J$ is $O(1)$. In addition, we will first set $\lambda=0$ and consider only the dynamics of the fractons and their composites. After analysing this sector, we will comment on the consequences of a non-zero $\lambda$, which would allow the $m^{(1)}$ particles to hop as well. The Hamiltonian pertinent for the following discussions is thus

[TABLE]

with the perturbation strength $\Lambda\ll 1$. We note that while we are focusing on the specific example of the X-Cube model, the results presented here hold broadly for all type I fracton phases.111 The X-Cube model lacks an exact electromagnetic duality, present in other type I models e.g., the Majorana cube or Checkerboard models, which have two distinct fracton excitations, $e^{(0)}$ and $m^{(0)}$. However, we can always choose to initialize the system with only one topological excitation or perturb the system only by terms that allow one of the excitations to hop.

III.1 Type I Fractons at Finite Energy Density

We begin our discussion of dynamics in fracton models by considering the X-Cube model Eq. (3) at finite energy density. Since we have switched off the term that would allow $m$’s to hop, we have three kinds of excitations with dynamics in our system: neutral composites, which, being fully mobile and created by local terms, are three-dimensional bosons; dim-2 excitations $e^{(2)}$, which are two-dimensional bosons Vijay et al. (2016); and the topologically charged fractons. A word on notations—since both the neutral composites and $e^{(2)}$’s are bosonic, we will henceforth refer to the former as the composite ($c$) sector and the latter as the bosonic ($b$) sector, with fractons ($f$) sometimes also referred to as the (topologically) charged sector.

From the preceding general discussion of the model, it is clear that each fracton hop is accompanied by the creation of two additional fractons, and so energetically costs an amount $W=4$. We will henceforth refer to this charge gap simply as $W$, since the analysis applies also to generalizations of Eq. (3) that are in the same phase, but perhaps with a different charge gap. A single fracton hop is depicted in Fig. 3: starting with an isolated fracton, we can move this over by one site by acting by a single $\sigma_{z}$ operator on the link adjacent to the fracton. This creates two additional topological excitations; however, this pair is a dim-2 $e^{(2)}$ excitation and can be moved off to infinity at no additional energy cost. Thus, each hop takes the system off energy shell. To fully understand the relaxation in the fracton sector, we must thus take into account the fully mobile neutral composites and the dim-2 excitations. These sectors act as a thermalizing bath for the fractons as rearrangements within these sectors allow the system to come back on energy shell after each hop.

Before studying the non-equilibrium dynamics, we consider the system—comprised of fractons, $e^{(2)}$’s, and neutral composites—in equilibrium, at some temperature $T\ll W$. The local energy scale in the topologically charged sector is $W$, where $W$ is the charge gap, i.e., the cost of creating two fractons. Further, the fracton sector is coupled to a dense, high temperature bath of neutral composites which hop at a rate $\Lambda\ll W$. Here each species is gapped, with a gap of $W/2$, $W$, and $2W$ for creating a single fracton, $e^{(2)}$, and composite respectively. Additionally, since the composites and $e^{(2)}$’s are neutral, with the fractons carrying only a $\mathbb{Z}_{2}$ charge, the density and temperature of these excitations cannot be controlled independently; rather, the equilibrium temperature uniquely determines the density of each species,

[TABLE]

where $n_{f}$ is the density of fractons, $n_{b}$ is the density of dim-2 bosonic excitations, and $n_{c}$ is the density of the composites.

Due to the form of the perturbation, $\sigma_{z}$, the only processes that lead to an exchange of energy between the three sectors are those where the total $\mathbb{Z}_{2}$ charge along each row and each column of the cubic lattice is conserved. For instance, in Fig. 3, the initial and final topological charges along each column and each row are preserved (modulo 2). From all possible on-shell processes by which the three sectors can exchange energy, we will consider only two body processes that are up to second order in $\Lambda$ (shown in Figs. 12 and 5) since all others will be further suppressed either by density factors or by the perturbation strength. As an illustration, Fig. 12 depicts the processes where a composite breaks into two $e^{(2)}$’s and where two $e^{(2)}$’s combine into a composite.

Within the thermalizing bath, the composites and $e^{(2)}$’s can either hop at a rate $\Lambda$ or scatter off each other while remaining locally on-shell, at a rate $\Lambda^{2}/W$. Since we are interested in the regime where $\Lambda\ll W$, the heat bath has a narrow local bandwidth $\sim\Lambda$, determined primarily by the hopping. As discussed earlier, a single fracton hop takes the system off energy shell by an amount $W$, and hence a single rearrangement within the bath cannot place the system on shell. The traditional analysis of localized systems coupled to narrow bandwidth baths Gopalakrishnan and Nandkishore (2014); Parameswaran and Gopalakrishnan (2017); Nandkishore and Gopalakrishnan (2016) makes use of many body rearrangements in the charge sector to obtain a relaxation rate that is power law slow in the bandwidth of the bath. However, in that setting the charge sector admits local re-arrangements that are “uphill” or “downhill,” which may be combined into a many body re-arragement that is off shell by much less than $W$. The present situation differs in that every local re-arrangement in the charge sector is “uphill” in energy, so the dominant relaxation mechanism from Ref. Gopalakrishnan and Nandkishore (2014); Parameswaran and Gopalakrishnan (2017); Nandkishore and Gopalakrishnan (2016) does not apply. Since the maximum energy the bath can provide is $\tilde{\Lambda}=\min(\Lambda,T)$, the bath must be probed $n\sim W/\tilde{\Lambda}$ times for the fractons to borrow enough energy to perform a single hop. This leads to a relaxation rate in the charge sector that is exponentially slow in $W/\tilde{\Lambda}$, scaling as

[TABLE]

Alternatively, we could make use of the extensive nature of the many body bandwidth in the bath and could use an $n\sim W/\tilde{\Lambda}$ particle rearrangement in the bath to place fracton rearrangements on shell. However, this yields Abanin et al. (2015) an exponentially slow relaxation of the same form as Eq. (5), up to subleading prefactors.

The fractons, however, have an additional channel through which they can hop in an on-shell manner. This process, depicted in Fig. 4, requires the presence of an $e^{(2)}$ excitation in the vicinity of the fracton. During this process, the fracton hops once, destroying the neighbouring $e^{(2)}$, and hops again, thereby returning the system on-shell by creating an $e^{(2)}$ particle. Since this process is mediated by the $e^{(2)}$’s, the rate at which it proceeds is additionally suppressed by the density of the bosons $n_{b}$

[TABLE]

Thus, the fractons have two possible hopping channels—either by borrowing energy from the bath or by using the bosons as intermediaries. The faster rate will dominate and so

[TABLE]

This behaviour is in marked contrast with the usual case of activated transport. In typical gapped phases, the relaxation displays an Arrhenius law $\Gamma_{A}\sim n\sim e^{-\Delta/T}$, where $n$ is the density of charge carries and $\Delta$ is the charge gap. This exponential slowness of the relaxation is governed primarily by the exponential rarity of charge carriers, since the mobility of these excitations is generically $\sim O(1)$. In type I fracton models, while the relaxation is suppressed by the density of fractons $n_{f}$, which is exponentially small in the inverse temperature, it displays an additional suppression due to mobility, which is also exponentially small in the inverse temperature. This is conceptually different from Arrhenius behaviour. Additionally, we observe that insofar as $T$ is the potentially ergodicity restoring parameter here and the relaxation rate is exponentially small therein, Eq. (7) corresponds to “asymptotic many body localization” in the taxonomy of Ref. De Roeck and Huveneers (2014a, b).

We note that hole burning, a scenario where, after a certain number of fracton hops, there is a depletion in the density of composites causing the remaining fractons to essentially be frozen until the composite sector thermalizes, does not arise here. This is due to the exponentially long time scales, $t\sim e^{W/T}$, over which the composite bath is probed, allowing it to effectively thermalize between probing attempts. Hole burning would only occur if the energy “consumption” $W\Gamma$ is larger than the energy “supply” $T$

[TABLE]

which does not hold true for a low temperature bath.

We note that we have thus far neglected the “back action” of the fractons on the composite sector. Insofar as $W/T\gg 1$, the back action is strong, and it thus may be tempting to postulate that the composite sector itself could be localized by the coupling to the charge sector, in a form of MBL proximity effect Nandkishore (2015); Hyatt et al. (2017); Nandkishore and Gopalakrishnan (2016). If the bath gets localized, then the relaxation rate will be zero. However, there will inevitably be regions where the density of fractons is low enough for the composite sector to be locally ergodic, and this argument will thus inevitably run into the rare region obstruction to perturbative constructions that have derailed previous attempts to establish translation invariant MBL Schiulaz et al. (2015). We therefore do not pursue this particular line of reasoning further, noting only that Eq. (5) should properly be understood only as an upper bound on the relaxation rate, which could be slower because of back action on the bath.

The situation considered above is one where the three sectors are in equilibrium. We could, in principle, prepare the system out of equilibrium with the bath comprised only of the composites and not the dim-2 $e^{(2)}$’s. In this case, the hopping would proceed at a rate

[TABLE]

since the channel where hopping proceeds through an intermediate boson will be unavailable. Hole burning will again be avoided here as we have a narrow bandwidth bath at high temperature. However, as we show in App. A.1, there will be a rapid equilibration between the composite and bosonic sectors, and upon equilibration, the relaxation will revert to the rate $n_{f}\,e^{-W/T}$.

III.2 Glassy dynamics in the approach to equilibrium

We now consider the non-equilibrium dynamics of type I fracton models prepared in their ground state (i.e., at vanishing energy density) in contact with a low but finite temperature bath comprised of neutral composites (the fully mobile bosons) and dim-2 excitations (the two-dimensional bosons). Let us begin our discussion more generally, by initializing the fractons at a temperature $T_{f}^{(0)}$ such that the density of fractons $n_{f}\sim e^{-W/2T_{f}^{(0)}}$, where $W$ is again the local energy scale in the fracton sector. The neutral composites are prepared at a temperature $T_{c}^{(0)}$ and the $e^{(2)}$’s are prepared at $T_{b}^{(0)}$, with corresponding densities $n_{c}\sim e^{-2W/T_{c}^{(0)}}$ and $n_{b}\sim e^{-W/T_{b}^{(0)}}$. Furthermore, since all three species are gapped here, every species has an exponentially small heat capacity, $C_{i}\sim n_{i}$, for each species $i=b,c,f$. We note that the usual low-temperature $T^{3}$ heat capacity that we expect for three-dimensional bosons holds only when the bosons are gapless and thus does not apply to the neutral composites. Similarly, for the dim-2 bosons the usual $T^{2}$ behaviour which follows from the Stefan-Boltzmann law does not hold here since that applies only to gapless bosons. Hence, here we expect that

[TABLE]

The regime of interest is $T_{i}^{(0)}\ll\Lambda$ where $i=c,b,f$, and we henceforth work in this regime. We further show in the Appendix A that even if we start with $T_{c}\neq T_{b}$ the boson and composite sectors rapidly equilibrate, so we henceforth assume that $T_{c}=T_{b}$. That is, we assume that the bath is itself in equilibrium, and examine the equilibration of the fracton sector with the bath.

III.2.1 Equilibration between Fractons and bath

We now consider the dynamics of the fractons when placed in contact with the thermal sector which contains composites and $e^{(2)}$’s at some initial temperature $T_{b}^{(0)}$. The fractons are prepared at a low temperature, $T_{f}^{(0)}\ll T_{b}^{(0)}$, as we are interested in the dynamics of fractons prepared in their ground state. Further, we will consider only two-body (processes involving more excitations will be further suppressed by density factors) on-shell processes (upto second order in perturbation theory) that lead to an exchange of energy between the bath and the $\mathbb{Z}_{2}$-charged fracton sector. Since we are assuming that the composites and bosons have already equilibrated, the dominant processes through which the bath and the $\mathbb{Z}_{2}$-charged fracton sector exchange energy will involve only the bosons and the fractons. This is because the density of neutral composites $\sim e^{-2W/T_{b}}$ while that of the bosons $\sim e^{-W/T_{b}}$, which dominates in the regime of interest $T_{b}\ll W$. The dominant processes, shown in Fig. 5, are then

Boson + Boson $\leftrightarrow$ 4 Fractons, 2. 2.

Boson + Boson $\leftrightarrow$ Boson + 2 Fractons, 3. 3.

Boson + Fracton $\leftrightarrow$ 3 Fractons.

Two fractons cannot convert into a single boson due to the form of the perturbation, and all other two body processes will involve at least one composite and are hence suppressed compared to the others. Processes where a single composite converts into a pair of bosons which then convert into four fractons or a boson and two fractons occur at third order in perturbation theory and can hence safely be ignored. The above three channels all occur at second order and their rates are controlled by the densities of the excitations involved. For instance, channel 1 (boson + boson) occurs at a rate $\Gamma\sim\frac{\Lambda^{2}}{W}n_{b}^{2}$. Since the bath lends energy $2W$ during this process, $dE_{b}/dt\sim-2W\Gamma=-2\Lambda^{2}n_{b}$. Following this example, we can establish a detailed balance between the thermal and the charged sector,

[TABLE]

Here, the integer coefficients’ magnitude should not be taken seriously, since they are strongly model dependent. However, the conclusions we draw below will not depend on these coefficients, lending our results broader applicability across type I fracton phases.

Since the heat capacities are given by Eq. (III.2), the detailed balance leads to the rate equations

[TABLE]

We can analytically study the equilibration process in two regimes—when the fractons are prepared in their ground space ($T_{f}\sim 0$) and when the system nears equilibration ($T_{b}\sim T_{f}$).

In the case where the fractons are initially prepared at a vanishing energy density, we are in the regime where $T_{f}^{(0)}\ll T_{b}^{(0)}$. Here, we find that (see App. A for details)

[TABLE]

where $b=\exp\left(\frac{2W}{T_{b}^{(0)}}-\frac{W}{2T_{f}^{(0)}}\right)$, and

[TABLE]

such that the fractons display logarithmically slow heating and the bath correspondingly cools logarithmically slowly. This behaviour persists until the fractons are close to equilibration, i.e., over an exponentially long time scale

[TABLE]

Until this time scale, the dominant processes are those of channels 1 and 2, where two bosons combine to pump energy into the fracton sector. Beyond this, however, channel 3 is activated since the fractons are close to equilibration ($T_{f}(t^{*})\sim T_{b}^{(0)}/2$) while the bath’s temperature has only changed slightly ($T_{b}(t^{*})\sim T_{b}^{(0)}$), such that $n_{f}(t^{*})\sim n_{b}(t^{*})$. Hence, we can no longer ignore processes where a boson and a fracton convert into three fractons and the behaviour of the bath is modified at this time scale,

[TABLE]

Thus, the bath first cools logarithmically over an exponentially long time scale, and then rapidly equilibrates with power-law behaviour, since the logarithmic heating of the charged fracton sector establishes a finite density of fractons at $t\sim t^{*}$. This behaviour can be explicitly seen by numerically solving Eq. (III.2.1) (see Fig. 6) and matching the analytic solutions with the numerical curves. As can be seen in Fig. 7, the fractons display logarithmic heating essentially until equilibration while the bath displays logarithmic cooling followed by rapid, power-law cooling, as shown in Fig. 8.

Thus, placing a type I fracton model prepared in its ground space222We note that while we are primarily interested in the limit $T_{f}^{(0)}\ll T_{b}^{(0)}$, due to the form of Eq. (III.2.1), our discussion holds more generally as long as $T_{f}^{(0)}<T_{b}^{(0)}/2$. in contact with a finite temperature heat bath leads to equilibration that exhibits glassy behaviour (i.e. a $\log t$ approach to equilibrium) over an exponentially wide window of intermediate time scales,

[TABLE]

for some positive constant $c$.

At long times, the system eventually equilibrates, $T_{b}=T_{f}$ and close to equilibration, we recover the standard exponential relaxation expected from Newton’s law of cooling,

[TABLE]

where $\delta T=T_{b}-T_{c}$. Thus, over exponentially long time scales, governed by $W/T_{b}^{(0)}$, the relaxation displays logarithmic (glassy) dynamics but close to equilibrium we recover an exponential relaxation with a relaxation rate that is exponentially small in $W/T_{b}^{(0)}$.

III.2.2 Fracton dynamics in an open quantum system

In principle, we can also consider an open quantum system such that the temperature of the composites (and so also the dim-2 bosons) is pinned to that of an external heat bath. In this situation, we are interested in the dynamics of the charged fracton sector which, following the discussion in the previous system, we expect will again demonstrate logarithmically slow heating over exponentially long time scales.

We consider initializing the charged fracton sector at a temperature $T_{f}^{(0)}$, with the neutral composites and dim-2 bosons coupled to an external bath at a temperature $T\gg T_{f}^{(0)}$. Since the fracton sector only exchanges energy with the composites and bosons, most of the discussion from the previous section holds i.e., the dominant processes by which the fracton sector exchanges energy are unchanged. Hence, the rate equation for the fractons is

[TABLE]

where $n_{b}=e^{-W/T}$ is set by the external heat bath and $n_{f}=e^{-W/2T_{f}}$. In the regime of interest, $T_{f}^{(0)}\ll T$, the dynamics are initially governed by channels 1 and 2, such that the fractons display logarithmically slow heating (see App. A for details)

[TABLE]

where $b=\exp\left(\frac{2W}{T}-\frac{W}{2T_{f}^{(0)}}\right)$. This behaviour persists over an exponentially long time scale,

[TABLE]

controlled by the temperature of the heat bath, $T\ll W$. Around $t=t^{*}$, however, channel 3 is activated since a finite density of fractons has been established, and the behaviour of the fractons is modified,

[TABLE]

Hence, the fractons display logarithmically slow heating over an exponentially long time scale, followed by rapid power-law behaviour until they are close to equilibration (see Fig. 9). Near equilibration, $T_{f}\sim T$, and the fractons follow the usual exponential behaviour expected from Newton’s law

[TABLE]

where $\delta T=T-T_{f}$. Thus, even in an open quantum system, fractons display glassy behaviour over exponentially long time scales when prepared in their ground state, where the time scales are controlled by $W/T$ i.e., by the initial temperature of the heat bath. In principle, our discussion applies as long as $T_{f}^{(0)}<T/2$, in which case the width of the $\log(t)$ behaviour will depend on the initial temperature of the fractons, but the situation considered here $T_{f}^{(0)}\ll T$ is of more interest.

While we have focused on the example of the X-Cube model Eq. (1) here, we expect that our results should hold in general for Type I fracton models. In particular, our results here are consistent with the glassy behaviour demonstrated for the related CBLT model Chamon (2005), where a Type I fracton system initialized with an isolated set of fractons and coupled to an external bath was shown to display logarithmic relaxation. While that work considered fractons which were coupled directly to the bath, in contrast we prepare our system with a finite density of fractons which are coupled to neutral composites and dim-2 excitations held at a fixed temperature $T$.

III.3 Role of $m^{(1)}$ excitations

Type I models, such as the X-Cube model Eq. (1), may host additional excitations created in pairs at the ends of Wilson lines of $\sigma_{x}$ operators. These excitations are free to move in one-dimension and hence referred to as $m^{(1)}$ excitations. In order to study the X-Cube model in full generality, we must thus consider the model in the presence of transverse fields, Eq. (2)

[TABLE]

where the term $\lambda\sigma_{x}$ introduces dynamics for the $m^{(1)}$ excitations. Importantly, since the $m^{(1)}$ excitations are created by $\sigma_{x}$ operators, they do not inter-convert with the $e$ excitations (composites, dim-2 bosons, and fractons). Thus, the only coupling between the $e$’s and $m$’s will be through the exchange of energy between the two sectors.

Let us first consider the case where $J=1$ with $\lambda,\Lambda\ll W$. Acting with a single $\sigma_{x}$ operator violates four vertex terms, $B^{(k)}$, such that the gap for creating a single $m^{(1)}$ is equal to $W$, which is the same as that of creating a dim-2 boson. For a system in equilibrium at some temperature $T$, the density of $m^{(1)}$’s will thus be $n_{m}\sim e^{-W/T}$. In addition to the relaxation channels considered in Sec. III.1, the finite density of $m^{(1)}$’s will then provide an additional relaxation channel for the fractons, since they can now borrow energy from the one-dimensional bath of $m^{(1)}$’s with a bandwidth controlled by the hopping rate $\lambda$.

Since the “magnetic” bath can supply a maximum energy $\tilde{\lambda}=\min(\lambda,T)$, the relaxation rate of the fractons due to coupling to this channel will proceed at a rate $\Gamma\sim n_{f}e^{-W/\tilde{\lambda}}$. If $T<\lambda$, then the relaxation will proceed at a rate controlled by the temperature $T$ and $\lambda$ will be rendered irrelevant. On the other hand, if $T>\lambda$, then it is more efficient for the fractons to couple to the “electric” bath, since boson mediated hopping (see Fig. 4) proceeds at a faster rate $\sim n_{f}e^{-W/T}$. Hence, including the dim-1 excitations does not effect the relaxation of fractons, at least in equilibrium, and can safely be ignored insofar as $\lambda$ is weak enough to not destroy the fracton topological order.

As long as $J=1$, the gap for creating an $m^{(1)}$ and an $e^{(2)}$ is the same. Thus, even if we were to consider a system in non-equilibrium, there will be an efficient equilibration between the dim-2 bosons and the dim-1 magnetic excitations, since processes where two $m^{(1)}$’s are destroyed and two $e^{(2)}$’s are simultaneously created will proceed in an on-shell manner, with the rate of equilibration controlled by the relative strength of the transverse fields $\lambda,\Lambda$. Once these sectors have rapidly equilibrated, we can again ignore the $m^{(1)}$’s since the most dominant processes through which the fractons equilibrate will be those considered in Sec. III.2.

If $J\neq 1$, however, then there will exist an imbalance between the electric and magnetic gaps and the equilibration between the $e^{(2)}$’s and $m^{(1)}$’s will no longer proceed in an on-shell manner. Since the case where $J\neq 1$ may also de-stabilise the fracton topological order, depending on the perturbation strengths, we leave the dynamics of the X-Cube model in this case as an open question, to be studied once the phase diagram of this model is better understood Vijay (2017); Ma et al. (2017).

IV Type II Fracton Models

We now turn to type II fracton models, such as Haah’s code Haah (2011). The fundamental difference between type I and II models is the lack of any local string-like operator that allows topological excitations to move in the latter. As we saw in the type I case, pairs of fractons form dimension-1 topological excitations, and as pairs are created at the ends of Wilson-lines, they can move along the Wilson lines without creating any further excitations. However, in type II models, there are no mobile subdimensional excitations.

As a specific example of a type II fracton phase, let us consider Haah’s cubic code model. This is an exactly solvable lattice model defined on a three-dimensional cubic lattice where each site has two spins (or qubits) placed on it. The Hamiltonian is given by the sum over all cubes of two eight-spin interaction terms

[TABLE]

with $G_{c}^{X}$ and $G_{c}^{Z}$ defined in Fig. 10. In analogy with the X-Cube model, Eq. (1), we set $J=1$ here, noting that once we perturb the system, our analysis will hold for $J$ being order unity. The pure Cubic Code model has a topological ground state degeneracy since the ground states cannot be distinguished by any local operator. On a three-torus of length $L$, the ground state degeneracy is $2^{k(L)}$ for some integer $2\leq k(L)\leq 4L$; see (Haah, 2013, Corollary 9.3) for an explicit formula. There exist two kinds of excitations in this model: $e$ type (violation of the $G_{c}^{X}$ term) and $m$ type (violation of the $G_{c}^{Z}$ term). This model has an exact duality between the two types as they are lattice inversions of each other, so it suffices to consider only one sector. We can see that these topological excitations are immobile since a single $\sigma_{z}$ operator on a link creates four cube excitations (see Fig. 11i) and there exists no local string-like operator that can create just a single pair of fractons. Here, single fractons are created at the ends of fractal operators, as shown in Fig. 11ii. Importantly, the “no-strings” rule Haah (2011) implies that any cluster of defects with a non-trivial topological charge must be immobile.

Unlike type I fracton phases, where there is a finite and constant energy barrier for topological charges to move, in type II fracton phases there exists an extensive logarithmic energy barrier preventing topologically charged excitations from diffusing. In particular, as proved in Bravyi and Haah (2011), there exists an energy barrier $\sim c\log(R)$ for creating an isolated fracton with no defects within a distance $R$, for some constant $c$.

Let us consider a type II fracton system in contact with a narrow bandwidth heat bath of composites, with bandwidth $\Lambda\ll W$ and temperature $T$. Let us consider a (non-equilibrium) initial state containing a single fracton. Following our logic in the type I case and replacing the energy scale $W$ in the charged sector by the new energy scale $Wc\log R$, we conclude that to move a distance $R$ will take a time

[TABLE]

Inverting this relation we obtain $R\sim t^{\tilde{\Lambda}/cW}$, which for $\tilde{\Lambda}\ll W$ implies strongly subdiffusive behaviour in a translation invariant three-dimensional model. This subdiffusive behaviour will persist up to the relaxation time, which (as we will shortly show) is super-exponentially long at low temperatures, scaling as $t_{relax}\sim\exp(+c^{\prime}\frac{W^{2}}{T^{2}})$.

Of course, a type II system left in contact with a heat bath will eventually equilibrate to have a non-zero density of fractons—fractons will be created in groups of four (borrowing the energy to do this from the heat bath), and will then be redistributed over the system. However, to achieve an equilibrium distribution with a thermal density of fractons $n_{f}\sim\exp(-W/2T)$, it will be necessary to move fractons over a lengthscale $a\sim\exp(2W/3T)$ (given that fractons can only be created in groups of four). Substitution into the above equation then leads us to conclude that the equilibration time will follow the “super-Arrhenius” scaling

[TABLE]

when $T<\Lambda$, where $c^{\prime}$ is an $O(1)$ numerical constant. This is consistent with a lower bound on relaxation rates derived in Bravyi and Haah (2013). This “super-Arrhenius” scaling is reminiscent of various “near MBL” models e.g. Gopalakrishnan and Nandkishore (2014); Gornyi et al. (2016a) and provides another example of the provocative connections between fracton dynamics and MBL.

We note that classical spin systems, such as the Newman-Moore model Newman and Moore (1999), have been known to display similar phenomenology, in that there exists a logarithmic energy barrier for transporting defects, leading to glassy behaviour. However, unlike type II fracton models, these classical spin models do not have a topologically ordered ground-state subspace and hence lack the robustness against generic local perturbation inherent in quantum fracton models, such as Haah’s code Bravyi et al. (2010); Bravyi and Hastings (2011).

We now compare our results to Siva and Yoshida (2016), who argued that the memory of the initial conditions in the cubic code relaxes for long enough times at finite temperature, with a relaxation time that scaled similarly to Eq. (27). Our broad conclusions are in agreement with Siva and Yoshida (2016). However, our analysis differs in that Siva and Yoshida (2016) examined dynamics towards Gibbs states making use of the Lindblad formalism, which assumes decoherence, whereas we are discussing closed system quantum dynamics, with equilibration implicitly defined in terms of eigenstate thermalization Nandkishore and Huse (2015) and are making no assumptions about decoherence.

To conclude, our analysis reveals that type II models exhibit subdiffusion of fractons on time scales short compared to the relaxation time, and additionally the relaxation time diverges at low temperatures as a super-exponential function of inverse temperature. This is reminiscent of MBL, in that the closed system quantum dynamics of a zero temperature fracton sector coupled to a finite temperature neutral sector could preserve a memory of the initial conditions in local observables over extremely long times that, while finite, are super-exponentially long at low temperature.

V Discussion and Conclusions

We have shown that fracton models naturally demonstrate glassy dynamics. Type I fracton models at finite energy density exhibit a mobility that is exponentially small in the “temperature,” corresponding to a type of “asymptotic many body localization.” Meanwhile, the equilibration of type I fracton models at low temperatures involves a logarithmic approach to equilibrium over an exponentially wide range of time scales—another signature of glassy dynamics. In Type II models, individual fractons exhibit subdiffusive dynamics up to a relaxation time that is super-exponentially long at low temperature, reminiscent of various near-MBL systems, but involving a translation invariant three-dimensional Hamiltonian.

This work has provocative implications for numerous fields. Firstly, it opens up a new direction for the investigation of three-dimensional topological order, suggesting that certain three-dimensional topologically ordered phases (the fracton models) naturally exhibit a rich glassy dynamics more conventionally associated with disordered systems. Can ideas from localization and disordered systems be fruitfully employed to understand three-dimensional topological order? Are there more surprises in store regarding the dynamical behaviour of such models? Given the novelty of these systems, it seems likely that there is more to be discovered.

Additionally, this work may have implications for quantum foundations, in that it identifies a class of quantum systems as being unusually robust to thermalization (in which coupling to a heat bath “observes” the system and “collapses” the wavefunction). This robustness is particularly strong for type II models, for which the time scale for thermalization diverges super-exponentially fast at low temperatures. These models could also have important technological implications, insofar as they display a long lived memory of the initial conditions.

There may also be a connection to the “quantum disentangled liquids” (QDL) introduced in Grover and Fisher (2014). The hallmark of the QDL is that it contains two species of particles, and the many body eigenstates have volume law entanglement entropy, but after performing a projective measurement on the more mobile species of particles, the resulting wavefunction has only area law entanglement entropy. In the fracton models at non-zero temperature, there are indeed two species of particles (fractons and neutral composites), and the neutral composites are in a thermal state, so the eigenstates will indeed have volume law entanglement entropy. However, the fractons themselves can only hop by coupling to the neutral sector, and after performing a projective measurement on the neutral sector, the “fracton only” portion of the wavefunction may well have sub-volume law entanglement entropy. If so, then finite temperature fracton models would provide a three-dimensional realization of a quantum disentangled liquid. Unfortunately, a direct numerical test of this scenario seems difficult, since these models are intrinsically three-dimensional, and numerical tools capable of dealing with three-dimensional glassy many body systems are severely limited.

Thus far we have worked with models with only a $Z_{2}$ charge. However, fracton models should admit of generalizations with $U(1)$ conserved charge Pretko (2016a, b) allowing us to define a charge conductivity. It should then follow, through reasoning analogous to our previous discussion of type II fracton models, that at zero temperature (but with the neutral composites prepared at finite density), that these models should realize a phase that is a thermal conductor but a charge insulator.

Finally, there are the implications for the study of MBL and glass physics. We have introduced a new class of translation invariant models which naturally exhibit glassy dynamics in three dimensions. This new class of models may well provide a new line of attack on important unsolved problems such as MBL in translation invariant systems Kagan and Maksimov (1984); Grover and Fisher (2014); Schiulaz et al. (2015); Yao et al. (2016); Papic et al. (2015); De Roeck and Huveneers (2014a, b); Hickey et al. (2016) and in higher dimensions Chandran et al. (2016); De Roeck and Huveneers (2016). Finally, while the fracton models are defined on lattices, analogous phenomenology can also be obtained in the continuum by making use of higher rank gauge theories. Given the interest in MBL in the continuum Nandkishore (2014); Gornyi et al. (2016b), these too may be worthy of investigation. We leave exploration of these various issues to future work.

Acknowledgements.

We acknowledge useful conversations with Liang Fu, Mike Hermele, Han Ma, and Michael Pretko. RN is supported by the Foundation Questions Institute under grant number FQXi-RFP-1617.

Appendix A Details on Type I Relaxation Dynamics

Here, we discuss in detail the solutions of the relevant differential equations encountered in the main text. All of the relevant equations can be brought to the form

[TABLE]

Since this is a separable equation, this is equivalent to

[TABLE]

Changing variables to $x=\pm\frac{1}{z}$, we get

[TABLE]

Thus, the solution is

[TABLE]

A.1 Equilibration within the bath

First, we consider the dynamics within the bath, ignoring the fractons. As depicted in Fig. 12, a composite can decompose into two bosons while remaining on-shell. Being a first order process (in $\Lambda\ll W$), this occurs at a rate $\sim\Lambda n_{c}$. Similarly, the reverse process where two bosons combine to form a composite occurs at a rate $\sim\Lambda n_{b}^{2}$. To leading order, these are the only two processes which contribute to the equilibration between these two sectors, as all other processes require additional composites/bosons and are suppressed by additional density factors. Considering only the leading order processes, a detailed balance is established between these two sectors,

[TABLE]

Putting in the heat capacities Eq. (III.2), we get

[TABLE]

In the regime where $T_{b}^{(0)}\ll T_{c}^{(0)}$, the bath is initially comprised primarily of neutral composite excitations, and in this limit the above equations are approximately

[TABLE]

For $T_{c}(t)$, with initial condition $T_{c}(0)=T_{c}^{(0)}$, we find

[TABLE]

Since we are working in the limit $T_{c}^{(0)}\ll W$, $T_{c}(t)\sim T_{c}^{(0)}$ for a time $t\sim\frac{W}{\Lambda T_{c}^{(0)}}$. We can hence treat $T_{c}(t)$ as a constant over this time scale such that $T_{b}(t)$ satisfies

[TABLE]

This is equivalent to Eq. (28) with $y=T_{b}/W$ and $a=2\Lambda e^{-\frac{2W}{T_{c}^{(0)}}}$, and hence has the solution

[TABLE]

where $b=e^{\frac{2W}{T_{c}^{(0)}}-\frac{W}{T_{b}^{(0)}}}$. This behaviour is expected to hold until the bosons are close to equilibration $T_{b}\sim T_{c}^{(0)}/2$, i.e., over the window

[TABLE]

where $k$ is some positive constant.

The bosons thus display logarithmic heating over a time scale $0\leq t\lesssim\frac{k}{\Lambda}$, for some positive constant $k$. They reach the equilibrium temperature faster than the neutral composites, which cool down with power-law behaviour over a time scale $0\leq t\lesssim\frac{k^{\prime}W}{\Lambda T_{c}^{(0)}}$ governed by their initial temperature, $T_{c}^{(0)}$. The equilibration process is depicted in Fig. 13, where the approximate analytic solutions developed above are shown to be in good agreement with the numerical solutions for Eq. (A.1). In particular, as illustrated in Fig. 14, the approximate analytic solution Eq. (37) fits the numerical solution for Eq. (A.1) quite well, confirming the logarithmic heating of the bosonic sector.

Close to equilibration, we set $T_{c}=T+\frac{\delta T}{2}$ and $T_{b}=T-\frac{\delta T}{2}$, such that Eq. (A.1) leads to

[TABLE]

in the limit $\delta T\ll T$ and $T=T_{c}^{(0)}/2\ll W$. Thus, close to equilibration the system recovers an exponential relaxation rate since $\delta T(t)=T_{c}(t)-T_{b}(t)$ behaves according to Newton’s law of cooling: $\delta T(t)\sim e^{-\Lambda t}$.

The above discussion demonstrates that while we could, in principle, initialize our system without any bosonic excitations, the presence of a finite density bath of neutral composites will quickly establish a finite density of bosons in equilibrium with the composites (see Fig. 13). Following this discussion, it is easy to show that starting in the opposite limit i.e., with $T_{c}^{(0)}\ll T_{b}^{(0)}$ will result in the same equilibrium configuration, except with the dynamics reversed: the bosons will cool with a power law while the composites will display a $\log(t)$ behaviour. When considering the dynamics of the fracton sector, it is therefore reasonable to assume that the composite and bosonic sectors have already equilibrated. Thus, we make this assumption in the main text and consider a heat bath with both the composites and bosons at an initial temperature $T_{b}^{(0)}$.

A.2 Equilibration between bath and fractons

We now consider equilibration between the charged fractons and the bath (comprised of $e^{(2)}$’s and composites). The equilibration process is governed by Eq. (III.2.1),

[TABLE]

In the regime of interest, $T_{f}^{(0)}\ll T_{b}^{(0)}$, such that the above equations reduce to

[TABLE]

The first equation is equivalent to Eq. (28) with $y=T_{b}/W$ and $a=-3\Lambda^{2}/W$ and hence has the solution

[TABLE]

establishing the logarithmic cooling displayed by the bath. To understand the heating of the fractons, we first note that the bath temperature stays roughly constant, $T_{b}(t)\sim T_{b}^{(0)}$, for an exponentially long period of time $t\sim\frac{W}{3\Lambda^{2}}e^{W/T_{b}^{(0)}}$. We can hence treat $T_{b}$ as a constant, so that the behaviour of the fractons is governed by

[TABLE]

which is again of the form Eq. (28) and has the solution

[TABLE]

where $b=e^{\frac{2W}{T_{b}^{(0)}}-\frac{W}{2T_{f}^{(0)}}}$. Since the bath cools logarithmically over an exponential time scale, this behaviour of the fractons holds until $T_{f}(t)\sim T_{b}^{(0)}/2$ i.e., over an exponentially long time scale

[TABLE]

Thus far, we have only considered processes where two bosons combine to pump energy into the fracton sector, and have neglected processes where a boson and a fracton convert into three fractons. This is no longer accurate once the fracton temperature $T_{f}\sim T_{b}^{(0)}/2$. At this time scale $t\sim\frac{W}{6\Lambda^{2}}e^{W/T_{b}^{(0)}}$, the fracton density $n_{f}\sim n_{b}$, since the temperature of the bath has yet to significantly deviate from its initial temperature $T_{b}^{(0)}$ while the fracton temperature has almost reached its equilibrium value, $T_{b}^{(0)}/2$. Once the fractons are close to equilibration, the behaviour of the bath is hence governed by

[TABLE]

where we are dropping terms of $O(n_{f}^{2})$ and higher. The second term on the r.h.s $\sim n_{f}$ corresponds to channel 3, which is now activated. Treating $n_{f}$ as a constant (since the fracton sector is close to equilibration), we find that the bath’s behaviour is now modified,

[TABLE]

To study the behaviour close to equilibration, we set $T_{b}=T+\frac{\delta T}{2}$ and $T_{f}=T-\frac{\delta T}{2}$, which leads to

[TABLE]

in the limit $T=T_{b}^{(0)}/2\ll W$ and where we are interested in the infinite time i.e, $\delta T\ll T$ behaviour. Here, we recover the standard exponential relaxation expected from Newton’s law of cooling

[TABLE]

A.3 Equilibration between fractons and external heat bath

Finally, we consider the situation where the fractons are prepared at an initial temperature $T_{f}^{(0)}\ll T$, where $T$ is the temperature of an external heat bath, to which the temperatures of the composites and dim-2 bosons are pinned. In this scenario, the behaviour of the fractons is governed by Eq. (19)

[TABLE]

where $n_{b}=e^{-W/T}$ and $n_{f}=e^{-W/2T_{f}}$. As we are interested in the dynamics of fractons prepared in the ground state, in the regime of interest $T_{f}^{(0)}\ll T$ we find that

[TABLE]

Thus, the fractons display logarithmic heating

[TABLE]

where $b=e^{\frac{2W}{T}-\frac{W}{2T_{f}^{(0)}}}$. Note that this is the same behaviour encountered in the previous section. However, once $T_{f}(t)\sim T/2$, we can no longer ignore processes where a boson and a fracton convert into three fractons (channel 3) since $n_{b}^{2}n_{f}\sim n_{b}$ at this point. This is in contrast to the previous section, where the fractons had almost equilibrated by the time channel 3 was activated, leading to a logarithmic heating of fractons essentially over the entire equilibration time scale.

Hence, when the composites and dim-2 excitations are coupled to an external heat bath, the logarithmic behaviour of the fractons persists until $T_{f}(t)\sim T/2$, i.e., over an exponentially long time scale

[TABLE]

controlled by the temperature of the heat bath, $T\ll W$. Beyond this time scale, however, channel 3 is active, and the dynamics of the fractons are governed by

[TABLE]

with the solution

[TABLE]

Thus, the fractons first heat up logarithmically slowly, but once a finite density of fractons is established, they display power law heating until they are close to equilibration. Near equilibration, $T_{f}=T-\delta T$ with $\delta T\ll T$, and the fractons then heat according to

[TABLE]

such that the fractons follow the usual Newton’s law close to equilibration,

[TABLE]