SuppMat.pdf

# Machian fractons, Hamiltonian attractors and non-equilibrium steady states

Abhishodh Prakash<EMAIL_ADDRESS>(he/him/his) Rudolf
Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3PU,
United Kingdom  Ylias Sadki<EMAIL_ADDRESS>Rudolf Peierls
Centre for Theoretical Physics, University of Oxford, Oxford OX1 3PU, United
Kingdom  S.L. Sondhi<EMAIL_ADDRESS>Rudolf Peierls Centre
for Theoretical Physics, University of Oxford, Oxford OX1 3PU, United Kingdom

###### Abstract

We study the $N$ fracton problem in classical mechanics, with fractons defined
as point particles that conserve multipole moments up to a given order. We
find that the nonlinear Machian dynamics of the fractons is characterized by
late-time attractors in position-velocity space for all $N$, despite the
absence of attractors in phase space dictated by Liouville’s theorem. These
attractors violate ergodicity and lead to non-equilibrium steady states, which
always break translational symmetry, even in spatial dimensions where the
Hohenberg-Mermin-Wagner-Coleman theorem for equilibrium systems forbids such
breaking. While a full understanding of the many-body nonlinear problem is a
formidable and incomplete task, we provide a conceptual understanding of our
results using an adiabatic approximation for the late-time trajectories and an
analogy with the idea of ‘order-by-disorder’ borrowed from equilibrium
statistical mechanics. Altogether, these fracton systems host a new paradigm
for Hamiltonian dynamics and non-equilibrium many-body physics.

_Introduction:_ The notion of thermal equilibrium and the technology of
statistical mechanics are central to our understanding of macroscopic systems.
The idea of ergodicity bridges the intellectual gap between the unceasing
microscopic evolution of any system and the success of time-independent
statistical averages. If the system dynamics is ergodic, the properties of the
late-time states reached by starting from generic initial conditions should
agree 111We will refer to this as equilibration without the qualifying
“thermal”. and be described by statistical mechanics. For classical
systems—and this is a paper about those—the canonical picture of ergodicity is
that while the precise details of a particular trajectory depend sensitively
on initial conditions, a typical trajectory densely covers all phase space
available to it, consistent with conservation laws, by repeatedly revisiting
the vicinity of any allowed phase space point under dynamics222We are being
purists here. In practice, a system may be “ergodic enough for government
work” and given experimental times, this is not a distinction that can be
tested directly. We are also unaware of a usable definition of a system that
is “ergodic enough”..

The question of deciding whether a given Hamiltonian gives rise to ergodic
dynamics or not has a long and distinguished history. For macroscopic systems
one tends to assume that ergodicity is the norm unless the system is
explicitly integrable and that integrable systems are isolated points in
Hamiltonian space. In this paper, we describe a family of Hamiltonian systems
whose native physics violates this expectation and leads to a breakdown of
equilibration and statistical mechanics. These are systems of fractons which
have been the subject of a large volume of work in recent years in the quantum
mechanical setting [3, 4, 5, 6, 7, 8] but whose classical mechanics has only
recently been introduced and studied by two of us and Goriely [9] for small
numbers of particles. More precisely, we consider “ungauged” fractons, i.e.
particles whose dynamics conserve a consistent set of charge multipoles.
Symmetry and locality dictate that such particles obey “Machian” dynamics,
where their inertial response to forces depends entirely on their proximity to
other particles, unlike Newtonian dynamics, where it depends entirely on a
property (the mass) of the particle alone.

Figure 1:  (a)-(b): Position ($x_{a}$) and momentum ($p_{a}$) fractonic
trajectories of a 1d system of 20 particles starting with uniform density
$\rho=1.4$. (c)-(e): Position trajectories of a 2d system of 64 particles
starting from uniform density $\rho=2.56$. Shaded circles with radius $0.5$
are drawn around each particle for visual clarity. The directions of the
momenta are included as arrows in (e). All trajectories are generated by the
Hamiltonian in Eq. 1 with $g=0.3$. Both the 1d and 2d systems exhibit emergent
crystallization starting from random initial conditions.

Machian dynamics, in turn, dictates a remarkable set of properties for systems
with arbitrary numbers $N$ of particles. First, their motion converges at late
times to attractors. Naively this should be impossible in a Hamiltonian system
obeying Liouville’s theorem, but the attractors are in position-velocity space
instead of in phase space, and the relationship between velocities and momenta
is very different in Machian and Newtonian dynamics. Second, there are many
attractors for large $N$ and so the dynamics does not lead to late-time states
whose properties are governed solely by global conserved quantities. Instead,
we see the emergence of further conserved quantities at asymptotically late
times. Third, late-time states always break translation symmetry even in low
dimensions, where the naive invocation of the Hohenberg-Mermin-Wagner-Coleman
[10, 11, 12] theorem would forbid breaking of this continuous symmetry.
Perhaps most striking (see Fig. 1) is the frequent evolution of high-density
fractons systems into states with crystalline order!

As local Machian dynamics is inevitably nonlinear, a full enumeration of
attractors and their basins of attraction for arbitrary $N$ is a formidable
task beyond our current methods, which rely heavily on the computationally
exact solution of the equations of motion. However, we are able to make
considerable progress in understanding these results using three key ideas.
First, we are able to develop a separation into fast and slow variables on the
approach to a given attractor that becomes asymptotically well justified and
allows us to understand the existence of the attractor self-consistently.
Second, we show that the partition function for our systems is formally
divergent at any $N$ and temperature $T$ due to the non-compactness of energy
hypersurfaces, consistent with the breakdown of equilibration observed in
dynamics. Third, despite this divergence, statistical mechanical reasoning of
the kind used in “order by disorder” (OBD) [13, 14] discussions in ergodic
systems can be adapted to gain insight into the temporal evolution of our non-
ergodic system. Essentially, the divergences stem from zero modes in phase
space whose numbers depend on particle configurations in real space, and the
observed dynamics tends to maximize this number.

In the balance of this letter and the supplementary material, we provide the
technical content of the above assertions. Before we do that, we remark that
much recent work on quantum systems has focused on the breakdown of quantum
ergodicity—most closely in lattice fracton systems at low density in the
phenomenon termed “shattering” of Hilbert space and most famously in the
phenomenon of many body localization [15] in disordered systems. Although our
classical systems are very far from these strongly quantum systems with very
small local Hilbert spaces, it is nonetheless notable that we find analogs of
shattering in our multiple attractor dynamics and of localization-protected
quantum order [16] in the breakdown of translation invariance.

_Symmetries and Hamiltonians:_ We consider $N$ identical non-relativistic
point particles in $d$ spatial dimensions. The state of the system is
specified by $Nd$ coordinates, $\\{x^{\mu}_{j},p^{\mu}_{j}\\}$ where the Greek
superscript indices $\mu=1,\ldots,d$ denote the component and the Latin
subscript indices $j=1,\ldots,N$ denote the particle number. We will be
interested in two classes of symmetries. The first is spatial translation,
which acts on position coordinates as $x^{\mu}_{j}\mapsto
x^{\mu}_{j}+\alpha^{\mu}$ and leads to the conservation of the total momentum,
$P^{\mu}=\sum_{j}p^{\mu}_{j}$. The second is the conservation of the total
multipole moment $Q^{\mu}_{\ell}\equiv\sum_{j}(x^{\mu}_{j})^{\ell}$. We will
focus on $\ell=1$ for now, when $Q^{\mu}_{1}\equiv D^{\mu}$ denotes the dipole
moment. Dual to translations, $D^{\mu}$ generates rigid shifts of the momentum
coordinates [9] as $p^{\mu}_{j}\mapsto p^{\mu}_{j}+\beta^{\mu}$. A physically
sensible and local Hamiltonian compatible with both symmetries takes the form
[9]

$H=\sum_{a<b=1}^{N}\frac{\left(\vec{p}_{a}-\vec{p}_{b}\right)^{2}}{2}K(\left|\vec{x}_{a}-\vec{x}_{b}\right|)+\ldots$
(1)

where $K(x)$ is a positive ‘mobility function’ that imposes locality. The
ellipses in Eq. 1 indicate other local symmetric terms, including conventional
interactions, which we drop in this work for simplicity, as their effects do
not qualitatively modify those we report.

In this work, we require $K(x)$ to have a strictly compact support restricted
to $|x|\leq l_{M}$ 333We discuss at the end of the paper what happens when we
loosen this restriction.. It is useful to pick families of functions that
contain as a limit the indicator function on this interval. In this paper, we
use the following family with Machian length $l_{M}$ set to 1:

$K_{g}(x)=\begin{cases}1&\mbox{ $x^{2}\leq 1-2g$}\\\
1-\frac{1}{2g^{3}}(x^{2}-1+2g)^{3}&\mbox{ $1-2g<x^{2}\leq 1-g$}\\\
-\frac{1}{2g^{3}}(x^{2}-1)^{3}&\mbox{ $1-g<x^{2}<1$}\\\ 0&\mbox{ $x^{2}\geq
1$}.\end{cases}$ (2)

Equation 2 is continuous and differentiable and takes on the desired limiting
form

$\lim_{g\rightarrow 0}K_{g}(x)=\Theta(x+1)-\Theta(x-1).$ (3)

where $\Theta(x)$ is the Heaviside function.

_$N$ particle dynamics:_ It is clear from the form of Eq. 1 that the dynamics
is Machian. $H$ vanishes for isolated, immobile, particles and mobility is
restored only by the proximity of others within a Machian length $l_{M}$. In
this simplest of Hamiltonians, this effect is pairwise additive. The few-body
dynamics of Eq. 1 for $N\leq 6$ was studied in [9] where it was shown that
particles initialized in sufficient proximity generically separate into
multiple clusters. The centers of mass of the clusters become immobile and
behave as asymptotic conserved quantities, while particles within a cluster
with more than one particle exhibit oscillations. The position- velocity space
exhibits attractors in the form of stable fixed points and limit cycles, while
there are no attractors in phase space, in conformity with Liouville’s
theorem.

We now turn to the dynamics generated by Eq. 1 for the finite-density problem
of interest for macroscopic systems, i.e. the limit $N\rightarrow\infty$ and
volume $V=L^{d}\rightarrow\infty$ keeping $\rho=N/V$ fixed. Although we focus
on one dimension for simplicity, we find analogous phenomena in higher
dimensions. We typically explore a set of random initial conditions consistent
with a fixed energy. To this end, we distribute the particles roughly
uniformly in space and choose their momenta by a random walk in momentum
space, terminating when the desired energy density is obtained. The resulting
variation in total $P$ and $D$ is not important—one can always find equivalent
initial conditions by symmetry. We have also explored selected configurations
outside this set, such as “big bang” initial conditions in which all particles
start within $l_{M}$. Subsequently, we numerically solve the Hamilton
equations. For example, the plots in Fig. 1 were generated this way with
$K(x)$ in Eq. 2 with $g=0.3$. Our principal findings are as follows:

  1. 1.

For any particle number $N$ for a fixed $L$ one can always find special
initial configurations where $N-M$ particles are exactly stationary and the
energy is concentrated in an $M$ particle _active_ group whose maximum
potential extent $Ml_{M}$ excludes the locations of the stationary ones. This
class of initial conditions leaves stationary particles as they are, and
active particles evolve as discussed in [9] and splinter, leading to the
formation of multiple steady-state clusters.

  2. 2.

For low densities $\rho<l^{-1}_{M}$, generic random initial conditions lead to
locations of particles distributed as a Poisson process with mean nearest-
neighbor separation $\sim\rho^{-1}>l_{M}$ and results, with high probability,
in isolated particles lacking any neighbors within $l_{M}$. Again, all of the
energy resides in the relatively rare active groups that occur with a finite
probability given by the Poisson distribution, which again evolve as discussed
in [9] and splinter, leading to the formation of multiple steady-state
clusters.

  3. 3.

For high densities $\rho>l^{-1}_{M}$ the mean nearest-neighbor separation is
now less than $l_{M}$. Thus, most particles start off within a large active
group that potentially spans the system, seemingly favoring restoration of
ergodicity, starting with generic initial conditions. Indeed, quantum lattice
fractons [18, 19] exhibit such a restoration of ergodicity. Surprisingly, this
does not happen in our models. Instead, we continue to see ergodicity breaking
and the formation of clusters with $\approx\rho$ number of particles each, but
now spaced at regular intervals of distance $\approx l_{M}$. The distribution
of particles among the clusters fluctuates with initial conditions. The
trajectories of a high-density 1d system are shown in Fig. 1(a),(b).

  4. 4.

With big bang initial conditions at high density, the particles generically do
not go on to occupy all the position space but remain localized within a
finite number of clusters, each with a large number of members [20].

  5. 5.

These observations are generalized straightforwardly to higher dimensions, as
shown in Fig. 1(c)-(e).

Figure 2:  (a,b): Position and momentum trajectories of a 3 particle system
(solid colored lines) against adiabatic solutions for the ‘slow’ variables
$X,P$ (dashed lines) obtained from solving Eq. 7. The $X$ trajectory
accurately reproduces the motion of the cluster centers, whereas the $P$
trajectory deviates from it as errors propagate due to its divergent nature.
Dotted lines indicate $X\pm\epsilon$ where $\epsilon=0.28$ is fitted to match
the amplitude of the oscillating fast variable $x$. (c), (d):
$K_{\mathrm{eff}}(x)$ defined in Eq. 9 plotted as a function of $X$ where it
represents an effective mobility function and $x$ where it represents a
confining potential. (e): $\dot{x}=(\dot{x}_{1}-\dot{x}_{2})/2$ and
$2p=p_{1}-p_{2}$ from the exact trajectory, predicted to match in the
adiabatic solution Eq. 8.

_Broken Ergodicity and the Hohenberg-Mermin-Wagner-Coleman theorem:_ The above
observations can be synthesized into the summary that at any density the state
of the fractons, with global conserved quantities fixed, converges to one of a
large number of attractors, all of which spontaneously break translation
symmetry. This occurs even in $d=1,2$ where the theorem of Hohenberg, Mermin,
Wagner, and Coleman (HMWC) [10, 11, 12] forbids the breaking of continuous
symmetries in classical systems at nonzero energy densities 444Recent work
[42, 43] has shown that the HMWC analysis needs to be modified when applied to
multipole symmetry breaking. As far as we can tell, this does not alter the
conditions for translation symmetry breaking.. Clearly, the theorems are
evaded because they assume the applicability of statistical mechanics (or
equivalent Euclidean quantum mechanics) while the breaking of ergodicity
contradicts that assumption. Interestingly, a similar way around the theorems
was discovered in [16] that involved one of us. There it was shown that for
quantum many-body systems with strong quenched disorder that break ergodicity
and exhibit many-body localization (MBL), discrete symmetries [22, 23] can be
broken in highly excited states with finite energy density even in $d=1$—a
phenomenon termed localization-protected quantum order (LPQO). In our case,
the mechanism leading to broken ergodicity is entirely different, so even a
continuous symmetry can be broken.

We now turn to providing a conceptual understanding of the above results.
Immediately we will explain how dynamical evolution leads to cluster formation
by providing a self-consistent treatment at late times. Thereafter, we will
turn to a statistical mechanical perspective.

_Machian Schismogenesis:_ We now analyze the cluster formation for the N
particle problem, for which the three-particle system is a tractable
microcosm. As seen in Fig. 2, particles that start within a single active
group generically splinter into two clusters. Once this sets in, there is a
seeming barrier to further inter-cluster exchanges. We will now show that this
is intimately related to the dynamics in the momentum space, shown in Fig.
2(d). As clustering sets in, the momenta branch out and evolve in such a way
that momentum differences between particles within a cluster are finite,
whereas across the clusters they diverge with time. The latter generates the
barrier observed in position-space dynamics. To see this, let us look at the
equations of motion for the three-body Hamiltonian Eq. 1,

$\displaystyle\dot{x}_{a}$ $\displaystyle=\sum_{b\neq
a}(p_{a}-p_{b})K(x_{a}-x_{b}),$ (4) $\displaystyle\dot{p}_{a}$
$\displaystyle=-\sum_{b\neq
a}\frac{(p_{a}-p_{b})^{2}}{2}K^{\prime}(x_{a}-x_{b}).$ (5)

The important observation is that the dynamics in Fig. 2 occurs along two time
scales, a slow one by the centers of clusters and a fast one within the two-
particle cluster where particles oscillate with a small amplitude and high
frequency. It is convenient to employ a corresponding decomposition of phase-
space variables to reflect this

$\displaystyle x_{1}$ $\displaystyle=X+x,~{}x_{2}=X-x,~{}x_{3}=-2X,$
$\displaystyle p_{1}$ $\displaystyle=P+p,~{}p_{2}=P-p,~{}P_{3}=-2P.$ (6)

The upper and lower case variables are the slow and fast degrees of freedom,
respectively. The phase-space variables are not all independent since the
total momentum and positions are conserved, and we have fixed both to zero
without loss of generality. When momentum branching sets in, we will show that
assuming $P>>p$ and $|x|<<1$ produce a self-consistent solution where the
equations of motion in Eq. 5 are simplified to [20]

$\displaystyle\dot{X}$ $\displaystyle\approx\frac{3P}{2}K_{\rm
eff}(X,x),~{}\dot{P}\approx-\frac{9P^{2}}{4}\frac{\partial K_{\rm
eff}(X,x)}{\partial x},$ (7) $\displaystyle\dot{x}$ $\displaystyle\approx
2p,~{}\dot{p}\approx-\frac{9P^{2}}{4}\frac{\partial K_{\rm eff}(X,x)}{\partial
x}.$ (8)

$K_{\rm eff}(X,x)$ is defined as

$K_{\rm eff}(X,x)\equiv K(3X+x)+K(3X-x).$ (9)

Now we analyze Eq. 8 in the adiabatic approximation [24, 20]. First, we treat
the slow variables $X,P$ as constants and solve the equations for $x,p$ whose
motion corresponds to an effective Newtonian particle with mass $0.5$ in an
external potential.

$H_{X,P}(x,p)\approx p^{2}+\frac{9P^{2}}{4}K_{\rm eff}(X,x).$ (10)

Assuming that $3X=1+\xi$ for some $1>\xi>0$, the potential takes the form
shown in Fig. 2. This strongly confines the particle within $|x|<|\xi|$ where
the particle oscillates rapidly with amplitude $\xi$. We now feed in the time-
averaged fast solution into the equations for $X,P$ simply by replacing
$x\rightarrow\xi$. This slow motion is generated by the Hamiltonian

$H_{\xi}(X,P)\approx\frac{3P^{2}}{4}K_{\rm eff}(X,\xi).$ (11)

Equation 11 was studied in [9] (see also [20]) and describes the dynamics of a
pair of fractons, here to be understood as the cluster centers. At late times,
the system in Eq. 11 reaches a steady state with $\dot{X}\rightarrow 0$, $P$
diverging and $X$ taking the smallest value so that $K_{\rm eff}(X)\rightarrow
0$. For the form shown in Eq. 9, this corresponds to $|3X|=1+\xi$, which self-
consistently supports the earlier assumption for fast motion. With time, we
see that the solution with the adiabatic approximation is increasingly valid:
(i) the cluster centers freeze out at positions $X$ and $-2X$ while the
particles within a cluster rapidly oscillate with amplitude $\sim\xi$. The
precise value of $|\xi|<1$ and the oscillation frequency depend on the initial
conditions. In Fig. 2 we compare the adiabatic solution for $X,P$ (broken
lines) with the actual dynamics (solid colored lines). We see that by fitting
$\xi$ to the amplitude of the fast oscillation (bounds of dotted lines), we
obtain excellent agreement with the motion of the centers of clusters for late
times. Furthermore, the dynamics within a cluster satisfies the Newtonian
relation $\dot{x}=2p$ as expected from Eq. 8. The solution for the momentum
$P$ deviates substantially from the adiabatic solution, as the errors are
compounded due to its divergent nature. However, from our perspective, the
main quantitative physics is in position space, whereas the momentum-space
behavior is important only qualitatively, which the adiabatic approximation
nicely reproduces.

This calculation can also be distilled into a more intuitive understanding by
tracking how energy is distributed. Once clustering sets in, the energy is
carried mainly by the active pair $1-2$,
$E=\frac{1}{2}(p_{1}-p_{2})^{2}=2p^{2}$ and does not depend on the large
values of $P$. However, when one of these particles, say $1$ approaches $3$,
the energy cost of the two entering each other’s range is $\delta
E\sim\frac{1}{2}(3P+p)^{2}$. Thus, $1$ senses a large energy barrier and is
repelled, while $2$ reverses its motion to conserve the center of mass, and
the story repeats. With time, as $P$ increases, so does this energy barrier to
cluster restructuring, and the particles are confined to their clusters.
Although all _physical_ attributes, such as positions and velocities, are
comparable for all three particles, irreconcilable momentum differences make
cluster identities asymptotically immutable. We term this _Machian
schismogenesis_ , after a similar social phenomenon [25].

The generalization to larger number of clusters and higher dimensions is
straightforward. We postulate that motion can always be decomposed into fast
and slow modes. The slow modes, positions of cluster centers, are adiabatic
invariants, which settle down to maximize separation between them just out of
Machian reach and retain fractonic behavior. The fast modes representing
relative motion within each cluster lose their fractonic character and behave
as regular interacting particles within a strong confining external potential
generated by the cluster centers and their divergent momenta. For higher
dimensions where momenta have a larger space to branch out, this naturally
leads to a nearly regular, close-packing arrangement with small deviations
from regularity given by $\xi$. Clustering also results in alignment of the
direction of momenta within each cluster and is visualized by attaching an
arrow corresponding to the direction of the momentum to each particle in Fig.
1(e).

The various clustering choices are attractors [9] in the position-velocity
space of solutions. To see this, notice that from the above calculation for
three particles, keeping the essential dynamics for the fast coordinates $x,p$
fixed i.e. leading to the same amplitude $\xi$, we see that various initial
configurations for the slow variables $X(0),P(0)$ all lead, at late times to
$\dot{X}\rightarrow 0$ and $3X=\pm\left(1+\xi\right)$. This generalizes to
arbitrary numbers of particles. The space of the attractors, which is an
unbounded continuous space [20] can be classified by the locations and
membership of the clusters.

_Failure and success of statistical mechanics:_ We can study the structure of
phase space explored by the fracton system and how the breaking of ergodicity
occurs from the point of view of statistical mechanics. Let us begin by
writing down the partition function in the canonical prescription for the one-
dimensional Hamiltonian in Eq. 1 with conservation laws imposed,

$\mathcal{Z}=\int\prod_{j=1}^{N}dx_{j}\delta(\sum_{j}x_{j}-X_{\rm tot})\\\
\int\prod_{j=1}^{N}dp_{j}\delta(\sum_{j}p_{j}-P_{\rm tot})e^{-\beta H}.$ (12)

Since $H$ is conveniently quadratic in momenta, we can consider integrating
them out to generate a statistical probability for the positions of the
particles

$P(x_{1},\ldots,x_{N})=\int\prod_{j=1}^{N}dp_{j}\delta(\sum_{j}p_{j}-P_{\rm
tot})e^{-\beta H}\\\ =\sqrt{\frac{(2\pi)^{N}}{\beta\det^{\prime}L}}.$ (13)

where, we have expressed the Hamiltonian as

$H=\sum_{j<k}\frac{(p_{j}-p_{k})^{2}}{2}K(x_{j}-x_{k})=\frac{1}{2}\sum_{j,k}p_{j}L_{jk}p_{k}.$
(14)

The probability distribution depends on the nature of the eigenvalues of $L$
that assumes a nice form if we consider the limiting form of the mobility
function shown in Eq. 3. Now, the system can be given the interpretation of an
undirected simple graph $G$, where the particles $1\ldots N$ label the
vertices of the graph, $V(G)$ and the edges $E(G)$ correspond to pairs $(i,j)$
such that $K(x_{i}-x_{j})=1$. The matrix $L$ in Eq. 14 is the Laplacian of $G$
[26],

$L(G)=D(G)-A(G).$ (15)

where, $D(G)$ is the degree matrix of the graph with only diagonal elements
$D_{ii}$ containing the degree of the vertex $i$ and $A(G)$ is the adjacency
matrix of the graph. The disconnected components of the graph $G$ correspond
to clusters. A well-known result [26] states that the number of connected
components of $G$ equals the dimensionality of the nullspace, that is, the
number of zero eigenvalues of $L$. Note that $L$ always has at least one zero
eigenvalue even when $G$ has a single component, which is eliminated by
$\delta(\sum_{j}p_{j}-P_{\rm tot})$ which imposes momentum conservation. The
expression $\det^{\prime}L$ in Eq. 13 denotes the product of all other
eigenvalues that may or may not be zero. For a connected graph, that is, when
all particles that form a single cluster, $P(x_{1},\ldots,x_{N})$ is finite.
_Any_ other configuration that leads to a graph with multipole components
results in a Laplacian $L$ with zero eigenvalues and thus a divergent Eq. 13.
Hence, statistical mechanics fails which is consistent, morally, with its
breakdown viewed from dynamics.

However, not all infinities are the same and we can extract useful guidance
from statistical mechanics by splitting the $N-1$ dimensional space of
positions $x_{1},\ldots,x_{N}$ with a fixed center of mass into different
sectors depending on the connectivity of the graph and postulating that
sectors with more zero modes will appear more frequently in time as the system
evolves. This works surprisingly well, as shown in Fig. 3(a) for a typical
trajectory. However, one cannot use this reasoning to confidently predict the
final state—else the high density big-bang state would necessarily lead to a
system spanning crystal and we have already noted that it does not. The
reasoning also generalizes to higher dimensions, where, although the phase-
space variables are vector-valued, the graph-theoretic interpretation is the
same.

At this point, the reader may already have recalled the phenomenon of order-
by-disorder (OBD) (see [27, 28, 29] and especially [13, 14]) where geometric
frustration leads to a large manifold of ground states, but entropy coming
from integration in orthogonal directions selects configurations that host the
maximum number of soft modes leading to unexpected order. While the family
resemblance to our rationalization of the selection of attractors is very
compelling, it is important to note that OBD is invoked in ergodic systems and
cannot lead to a violation of the HMWC theorem. In more detail, while the
dynamics absolutely takes advantage of the unbounded energy hypersurfaces in
phase space there is no sense in which it is ergodic on them that would
justify the entropy counting within the traditional ergodic framework.

Figure 3:  (a): Histogram of the eigenvalues of the $L$ matrix defined in Eq.
14 for the same simulation as Fig. 1(a). A large peak at $\lambda=0$ is
observed for late times corresponding to the formation of clusters. (b):
Particle densities binned over the same time windows. At late times, the
density is peaked in each cluster, indicating translation symmetry breaking.

_Higher multipole conservation:_ We now generalize the Hamiltonian in Eq. 1
which is invariant under translations and dipole symmetry to those with
translations and multipole symmetry. We keep to one dimension for simplicity,
where the conserved multipole moment is

$Q_{\ell}\equiv\sum_{j}x_{j}^{\ell}.$ (16)

To begin, let us note that the Poisson bracket between $Q_{\ell}$ and the
total momentum $P=\sum_{j}p_{j}$ is non-vanishing and symmetry generators
satisfy the classical multipole algebra [30],

$\\{Q_{\ell},P\\}=\ell Q_{\ell-1},~{}\\{Q_{\ell},Q_{\ell^{\prime}}\\}=0.$ (17)

From Jacobi’s identity, we see that conservation of $Q_{\ell}$ and $P$ imposes
conservation of all $Q_{1},\ldots,Q_{\ell}$.

$\\{\\{H,Q_{\ell}\\},P\\}+\\{\\{P,H\\},Q_{\ell}\\}+\\{\\{Q_{l},P\\},H\\}=0\\\
\implies\\{H,Q_{\ell-1}\\}=0.$ (18)

A Hamiltonian with all these symmetries is

$H=\sum_{a_{1},\ldots,a_{\ell+1}=1}^{N}L^{2}(a_{1},\ldots,a_{\ell+1})K(x_{a_{1}},\ldots,x_{a_{\ell+1}})$
(19)

where $L(a_{0},\ldots,a_{\ell})$ is an $\ell+1$ body term defined as

$L(a_{0},\ldots,a_{\ell})=\sum_{\alpha_{k}\in
a_{0},\ldots,a_{\ell}}\epsilon_{\alpha_{0}\ldots\alpha_{\ell}}~{}p_{\alpha_{0}}\prod_{k=1}^{\ell}x_{\alpha_{k}}^{\ell-k}.$
(20)

$\epsilon_{\alpha_{0},\ldots,\alpha_{\ell}}$ is the Levi-Civita tensor whose
elements are determined through its total antisymmetry property via a choice
for one of the elements, say $\epsilon_{a_{0}\ldots a_{\ell}}=+1$. This choice
does not matter because $L$ is squared in Eq. 19. $K(x_{0},\ldots,x_{\ell})$
is any translationally invariant term that imposes locality on the $\ell+1$
body term. A suitable form for $K$ is

$K(x_{1},\ldots,x_{k})=\prod_{a<b=1}^{k}K(x_{a}-x_{b}).$ (21)

It is easy to verify that Eq. 19 is invariant under translations and the
symmetries generated by $\\{Q_{1},\ldots Q_{\ell}\\}$ i.e.

$\displaystyle x_{a}$ $\displaystyle\mapsto x_{a}+\alpha,$ (22) $\displaystyle
p_{a}$ $\displaystyle\mapsto p_{a}+\sum_{k=0}^{\ell-1}\beta_{k}x_{a}^{k}.$
(23)

For concreteness, let us write down the form of $L$ corresponding to Eq. 19
for $\ell=2$

$L(a,b,c)=p_{a}(x_{b}-x_{c})+p_{b}(x_{c}-x_{a})+p_{c}(x_{a}-x_{b}).$ (24)

Equation 19 generates a more complex flavour of Machian dynamics. While in Eq.
1 and motion requires the presence of at least two proximate particles, Eq. 19
requires at least $\ell+1$ particles. We do not analyze Eq. 19 further but
mention that numerical analysis shows ergodicity breaking signatures that are
qualitatively the same as with dipole conservation [20].

_In closing:_ We have presented a novel, robust setting for ergodicity
breaking in classical systems where symmetries and locality lead to dynamical
non-equilibrium steady states governed by attractors in position-velocity
space that evade both Liouville and Hohenberg-Mermin-Wagner-Coleman theorems.
At the classical level, the next obvious task is to gauge the system and study
the resulting dynamics. This will require us to work in two- and higher-
dimensions. Quantizing our system and then comparing what we find with
available results on lattice quantum systems is another natural task
555Upcoming work. We noted that ergodicity breaking at high particle densities
is not observed in quantum lattice systems [32, 33, 34, 19, 18]. On a related
note, the reader might wonder what the effect is of relaxing the strict
compact nature of the mobility function in Eq. 2. A form with exponential
tails was studied in [9] and was shown to produce clustering which we have
numerically checked persists for a large number of particles as well. This is
again in contrast to quantum lattice models, where an equivalent modification
is expected to restore ergodicity. It would be interesting to further explore
the qualitative changes in physics when passing from a lattice to continuum,
consistent with the phenomenon of UV-IR mixing [35, 36, 37] known to occur in
fracton systems. Finally, it would be useful to make connections to realistic
systems where our results can potentially be observed. A promising setting is
the presence of strong tilted fields [38, 39, 40] and harmonic traps [41] that
may dynamically produce the conservation of multipole moments.

_Acknowledgments_ : We thank John Chalker, Siddharth Parameswaran, David
Logan, Sanjay Moudgalaya, Michael Knap, Frank Pollmann, Jonathan Classen-
Howes, Riccardo Sense, Rahul Nandkishore for helpful discussions and Alain
Goriely for collaboration on related work [9]. A.P. was supported by the
European Research Council under the European Union Horizon 2020 Research and
Innovation Programme, Grant Agreement No. 804213-TMCS and the Engineering and
Physical Sciences Research Council, Grant number EP/S020527/1. S.L.S. and Y.S.
were supported by a Leverhulme Trust International Professorship, Grant Number
LIP-202-014. For the purpose of Open Access, the authors have applied a CC BY
public copyright license to any Author Accepted Manuscript version arising
from this submission.

## References

  * Note [1] We will refer to this as equilibration without the qualifying “thermal”.
  * Note [2] We are being purists here. In practice, a system may be “ergodic enough for government work” and given experimental times, this is not a distinction that can be tested directly. We are also unaware of a usable definition of a system that is “ergodic enough”.
  * Chamon [2005] C. Chamon, Quantum glassiness in strongly correlated clean systems: An example of topological overprotection, Phys. Rev. Lett. 94, 040402 (2005).
  * Haah [2011] J. Haah, Local stabilizer codes in three dimensions without string logical operators, Phys. Rev. A 83, 042330 (2011).
  * Vijay  _et al._ [2015] S. Vijay, J. Haah, and L. Fu, A new kind of topological quantum order: A dimensional hierarchy of quasiparticles built from stationary excitations, Phys. Rev. B 92, 235136 (2015).
  * Nandkishore and Hermele [2019] R. M. Nandkishore and M. Hermele, Fractons, Annual Review of Condensed Matter Physics 10, 295 (2019), https://doi.org/10.1146/annurev-conmatphys-031218-013604 .
  * Pretko  _et al._ [2020] M. Pretko, X. Chen, and Y. You, Fracton phases of matter, International Journal of Modern Physics A 35, 2030003 (2020).
  * Gromov and Radzihovsky [2022] A. Gromov and L. Radzihovsky, Fracton matter, (2022), arXiv:2211.05130 [cond-mat.str-el] .
  * Prakash  _et al._ [2023] A. Prakash, A. Goriely, and S. L. Sondhi, Classical non-relativistic fractons, (2023), arXiv:2308.07372 [cond-mat.str-el] .
  * Hohenberg [1967] P. C. Hohenberg, Existence of long-range order in one and two dimensions, Phys. Rev. 158, 383 (1967).
  * Mermin and Wagner [1966] N. D. Mermin and H. Wagner, Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic heisenberg models, Phys. Rev. Lett. 17, 1133 (1966).
  * Coleman [1973] S. Coleman, There are no Goldstone bosons in two dimensions, Communications in Mathematical Physics 31, 259 (1973).
  * Moessner and Chalker [1998] R. Moessner and J. T. Chalker, Properties of a classical spin liquid: The heisenberg pyrochlore antiferromagnet, Phys. Rev. Lett. 80, 2929 (1998).
  * Chalker [2011] J. T. Chalker, Geometrically frustrated antiferromagnets: Statistical mechanics and dynamics, in  _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_, edited by C. Lacroix, P. Mendels, and F. Mila (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011) pp. 3–22.
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annual Review of Condensed Matter Physics 6, 15 (2015), https://doi.org/10.1146/annurev-conmatphys-031214-014726 .
  * Huse  _et al._ [2013] D. A. Huse, R. Nandkishore, V. Oganesyan, A. Pal, and S. L. Sondhi, Localization-protected quantum order, Phys. Rev. B 88, 014206 (2013).
  * Note [3] We discuss at the end of the paper what happens when we loosen this restriction.
  * Pozderac  _et al._ [2023] C. Pozderac, S. Speck, X. Feng, D. A. Huse, and B. Skinner, Exact solution for the filling-induced thermalization transition in a one-dimensional fracton system, Phys. Rev. B 107, 045137 (2023).
  * Morningstar  _et al._ [2020] A. Morningstar, V. Khemani, and D. A. Huse, Kinetically constrained freezing transition in a dipole-conserving system, Phys. Rev. B 101, 214205 (2020).
  * [20] See Supplementary Material for more details .
  * Note [4] Recent work [42, 43] has shown that the HMWC analysis needs to be modified when applied to multipole symmetry breaking. As far as we can tell, this does not alter the conditions for translation symmetry breaking.
  * Potter and Vasseur [2016] A. C. Potter and R. Vasseur, Symmetry constraints on many-body localization, Phys. Rev. B 94, 224206 (2016).
  * Prakash  _et al._ [2017] A. Prakash, S. Ganeshan, L. Fidkowski, and T.-C. Wei, Eigenstate phases with finite on-site non-abelian symmetry, Phys. Rev. B 96, 165136 (2017).
  * Landau and Lifshitz [1982] L. Landau and E. Lifshitz,  _Mechanics: Volume 1_, v. 1 (Elsevier Science, 1982).
  * Bateson [1935] G. Bateson, Culture contact and schismogenesis, Man 35, 178 (1935).
  * Chung [1997] F. R. Chung,  _Spectral graph theory_ , Vol. 92 (American Mathematical Soc., 1997).
  * Villain  _et al._ [1980] J. Villain, R. Bidaux, J.-P. Carton, and R. Conte, Order as an effect of disorder, Journal de Physique 41, 1263 (1980).
  * Shender [1982] E. Shender, Antiferromagnetic garnets with fluctuationally interacting sublattices, Soviet Journal of Experimental and Theoretical Physics 56, 178 (1982).
  * Chalker  _et al._ [1992] J. T. Chalker, P. C. W. Holdsworth, and E. F. Shender, Hidden order in a frustrated system: Properties of the heisenberg kagomé antiferromagnet, Phys. Rev. Lett. 68, 855 (1992).
  * Gromov [2019] A. Gromov, Towards classification of fracton phases: The multipole algebra, Phys. Rev. X 9, 031035 (2019).
  * Note [5] Upcoming work.
  * Pai  _et al._ [2019] S. Pai, M. Pretko, and R. M. Nandkishore, Localization in fractonic random circuits, Phys. Rev. X 9, 021003 (2019).
  * Sala  _et al._ [2020] P. Sala, T. Rakovszky, R. Verresen, M. Knap, and F. Pollmann, Ergodicity breaking arising from hilbert space fragmentation in dipole-conserving hamiltonians, Phys. Rev. X 10, 011047 (2020).
  * Khemani  _et al._ [2020] V. Khemani, M. Hermele, and R. Nandkishore, Localization from hilbert space shattering: From theory to physical realizations, Phys. Rev. B 101, 174204 (2020).
  * Gorantla  _et al._ [2021] P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, Low-energy limit of some exotic lattice theories and uv/ir mixing, Phys. Rev. B 104, 235116 (2021).
  * You and Moessner [2022] Y. You and R. Moessner, Fractonic plaquette-dimer liquid beyond renormalization, Phys. Rev. B 106, 115145 (2022).
  * Minwalla  _et al._ [2000] S. Minwalla, M. V. Raamsdonk, and N. Seiberg, Noncommutative perturbative dynamics, Journal of High Energy Physics 2000, 020 (2000).
  * Scherg  _et al._ [2021] S. Scherg, T. Kohlert, P. Sala, F. Pollmann, B. Hebbe Madhusudhana, I. Bloch, and M. Aidelsburger, Observing non-ergodicity due to kinetic constraints in tilted fermi-hubbard chains, Nature Communications 12, 4490 (2021).
  * Morong  _et al._ [2021] W. Morong, F. Liu, P. Becker, K. Collins, L. Feng, A. Kyprianidis, G. Pagano, T. You, A. Gorshkov, and C. Monroe, Observation of stark many-body localization without disorder, Nature 599, 393 (2021).
  * Guo  _et al._ [2021] Q. Guo, C. Cheng, H. Li, S. Xu, P. Zhang, Z. Wang, C. Song, W. Liu, W. Ren, H. Dong, R. Mondaini, and H. Wang, Stark many-body localization on a superconducting quantum processor, Phys. Rev. Lett. 127, 240502 (2021).
  * Bagchi  _et al._ [2023] D. Bagchi, J. Kethepalli, V. B. Bulchandani, A. Dhar, D. A. Huse, M. Kulkarni, and A. Kundu, Unusual ergodic and chaotic properties of trapped hard rods, (2023), arXiv:2306.11713 [cond-mat.stat-mech] .
  * Kapustin and Spodyneiko [2022] A. Kapustin and L. Spodyneiko, Hohenberg-mermin-wagner-type theorems and dipole symmetry, Phys. Rev. B 106, 245125 (2022).
  * Stahl  _et al._ [2022] C. Stahl, E. Lake, and R. Nandkishore, Spontaneous breaking of multipole symmetries, Phys. Rev. B 105, 155107 (2022).

See pages 1, of SuppMat.pdf See pages 0, of SuppMat.pdf