Hilbert space fragmentation produces a"fracton Casimir effect"

Fracton systems exhibit restricted mobility of their excitations due to the presence of higher-order conservation laws. Here we study the time evolution of a one-dimensional fracton system with charge and dipole moment conservation using a random unitary circuit description. Previous work has shown that when the random unitary operators act on four or more sites, an arbitrary initial state eventually thermalizes via a universal subdiffusive dynamics. In contrast, a system evolving under three-site gates fails to thermalize due to strong"fragmentation"of the Hilbert space. Here we show that three-site gate dynamics causes a given initial state to evolve toward a highly nonthermal state on a time scale consistent with Brownian diffusion. Strikingly, the dynamics produces an effective attraction between isolated fractons or between a single fracton and the boundaries of the system, in analogy with the Casimir effect in quantum electrodynamics. We show how this attraction can be understood by exact mapping to a simple classical statistical mechanics problem, which we solve exactly for the case of an initial state with either one or two fractons.


I. INTRODUCTION
A fracton is a kind of excitation in certain quantum systems that exhibits reduced or fractionalized mobility, such that no local operator can move the fracton without producing additional excitations [1][2][3][4][5][6][7][8] (see also Refs. 9 and 10 for reviews). Such reduced mobility can be framed in terms of higher-order conservation laws [6,11,12]. A paradigmatic example, which we focus on in this paper, is that of a one-dimensional system with both charge and dipole moment conservation. In this setting, a single charge can only move in a given direction if a dipole is simultaneously created from the vacuum that points in the opposite direction [6,13]. Such restricted mobility dynamics has led to intense recent interest in fracton systems as candidates for quantum many-body systems that fail to thermalize even at infinite time (e.g., Refs. [14][15][16][17][18][19]. An important test case for fracton systems is the limit of random dynamics (infinite temperature), in which one considers the equal-weight statistical average of all possible unitary evolutions. In a typical quantum many-body system, such random dynamics produces, in the longtime average, a thermal ensemble of all states consistent with the conserved quantities (the "symmetry sector") [20][21][22]. But in fracton systems the restricted mobility can lead to a "fragmentation" of the symmetry sector into disconnected subsectors that are not mutually accessible via unitary evolution [17,18]. Such fragmentation precludes thermalization when it occurs in a "strong" way, such that all subsectors occupy a small volume of the relevant Hilbert space. The question of which conditions produce or preclude thermalization in fracton systems continues to attract significant interest.
A useful method for describing random dynamics is using the random unitary circuit [23][24][25][26][27][28][29], as depicted in Fig. 1 (a), in which the system evolves via discrete gates, each of which is chosen randomly among the set of all possible operators that conserve charge and dipole moment. More precisely, one can consider a spin-1 chain, such that the "charge" at a given site i corresponds to the expec- An illustration of the attraction between fractons (red dots). Fractons move in a given direction by emitting dipoles in the opposite direction, or absorbing dipoles from the same direction. In a system with three-site gates, opposite-facing dipoles cannot pass through each other, and this causes inward-facing dipoles to return more quickly to the fracton that emitted them, biasing the motion of fractons toward each other. (c) An illustration of the mapping between fracton charge (left) and the corresponding height field (right). tation value of the operator S z i , which has eigenstates −1, 0, or +1. In this description the "charge density" at site i is given by S z i (where ... denotes the quantummechanical expectation value) and the total dipole moment of the system P = i x i S z i , where x i is the position of the ith site. Unfortunately, a direct simulation of the quantum evolution of a large system is difficult due to the exponentially large Hilbert space. However, a recently-proposed protocol called automaton dynamics (AD) [30][31][32][33] overcomes this limitation by focusing on operators that take any product state to another product state (multiplied by an overall phase). Such dynamics can be sampled using classical Monte Carlo, enabling the study of large system sizes and long times.
Previous work has shown that for random unitary dynamics with gate size larger than three, the fracton system is thermalized after a long enough time [17,[34][35][36][37]. Under such thermalizing dynamics, the expectation value of any observable (such as the charge density) evolves to that of a thermal ensemble consistent with the fixed charge and dipole moment. The approach to thermalization is described by a universal subdiffusive dynamics, for which a wave vector q is associated with a relaxation time scale τ ∝ 1/q 4 [34][35][36][37][38] (as compared to Brownian diffusion, which has τ ∝ 1/q 2 ). On the other hand, in a circuit with three-site gates, the system fails to thermalize even at infinite time [16,17]. That is, the system retains a memory of the initial state even after infinite time.
In this paper we examine such non-thermalizing dynamics, focusing particularly on unitary circuits with three-site gates. We restrict our attention to the evolution of a spin-1 chain from initial states with a small number of fracton charges. We show that both the final state and the time scale associated with approaching the final state are governed by very different rules than in the thermalizing case. One of our most striking results (first hinted at by numerical results presented in Ref. 16) is an effect reminiscent of the Casimir effect, in which random fluctuations of the "dipole field" produce an effective attraction between two initially separated fracton charges.
One could anticipate this attraction, in a qualitative sense, by imagining an initial state with two fractons in an otherwise empty system [as in Fig. 2(c)]. Each fracton can move (say, to the right) only by emitting a left-facing dipole to its left or absorbing a right-facing dipole from its right. In a common heuristic description (as, e.g., in Ref. 16), one imagines a dipole as a separate kind of particle, which diffuses freely through the system until it encounters a fracton and is absorbed. Importantly, however, under three-site gate dynamics two opposite-facing dipoles cannot pass through each other, and consequently the two fractons effectively reflect each other's emitted dipoles. This reflection causes the dipoles emitted inward to return more quickly than the dipoles emitted outward, leading to a bias in the random motion of the fractons that pulls them together [as depicted in Fig. 2(b)]. Below we make this heuristic picture more precise, and show how the attraction and its associated time scale can be described using an exact mapping to a simple classical problem involving the sliding of "blocks" in an area-preserving height field [depicted in Fig. 2(c)].
We emphasize that the attraction between fractons that we discuss in this paper is qualitatively different, and of different origin, than the "emergent gravity" between fractons that is discussed in Refs. [13,39,40]. We are considering random dynamics, which has no associated Hamiltonian or conserved energy, while these previous references discuss systems characterized by a well-defined Hamiltonian. The major point of our paper is to show that even completely random, non-Hamiltonian dynamics produces an effective attraction between fractons, but it requires the fragmentation of the Hilbert space that is associated with limited operator size.
The reminder of this paper is organized as follows. We begin in Sec. II by considering a random circuit with four-site gates, and showing that the dynamics leads to a thermalized state with a fracton charge density profile that is consistent with a simple maximum-entropy derivation. Section III considers three-site gate dynamics, and presents numeric results for the evolution of initial states having either one or two fractons. We present results for the final state and the time scale associated with its approach. In Sec. IV, we propose an exact mapping to a classical problem of randomly sliding blocks, and we use it to use it to explain our numeric results. We conclude in Sec. V with a summary and discussion of potential future work.

II. SIMULATION METHOD AND THERMALIZATION WITH FOUR-SITE GATES
The system we consider is a one-dimensional spin-1 chain with L sites. The basis states of the system can be written as strings of S z values, with each character in the string denoting the z-component of spin at the corresponding site, which takes one of the values +, 0, or −. The total fracton charge is defined by Q tot = i S z i , and the total dipole moment is P tot = i x i S z i , where x i is the coordinate of site i. We define our units and coordinate system such that x i = i, with the index i running from 1 to L. Each gate in the quantum circuit conserves both the charge and the dipole moment, so that both quantities are constant over time.
Recent work has shown that, under the dynamics of a random unitary circuit with four-site gates, the spin-1 chain thermalizes after a long enough time [17,[34][35][36][37]. Consequently, the entropy of a given observable is maximized within the symmetry sector. In fact, even under four-site gate dynamics, a given symmetry sector is fragmented into exponentially many subsectors [17,18]. But in the thermodynamic limit one of these subsectors becomes dominant, occupying nearly all of the volume of the Hilbert space of the symmetry sector. This "weak fragmentation" allows the dynamics to recover its ergodicity.
In this section we explicitly check the ergodicity of four-site gate dynamics, and we derive a simple relation for the charge density profile S z i (x i ) based on a maximum-entropy argument. We confirm this theoretical result using two different, complementary numerical simulations. First, we directly construct the maximum entropy result by summing over all basis states within the symmetry sector -we refer to this procedure as the "maximum-entropy simulation". Second, we simulate the four-site dynamics using the AD method. We find that the corresponding numerical results for S z i agree both with each other and with the theory.
The details of our AD simulation method are as follows. We start with a well-defined initial state |ψ(0) , which is a product state and represents a specific configuration of charge. We then apply a random ordering of unitary gates on the initial state, with each gate operating on a randomly chosen set of three nearest-neighboring spins. Each gate is chosen randomly from the set of operations that conserve both the charge and dipole moment (the full set of these operations is enumerated below for three-site gate dynamics, and for four-site dynamics it is summarized in Table I of Ref. 16). Under AD, the resulting state remains a product state after the operation of the gate [30,32]. A single time step is defined such that during one time step each spin has, on average, been affected by one gate (i.e., L/n gates constitute a time step for dynamics with n-site gates). A single realization j of this protocol produces a product state |ψ(t) j , where t denotes the number of time steps. To compute the expectation value of some physical operatorÂ, we average the value of ψ(t)| jÂ |ψ(t) j over many independent realizations j. Figure 2 shows the value of S z i as a function of the site index i for an example system with size L = 14, total charge Q tot = 2, and total dipole moment P tot = 7. The red circles correspond to the maximum-entropy simulation, while orange squares and blue crosses are taken from the AD simulation with different initial states. For the AD simulation, the data corresponds to the state of the system after 45000 time steps and is averaged over 5000 independent realizations. The initial state was chosen to have S z j = 1 at two specific sites j, with all other sites i having S z i = 0. From Fig. 2, we see that choosing j = 3, 4 and j = 2, 5 yield numerically indistinguishable results. The close equivalence of the simulation results implies that the system is thermalized after a sufficiently long time, and all information from the initial state is lost.
We can understand the curve S z (x i ) quantitatively using a simple entropy maximization argument. Consider that each site i has some probability p i (s) of having the S z value s = +, 0, or −, with s p i (s) = 1. In the limit where there are many basis states in a given symmetry sector, we can take the probabilities p i (s) to be independent for different sites i, so that the probability of a given string {s 1 , s 2 , ..., s L } is given by i p i (s i ). The corresponding Shannon entropy of the probability distribution is For a given symmetry sector, there are two constraints on the probabilities p i (s i ): We can extremize the value of H using the method of Lagrange multipliers, which gives Here, λ Q and λ P are Lagrange multipliers whose values correspond to the solutions of the nonlinear equations In the limit |Q tot | L and |P tot | L 2 , one can linearize these equations to arrive at In the same limit, the corresponding value of S z i = p i (+) − p i (−) is given by the simple linear equation This analytical result is plotted as the black line in Fig. 2 using the parameters from the simulations, Q tot = 2 and P tot = 7.

III. NONTHERMALIZING DYNAMICS WITH THREE-SITE GATES
As mentioned in the introduction, the dynamics of a random circuit with three-site gates is very different from that of a circuit with four-site gates. In the three-site gate case, the system fails to thermalize due to strong fragmentation of the symmetry sector [17], and the longtime behavior of observables cannot be anticipated using maximum-entropy arguments. Instead, different initial states within the same symmetry sector may evolve to produce qualitatively different values of a given observable. In this section we consider the dynamics of the charge density profile S z i under three-site gate dynamics, starting from initial states with only one or two fractons. We present numerical results here, and characterize the effective attraction between fractons and between fractons and the boundary of the system. Section IV provides an analytical description of these phenomenona via an exact mapping.
Under three-site gate dynamics, the set of operations is very limited by the symmetry constraints. Table I gives the full set. Only eight of the 3 3 possible three-site basis strings admit a nontrivial operation that conserves charge and dipole moment. In the heuristic language of the introduction, the top line of Table I may be viewed as single fracton hopping left (or right) while emitting a right-facing (left-facing) dipole. The bottom two lines can be viewed as the free translation of a dipole.

A. Single fracton initial state
A simple starting point for examining the dynamics of nonthermalizing states is to consider initial states for which there is only a single fracton charge (which, for definiteness, we take to be positive) located at site index i 0 (position x 0 ). In other words, the initial state is The AD simulation method [34] allows us to extract the charge density profile S z i as a function of time. The result is shown in Fig. 3 for the case where i 0 corresponds to the center of the system. As time increases, the delta-function peak of charge spreads outward, as one would expect. However, rather than spreading to uniformly fill the system, as one would FIG. 3. The evolution of the charge density starting from a single fracton peak. The data corresponds to AD simulations with three-site gates for a system with size L = 51. At t = 0, we have a single fracton peak in the middle of the system. As time progresses, the peak gradually spreads outward, before eventually accumulating at the system boundaries.
naively anticipate from the maximum entropy state, at late times the fracton charge accumulates at the boundaries of the system. After a very long time, the fracton charge is completely concentrated into two delta-function peaks, one at either boundary. The weight of the two delta-function peaks is such that the charge and dipole moment from the initial state are preserved [i.e., the right edge of the system has charge (x 0 − 1)/(L − 1) and the left edge has charge (L − x 0 )/(L − 1)].
In order to characterize the time scale associated with the approach to the final state, we define the quantity which describes the effective width of the fracton peak as a function of time. We consider the time scale τ such that r(τ ) = L/4, given an initial state with x 0 = L/2. Figure 4 shows the value of this time scale as a function of system size, along with a power law fit. The fitted exponent (τ ∝ L 2.03 ) suggests a quadratic dependence of the dynamical time scale on the system size, which stands in contrast to the universal dynamics τ ∼ L 4 of thermalizing fracton systems [34][35][36][37][38] .

B. Attraction between fractons
We now consider the case where the initial state consists of two fractons: |ψ(0) = |+ i1 ⊗|+ i2 ⊗ i =i1,i2 |0 i , where i 1 and i 2 are the indices of the two initial fracton positions. Figure 5 shows numerical results for the initial-time and late-time charge density profiles for an initial state with i 1 = 16 and i 2 = 36 and system size L = 51, as given for an initial state with two fractons that evolves under threesite-gate dynamics. The effective attraction between fractons is manifest in a persistent peak of charge density that forms at the midpoint between the initial fracton positions after a long time (red curve). In this example the system size L = 51 and the two initial fracton positions are i1 = 16 and i2 = 36 respectively (blue curve). by an AD simulation. As in the single-fracton case, the boundary of the system exhibits delta-function peaks of charge density. More striking, however, is the appearance of a persistent, localized peak of charge density at the midpoint between the two initial fracton positions. This peak persists even in the limit of infinite time, and has a width significantly smaller than the distance between the two initial fracton positions.
Since the dynamics in our system is explicitly re-versible, this effective attraction between the two fracton peaks can be said to have an entropic origin. That is, within the subsector that includes the initial state, basis states tend to have positive charge near the midpoint between the two initial fracton peaks. An intuitive rationale for this peak is similar to the one given in the introduction: a single fracton moves by emitting dipoles, but under three-site gate dynamics these dipoles cannot pass through each other or through another fracton charge. In this sense each fracton serves like a "boundary condition" for the dipoles emitted by other fractons, in loose analogy with the fluctuations of the electromagnetic field that drive the Casimir effect, and the total configurational entropy of the emitted dipoles is maximized when the two fractons are pushed together. Below we provide a more precise way of describing this effect, which does not rely on an artificial separation of the system's + and − charges into "fractons" and "dipoles." In order to study the time scale associated with the attraction between the two fracton peaks, we define the quantity r Starting from an initial state where the two fracton peaks are equidistant from the center of the system, this quantity represents the average distance between charge at x < L/2 and x > L/2. In the definition of r (t) we neglect the charge peaks at the two boundaries. We define a time scale τ based on the time at which the two fracton peaks approach each other to half their initial separation: r (τ ) = ∆/2, where ∆ is the initial distance between the two fractons. Figure 6 shows the behavior of this time scale as a function of the distance ∆, with the ratio ∆/L kept constant. As in the single fracton case, the time scale τ increases quadratically with the length scale ∆, suggesting that the dynamics is governed by Brownian diffusion rather than by subdiffusion.

IV. AN EXACT MAPPING SOLUTION
In the previous section we showed, numerically, two dynamical phenomena associated with non-thermalizing dynamics that are starkly different from the thermalizing case. First, we showed that fractons exhibit an effective attraction both to each other and to the boundaries of the system. Second, we showed that the time scale associated with the evolution of the fracton charge density is diffusive, τ ∝ L 2 , rather than subdiffusive (τ ∝ L 4 ) as in the thermalizing case. In this section we provide a derivation for these behaviors by mapping to a simple problem of classical statistical mechanics.
The key idea is that the fracton charge as a function of position can be mapped to a height field, which evolves FIG. 6. The time scale τ associated with attraction between two fractons under three-site-gate dynamics is plotted as a function of their initial separation ∆. Blue points correspond to AD simulations with fixed ratio ∆/L = 1/2. The red curve is a power-law best fit, which gives an exponent τ ∝ ∆ 1.97 . The error bars for each point are smaller than the symbol size.
according to simple dynamical rules. For a charge distribution S z i (t), we define the "height" so that h 0 = 0 and h L is equal to the total charge Q tot . In this mapping, the charge at site i corresponds to the discrete derivative Importantly, the total area under the height field is constant as a function of time. Thus, the fracton dynamics corresponds to area-preserving deformations of the height field h i (t). We note that the same height field mapping is used in Ref. 37 to study thermalizing dynamics, and to derive a generalized hydrodynamic relation. The entropy of height fields with fixed area has also been studied in the mathematical physics literature (e.g., Refs. 41 and 42). These previous results, however, correspond to the ergodic situation, where the full set of height fields (with fixed area and fixed end points) can be explored. As we now show, under three-site-gate dynamics only a restricted set of height fields is explored. We mention also that Ref. 19 uses a different, but equivalent classical mapping to label the subsectors of three-site-gate dynamics.
In the language of the height field, the three-site gate operations listed in Table I correspond to simple operations in which a single "block" with unit height and width is slid either to the right or left. This correspondence is outlined in Table II. In this language, the height field can be viewed as a collection of blocks, stacked on top of each other to form some profile h i (t). Under the action of the random circuit these blocks are moved around according to the following simple rules: • A three-site gate corresponds to a single block sliding either to the left or to the right by one unit. Blocks may not slide through each other.
• Any move that would cause the height field to change by two units at one position i is prohibited.
We now show that this simple description of sliding blocks permits an exact solution for both the singlefracton and two-fracton problems.

A. Single fracton block mapping
We first consider the case where the initial state contains only a single fracton at position p. In this situation the height field at time zero satisfies This mapping is illustrated by Fig. 7, with the blue area corresponding to L − p contiguous blocks. At the initial time, the only allowed move is for the leftmost block to slide by one unit to the left. Over time, vacancies open up in the chain of blocks, allowing others to move, and the distribution of blocks becomes more uniform. This process of diffusion with excluded volume is described mathematically by the so-called simple exclusion process [43]. In the continuous limit (L, t 1), the physics is controlled by the usual one-dimensional diffusion equa- tion [43] ∂ρ(x, t) ∂t = a ∂ 2 ρ(x, t) ∂t 2 .
Here ρ(x, t) is the average block density and a (the diffusion constant) is an order-1 parameter determined by the details of the circuit. Although an exact analytical solution for ρ(x, t) is mathematically complicated due to the boundary conditions in our problem, Eq. (16) gives an immediate explanation for the diffusion-like time scale τ ∝ L 2 shown in Fig. 5. This mapping also gives a simple explanation for the delta-function peaks of charge density that appear at the boundaries of the system after a long time. In the block description, it is clear that after a very long time the block density becomes spatially uniform, ρ i = (L − p)/(L − 1). Thus, the height field adopts an average value h i = (L − p)/(L − 1) for all i in the interval 1, ..., L − 1. On the other hand, h 0 = 0 and h L = 1 are fixed, so that taking the discrete derivative [Eq. (13)] gives This expression is consistent with the numeric results shown in Fig. 3, and one can check that it preserves the total charge Q tot = 1 and dipole moment P tot = p.
It is worth noting that this simple reasoning can similarly be applied to any initial state containing an alternating pattern of + and − charges (with any arrangement of 0's between them). Hamiltonian dynamics of such states was considered in detail in Ref. 19. In the language of the height field, such states correspond to a height field that remains always between 0 and 1 (or between 0 and −1), so that it can be viewed as a singletiered arrangement of some number of blocks. Thus, the final state similarly exhibits only two delta-function peaks of charge, whose heights can be determined simply by counting the number of blocks.

B. Double fracton block mapping
A more complicated case is the initial state with two fractons. In this case mapping to the height field gives two layers of blocks, as shown in Fig. 8. In principle, one can directly apply the dynamical rules listed at the beginning of this section to both layers of blocks in order to understand the dynamics. However, a further mapping makes the dynamics easier to understand. We introduce the modified height field which is defined such that the reference state (h = 0) is a filled first layer of blocks (h = 1). We can then discuss two different kinds of excitations in the modified height field: "particles" (denoted by blue blocks in Fig. 8) that correspond to h = 1 and "holes" (red blocks in Fig. 8) that correspond to h = −1. The dynamical evolution of the system consists of the particles and holes diffusing to fill the empty middle region of the system, like two gases expanding into a region of vacuum between them. One subtlety is that, since the height field may not change by two units at one site, the leftmost particle cannot be immediately adjacent to the rightmost hole. Instead, these two must be separated by at least one space, which we describe by a freely-sliding "piston" (gray block in Fig. 8) that separates the two "gases". Since both gases expand outward via the diffusion equation, the observed time scale τ ∝ ∆ 2 [see Fig. 6] is natural. The final state is also easy to understand: it corresponds to a uniform spatial distribution of blocks (red and blue combined). That is, ρ red (x, t) + ρ blue (x, t) becomes approximately constant as a function of position at late times (with the only caveat being the single piston block). The localized peak in fracton charge S z i is dictated by the statistics of the interface between the two gases. Specifically, since blue particles correspond to positive h and red holes to negative h , the charge density is given by The charge density therefore has a peak at the midpoint between the two gases, where ρ red is declining with increased index i while ρ blue is increasing. The width of this peak is of order √ ∆, which is much smaller than the initial distance ∆ between the two fractons.
A full derivation of the charge density S z i at late times is presented in Appendix A. The resulting analytical curve is plotted in Fig. 9 (blue curve), and it very closely matches the numerical result.

V. CONCLUSION
In this paper we have explored the dynamics of nonthermalizing fracton systems, focusing on the case of a spin-1 chain operated on by a random unitary circuit that conserves charge and dipole moment. In the case where the circuit contains four-site gates, the dynamics is thermalizing. We have shown, in particular, that in this case the distribution of charge density after a long time can be captured by a simple maximum-entropy argument, and there is no need to know the specific initial state. On the other hand, the dynamics with only threesite gates is highly nonthermalizing, such that different initial states within the same symmetry sector evolve to produce strongly different distributions of charge density after a long time.
Most striking among our results is that initial peaks of fracton charge density exhibit an effective attraction toward each other and toward the boundaries of the system. The tightly-localized peaks of charge density which appear at late times are in strong contrast to the thermalizing case, for which early-time peaks of charge density tend to spread uniformly across the system. This difference arises fundamentally from the strong fragmentation of the Hilbert space that occurs in the nonthermalizing case.
We explain the dynamics in the nonthermalizing case using a simple, exact mapping to a classical statistical mechanics problem of sliding blocks in a height field. The mapping gives a natural explanation for both the appearance of localized charge peaks at infinite time and for the diffusion-like time scale τ ∝ L 2 associated with the system's approach to the steady state.
In principle, the mapping we present gives a natural way to describe the dynamics of arbitrary initial states. So far, however, we have only focused on two simple cases, which correspond to the height field having either one or two tiers. Generalizing our approach to a more generic set of initial states may be a promising direction for future work.
Unfortunately, the height field mapping does not imply an obvious generalization to higher dimensions. But the question of whether a similarly simple statistical mechanics description can capture the dynamics of twodimensional or higher-dimensional nonthermalizing fracton systems is an interesting one.
In this Appendix we present the full derivation of the late-time charge distribution for an initial state with two fractons. As mentioned in Sec. IV B, we map the problem onto the problem of two kinds of blocks, red and blue, diffusing to fill the initially empty space between them. The average densities of the two kinds of blocks are given by ρ red (x, t) and ρ blue (x, t), respectively.
Separating the two kinds of blocks is a movable 'piston', which can be viewed as a special boundary between them. The densities of red and blue blocks can be written as Here ξ labels the position of the piston and ρ(x, t; ξ) means the average density for those states with the piston located at ξ at time t. W(ξ, t) is a weight function that represents the probability distribution of ξ at time t.
We expect that after the final state is reached, the weight function no longer changes with time and can be written as W(ξ). Also, since the densities of blocks no longer change, which lead to a simple solution This equation means that a static solution is reached for red and blue blocks with a specific ξ. Recalling that the diffusion equation is satisfied for red and blue blocks respectively, a static solution should be a constant function or a linear function. The latter is forbidden since a nonzero slope means the change of block density on the boundaries. Thus we conclude that for the final state, Here the origin of the x axis is chosen to be the middle of the system. ∆ is the distance between the two initial fracton peaks, and for the initial state, we have chosen the two peaks to be symmetric about the middle of the system. The position positions obeys −∆/2 < ξ < ∆/2.
This condition is consistent with the conservation of the block area.
In order for the piston to jump left or right, the target site must be empty of blue or red blocks. Since the piston has same probability to jump to left or right, P(ξ → ξ + 1) ∝ 1 − ρ blue (ξ + 1, ξ) (A11) P(ξ → ξ − 1) ∝ 1 − ρ red (ξ − 1, ξ) Using the fact S z (x) −dρ red /dx + dρ blue /dx (here the discrete derivative in the height mapping has been approximated with a continuous derivative), we get the charge distribution This solution gives a localized peak of charge density in the center of the system in the limit of infinite time, as shown in Fig. 9. The width of this peak is ∼ √ ∆, much smaller than the initial fracton separation when ∆ 1. For the charge on the boundaries of the system (x = ±L/2), we point out that two additional conditions, are needed to make the block system consistent with the fracton system. Using these equations, the charge on boundary can be obtained (A21) Equation (A18), together with the value for the boundaries from Eq. (A21), is plotted as the blue line in Fig. 9.