Machian fractons, Hamiltonian attractors and non-equilibrium steady states

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.


I. 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 [1] 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 dynamics [2].
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.
Machian dynamics, in turn, dictates a remarkable set of properties for systems of N 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 positionvelocity 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!To characterize the non-linear dynamics of our system, we perform a stability analysis on the solution space.We develop an asymptotically self-consistent separation into fast and slow variables that demonstrates the existence of attractors.Further, we show that partition functions for our systems are generically divergent due to the noncompactness of energy hypersurfaces, consistent with the breakdown of equilibration observed in dynamics.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 this paper, 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.

II. SYMMETRIES AND HAMILTONIANS
We consider N identical non-relativistic point particles in d spatial dimensions.The state of the system is specified by N d coordinates, {x µ j , p µ j } where the Greek superscript indices µ = 1, . . ., d denote the component and the Latin subscript indices j = 1, . . ., 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 µ j → x µ j + α µ and leads to the conservation of the total momentum, P µ = j p µ j .The second is the conservation of the total multipole moment Q µ ℓ ≡ j (x µ j ) ℓ .We will focus on ℓ = 1 for now, when Q µ 1 ≡ D µ denotes the dipole moment.Dual to translations, D µ generates rigid shifts of the momentum coordinates [9] as p µ j → p µ j + β µ .A physically sensible and local Hamiltonian compatible with both symmetries takes the form [9] 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| ≤ l M where the Machian length l M is a microscopic length scale that characterizes dynamics [17].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: (2) Equation ( 2) is continuous and differentiable and takes on the desired limiting form where Θ(x) is the Heaviside function.

III. 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 .The few-body dynamics of Eq. ( 1) for N ≤ 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 → ∞ and volume V = L d → ∞ keeping ρ = N/V fixed.Although we first focus on one dimension for simplicity, we find analogous phenomena in higher dimensions.
We focus on random initial conditions at fixed energy.Particles are distributed uniformly in space, with momenta chosen by a random walk in momentum space, terminated when the desired energy is obtained.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.For low densities ρ < l −1 M , generic random initial conditions lead to locations of particles distributed as a Poisson process with mean nearest-neighbor separation ∼ ρ −1 > l M and results, with high probability, in isolated particles lacking any neighbors within l M .All of the energy resides in relatively rare active groups, which splinter and form multiple steady-state clusters as discussed in [9].

For high densities ρ > 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 for (1) ( FIG. 3. (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 ± ϵ where ϵ = 0.28 is fitted to match the amplitude of the oscillating fast variable x. (c), (d): K 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): ẋ = ( ẋ1 − ẋ2)/2 and 2p = p1 − p2 from the exact trajectory, predicted to match in the adiabatic solution Eq. ( 8).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 ≈ ρ number of particles each, but now spaced at regular intervals of distance ≈ 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).
3. 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 (see Appendix C).
4. These observations are generalized straightforwardly to higher dimensions, as shown in Fig. 1(c)-(e).

IV. BROKEN ERGODICITY AND THE HOHENBERG-MERMIN-WAGNER-COLEMAN THEOREM
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 [20].Evidently the theorems are evaded by breaking ergodicity and hence the assumption of validity of statistical mechanics (or equivalent Euclidean quantum mechanics).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 [21,22] 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.First we provide a conceptual understanding of the above results in terms of a self-consistent treatment at late times.Thereafter, we will turn to a statistical mechanical perspective.

V. 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. 3, 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. 3(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 positionspace dynamics.To see this, let us look at the equations of motion for the three-body Hamiltonian Eq. ( 1), The important observation is that the dynamics in Fig. 3 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 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 (see Appendix A): K eff (X, x) is defined as Now we analyze Eq. ( 8) in the adiabatic approximation (see [23] and Appendix A).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.
Assuming that 3X = 1 + ξ for some 1 > ξ > 0, the potential takes the form shown in Fig. 3.This strongly confines the particle within |x| < |ξ| where the particle oscillates rapidly with amplitude ξ.We now feed in the time-averaged fast solution into the equations for X, P simply by replacing x → ξ.This slow motion is generated by the Hamiltonian Equation ( 7) was studied in [9] (see also Appendix A) and describes the dynamics of a pair of fractons, here to be understood as the cluster centers.At late times, the system in Eq. ( 7) reaches a steady state with Ẋ → 0, P diverging and X taking the smallest value so that K eff (X) → 0. For the form shown in Eq. ( 9), this corresponds to |3X| = 1 + ξ, 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 ∼ ξ.The precise value of |ξ| < 1 and the oscillation frequency depend on the initial conditions.In Fig. 3 we compare the adiabatic solution for X, P (broken lines) with the actual dynamics (solid colored lines).We see that by fitting ξ 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 ẋ = 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 = 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 δE ∼ 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 [24].
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 ξ.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 ξ, we see that various initial configurations for the slow variables X(0), P (0) all lead, at late times to Ẋ → 0 and 3X = ± (1 + ξ).This generalizes to arbitrary numbers of particles.The space of the attractors, which is an unbounded continuous space (see Appendix B) can be classified by the locations and membership of the clusters.

VI. 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, Since H is conveniently quadratic in momenta, we can consider integrating them out to generate a statistical probability for the positions of the particles where, we have expressed the Hamiltonian as 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 . . .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 [25], 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 [25] 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 δ( j p j − P tot ) which imposes momentum conservation.The expression det ′ 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 , . . ., 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 , . . ., 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. 4(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 graphtheoretic interpretation is the same.
At this point, the reader may already have recalled the phenomenon of order-by-disorder (OBD) (see [26][27][28] 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.

VII. 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 To begin, let us note that the Poisson bracket between Q ℓ and the total momentum P = j p j is non-vanishing and symmetry generators satisfy the classical multipole algebra [29], From Jacobi's identity, we see that conservation of Q ℓ and P imposes conservation of all Q 1 , . . ., Q ℓ .
x a → x a + α, For concreteness, let us write down the form of L corresponding to Eq. ( 19) for ℓ = 2 As shown in Fig. 5, clustering and ergodicity breaking is observed in the dynamics of quadrupole-conserving fractons qualitatively similar to dipole conserving ones.Equation ( 19) generates a more complex flavour of Machian dynamics.While motion of Eq. ( 1) requires the presence of at least two proximate particles, Eq. ( 19) requires at least ℓ + 1 particles.We may ask whether an alternative Hamiltonian with fewer interacting particles may be found with the same symmetries.We now present a geometric argument to show that this is not so.The important observation is that when a Hamiltonian is built of local terms, each term should be independently symmetric, which places constraints on the available space to explore.Let us begin with dipole conservation and consider a k-body term that preserves it.That it, any dynamics induced by the local term preserves k a=1 x a = Q 1 (24) This tells us that dynamics occurs along a k − 1 dimensional hypersurface in R k preserving Eq. (24).It is only for k ≥ 2 that the dynamics can be non-trivial.For k = 1 for instance, we get a constraint x 1 = Q 1 with no additional freedom and thus no dynamics.Thus, we recover the fact that dipole-conserving Hamiltonians are built out of at least two-body terms.For systems that conserve dipole and quadrupole moments, a local k-body term needs to satisfy Eq. ( 24) as well as an additional constraint, We see that for Eqs. ( 24) and ( 25) describe a k − 2 dimensional hypersurface in R k where dynamics can occur.However, for k = 2, when a solution exists, it yields a discrete set of points with no freedom for dynamics.Thus, a local term generating non-trivial dynamics occurs for k = 3 i.e. a 3-body term.This generalizes to general ℓ.The local term needs to conserve ℓ conservation laws The only way for the system to have dynamics is if the set of equations in Eq. ( 26) are overdetermined i.e. the number of variables are ≥ ℓ + 1.The motion then occurs on the k − ℓ dimensional hypersurface formed by the intersection of Eq. ( 26).

VIII. 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.In our previous work [9], we showed that strong attractions can qualitatively change the two-particle dynamics.We have not found any similar change for the many-particle problem for generic initial conditions, as the effect of momentum divergence dominates the effect of any interaction at late times, leading to robust clustering properties.Nevertheless, we cannot rule out the existence of interactions that break this picture for fine-tuned, e.g.crystalline, initial conditions, and it would be useful to clarify this.Quantizing our system and then comparing what we find with available results on lattice quantum systems is another natural task.We noted that ergodicity breaking at high particle densities is not observed in quantum lattice systems [18,19,[30][31][32].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 [33][34][35] 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 [36][37][38] and harmonic traps [39] that may dynamically produce the conservation of multipole moments.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.We see that the regions 2 a −2 f and 3 a −3 f are unbounded, where P (x 1 , x 2 , x 3 ) diverges.Dynamically, we see that when the particle starts in 1 or 1 a −1 f , it is most likely to encounter 2 a −2 f which represents the steady states described in the main text where two particles form a cluster with the third out of reach.

Appendix C: Initial condition dependence
In the main text, we presented trajectories in 1 and 2 dimensions for systems starting from configurations of approximately uniform density (and ρ > 1).This dynamically leads to the emergent crystallisation, filling the volume.However this does not occur for some specific initial conditions.For example, starting a large number of particles in a small volume, such that all K(x i − x j ) = 1, does not expand to a uniform crystal.Shown in Fig. 11, an single cluster evolves into only 3 clusters at late times, with a large number of particles remaining in the central cluster.This is a general result for such initial conditions: even if the number of particles is taken to be very large, few clusters will form, as large energy barriers form between the individual clusters.

FIG. 1 .
FIG. 1. (a)-(b): Position (xa) and momentum (pa) fractonic trajectories of a 1d system of 20 particles starting with uniform density ρ = 1.4.(c)-(e): Position trajectories of a 2d system of 64 particles starting from uniform density ρ = 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.

FIG. 2 .
FIG.2.The locality function Kg(x) defined in Eq. (2) for various representative values of g.

FIG. 4 .
FIG. 4. (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 λ = 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.

FIG. 10 .
FIG. 10.The position space divided into different regions marked by the number of clusters.We see that 2a − 2 f and 3a − 3 f which have two and three clusters respectively are unbounded.