# Scaling and localization in kinetically constrained diffusion

Jung Hoon Han<EMAIL_ADDRESS>Department of Physics, Sungkyunkwan
University, Suwon 16419, South Korea  Ethan Lake<EMAIL_ADDRESS>Department of
Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139, USA  Sunghan Ro<EMAIL_ADDRESS>Department of Physics, Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139, USA

###### Abstract

We study diffusion in systems of classical particles whose dynamics conserves
the total center of mass. This conservation law leads to several interesting
consequences. In finite systems, it allows for equilibrium distributions that
are exponentially localized near system boundaries. It also yields an unusual
approach to equilibrium, which in $d$ dimensions exhibits scaling with
dynamical exponent $z=4+d$. Similar phenomena occur for dynamics that
conserves higher moments of the density, which we systematically classify
using a family of nonlinear diffusion equations. In the quantum setting,
analogous fermionic systems are shown to form real-space Fermi surfaces, while
bosonic versions display a real-space analog of Bose-Einstein condensation.

##### Introduction:

The phenomenon of diffusion is ubiquitous in physics, capturing the near-
universal tendency of many-body systems to relax at long times towards a
uniform steady state. Diffusion is typically modeled using Fick’s law by the
equation $\partial_{t}\rho=\nabla\cdot(D\nabla\rho)$, with $\rho$ the particle
density and $D$ the diffusion constant. The assumptions going into the
derivation of this equation are often very minimal, and the results are
applicable to a broad range of physical phenomena. It is therefore important
to understand situations in which conventional diffusive behavior breaks down.

An interesting question to ask in this direction is how diffusion is modified
when restrictions are placed on the ways in which particles can move. A
natural way of doing so is by the imposition of conservation laws that
partially constrain particle motion. An example which has attracted much
interest in the quantum dynamics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] and many-body
physics [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]
communities is provided by dynamics which conserves both the total particle
number $N=\int dx\,\rho(x,t)$ and the total dipole moment $Q_{x}=\int
dx\,x\rho(x,t)$, or equivalently the total center of mass $x_{\rm cm}=Q_{x}/N$
(working in 1d for simplicity). The requirement that both $N$ and $Q_{x}$ be
time-independent is well-known [11] to mandate a continuity equation with two
derivatives, viz.

$\displaystyle\partial_{t}\rho(x,t)=-\partial_{x}^{2}J(x,t).$ (1)

An expression for the current explored in recent hydrodynamically-motivated
studies [ledwith2017angular, 6, 14, 15, 27, glodkowski2023hydrodynamics] is

$J=\widetilde{D}\partial_{x}^{2}\rho,$ (2)

with the subdiffusion constant $\widetilde{D}$ having dimensions of
[length]${}^{4}\ $[time]-1. This expression for $J$ then leads to subdiffusive
behavior with dynamical exponent $z=4$.

Figure 1:  (a) Schematic of an elementary dipole-conserving hopping process,
where particles hop in pairs, and move in opposite directions. (b) Evolution
of the standard deviation of the density $\sigma(t)$. In both 1d (circles) and
2d (triangles), $\sigma(t)$ exhibits an excellent scaling collapse with
dynamical exponent $z=4+d$. (c) Scaled density profiles obtained at different
times for $d=1$ dipole-conserving diffusion ($z=5$). (d) Exponentially-
localized equilibrium density profile for 1d dipole-conserving diffusion on
the half-space $x>0$.

In this work, we show that the physics of dipole-conserving diffusion can be
much richer than what the above equation — which only describes small
fluctuations about a nonzero background density — might suggest. By deriving
an explicit lattice master equation for a natural class of dipole-conserving
dynamics (see Fig. 1(a)), we find a current which instead has the nonlinear
form

$J=D[\rho\partial_{x}^{2}\rho-(\partial_{x}\rho)^{2}]=D\rho^{2}\partial_{x}^{2}\ln(\rho),$
(3)

with $D$ having dimensions of [length]${}^{5}\ $[time]-1 111Note that ordinary
diffusion equation can be rewritten as $\partial_{t}\rho=\partial_{x}J$ with
$J=D\rho\partial_{x}\ln(\rho)$..

The extra power of $\rho$ in (3) yields a dynamic exponent of $z=5$, or
$z=d+4$ in $d$ dimensions. A further notable feature of this equation is that
$\rho(x)\propto e^{-x/\ell}$ is a steady-state solution for all $\ell$. As we
will see, the existence of exponentially-localized steady states is a simple
consequence of entropy maximization, indicating that the general form of (3)
is insensitive to the detailed nature of the microscopic dynamics. Both
localization and $z=4+d$ scaling are readily observed in microscopic
simulations, as seen in Fig. 1(b-d).

The analysis leading to (3) is readily generalized to dynamics conserving
higher multipole moments of the density. For quadrupole-conserving dynamics
conserving both $x_{\mathrm{cm}}$ and the standard deviation
$\sigma=(N^{-1}\int dx\,x^{2}\rho(x,t))^{1/2}$, the steady-state distributions
are Gaussians, with the conservation of $\sigma$ preventing the distribution
from relaxing.

The same guiding principles underlying the above results are operative when
the particles in question obey quantum statistics. Multipole-conserving
fermions turn out to form real-space Fermi surfaces, while bosons exhibit a
real-space analog of Bose-Einstein condensation, with a macroscopic number of
bosons condensing onto a single site.

##### Steady-state distribution:

We begin by understanding the origin of the exponentially-localized steady-
state shown in Fig. 1(d). In the absence of any conserved energy, the steady-
state density $\rho^{\rm eq}(x)$ can be obtained by extremizing the entropy
$S=-\int dx\,\rho(x)\ln\rho(x)$ subject to the conservation of $N=\int
dx\,\rho(x)$ and the center of mass $x_{\mathrm{cm}}=N^{-1}\int dx\,x\rho(x)$.
Carrying out the extremization with the use of Lagrange multipliers directly
yields an exponentially-localized density profile. In particular, if the
particles are confined to the half-line $x>0$, the equilibrium density is

$\rho^{\rm eq}(x)=\frac{N}{x_{\mathrm{cm}}}e^{-x/x_{\mathrm{cm}}}.$ (4)

The distribution of dipole-conserving particles is thus equivalent to the
Boltzmann distribution in the canonical ensemble, with $x_{\mathrm{cm}}$
playing the role of an effective temperature controlling the real-space width
of the particle distribution.

It is important to note that this approach of maximizing $S$ tacitly assumes
ergodic dynamics, an assumption which can break down in certain regimes. As an
extreme example, a single isolated particle is immobile due to the dipole
constraint, and its dynamics is completely non-ergodic. More generally, the
dynamics can be non-ergodic at low particle densities when only the strictly
local hopping processes are considered [2, 3, 4]. In what follows we will
avoid this issue by always working at typical particle numbers per site much
greater than unity, far above [29, 30] the critical density for the onset of
non-ergodicity.

##### Derivation of the diffusion equation:

We now consider a microscopic model within which (3) can be explicitly
derived. To proceed, consider a continuous-time lattice gas where a pair of
particles close to each other are randomly chosen and then displaced locally
in a way that conserves $x_{\mathrm{cm}}$. Denoting the number of particles at
site $x$ as $\rho_{x}$, the particles can “diffuse out” as

$(\Delta\rho_{x},\Delta\rho_{x+1},\Delta\rho_{x+2},\Delta\rho_{x+3})=(+1,-1,-1,+1)$
(5)

with a rate taken to be $rdt\rho_{x+1}\rho_{x+2}$ (here $\Delta\rho$ gives the
change in $\rho_{x}$ between $t+dt$ and $t$). Similarly, the particles may
“diffuse in” as

$(\Delta\rho_{x},\Delta\rho_{x+1},\Delta\rho_{x+2},\Delta\rho_{x+3})=(-1,+1,+1,-1)$
(6)

with a rate $rdt\rho_{x}\rho_{x+3}$ (these rates can be readily checked to
satisfy detailed balance [31]). The fact that these rates are density-
dependent is crucial for obtaining the nonlinear diffusion equation to follow.
This should be contrasted with a distinct class of models inspired by quantum
circuits [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], where the dynamics is driven by
randomly implementing dipole-conserving hops at a density-independent rate.

Using the rates computed above, we obtain a master equation for the evolution
of $\rho_{x}$ of the form 222See also Ref. [ledwith2017angular] for a similar
calculation in the context of transport in 2d Fermi liquids.

$\displaystyle\partial_{t}\rho_{x}$
$\displaystyle=r(\rho_{x-1}\rho_{x+2}+\rho_{x-2}\rho_{x+1}+\rho_{x+1}\rho_{x+2}+\rho_{x-2}\rho_{x-1})$
$\displaystyle-r\rho_{x}(\rho_{x-1}+\rho_{x+1}+\rho_{x-3}+\rho_{x+3}).$ (7)

We now perform a derivative expansion of the above equation, yielding

$\partial_{t}\rho=-D\partial_{x}^{2}\left[\rho\partial_{x}^{2}\rho-(\partial_{x}\rho)^{2}+(\cdots)\right],$
(8)

where we have defined $D\equiv 4ra^{4}$ with $a$ the lattice spacing, and
where $(\cdots)$ indicates terms with higher orders of derivatives. For
distributions with a characteristic length scale $\ell$ (e.g. $\ell=x_{\rm
cm}$ in (4)), these terms are suppressed in powers of $a/\ell$, and (3) is
then reproduced in the limit $a/\ell\ll 1$. More details on the validity of
this limit, as well as what happens when different types of elementary hopping
processes are considered, can be found in the SI [31]. Note that the nonlinear
dependence of $J$ on $\rho$ is a simple consequence of the fact that dipole-
conserving motion always involves pairs of particles.

The unusual scaling observed in the relaxation towards equilibrium follows
from the nonlinearity of the new diffusion equation (8). Indeed, the fact that
$\rho$ has dimensions of inverse length immediately yields a dynamic exponent
of $z=4+1$, in agreement with the measurement of the standard deviation of
density in Fig. 1(b) and dynamic scaling of the density profiles shown in Fig.
1(c). Nevertheless, small fluctuations about a nonzero background density
scale with $z=4$, in accordance with the hydrodynamic expectation (2). This is
seen by taking $\rho(x,t)=\rho^{\rm eq}(x)+\delta\rho(x,t)$, linearizing in
$\delta\rho$, and dropping all derivatives of $\rho^{\rm eq}(x)$, under which
(3) reduces to

$\partial_{t}\delta\rho\approx-D\rho^{\rm eq}(x)\partial_{x}^{4}\delta\rho,$
(9)

with $\widetilde{D}\equiv D\rho^{\rm eq}(x)$ constituting an effective
(sub)diffusion constant. We conclude that the density relaxes with different
dynamical exponents at different length scales.

We briefly note two extensions of this analysis. The first is to incorporate
the effect of noise into the deterministic equation (3). An application of
standard techniques [32, 33] to the present dipole-conserving case gives the
Langevin equation [31] (the structure of which can also be inferred on
symmetry grounds alone [27])

$\partial_{t}\rho(x,t)=-D\partial_{x}^{2}(\rho^{2}\partial_{x}^{2}\ln(\rho))+\sqrt{2D}\partial_{x}^{2}\eta(x,t),$
(10)

where the noise field $\eta$ satisfies

$\displaystyle\langle\eta(x,t)\eta(x^{\prime},t^{\prime})\rangle=\delta(t-t^{\prime})\delta(x-x^{\prime})\rho^{2}(x,t).$
(11)

Linearizing the Langevin equation in $\delta\rho$ allows us to derive the
structure factor [31]

$S(k)\sim\langle\delta\rho(k,t)\delta\rho(-k,t)\rangle=\rho^{\rm eq},$ (12)

which is the same result as for conventional diffusion.

The second extension is the generalization of (3) to $d$ dimensions, where the
dipole-constrained diffusion equation reads
$\partial_{t}\rho=\partial_{a}\partial_{b}J^{ab}$, with the current

$J^{ab}=-D\rho^{2}\partial^{a}\partial^{b}\ln(\rho).$ (13)

The resulting dynamical exponent is then

$z=4+d.$ (14)

We verify this numerically by repeating the scaling analysis for a system of
dipole-conserving particles on a 2d square lattice, with the results of Fig.
1(b) indeed demonstrating excellent $z=6$ scaling.

##### Quadrupole conserving diffusion:

Figure 2: (a) A schematic illustration of a 1d quadrupole conserving diffusion
process: two particles move two sites in one direction, while another particle
moves three sites in the opposite direction. (b) The time evolution of the
density field of particles following the quadrupole-conserving diffusion
process. The density is initially given as two widely separated Gaussian
packets, which merge over time (black $\to$ red $\to$ blue) into a single
Gaussian profile, with the standard deviation $\sigma$ remaining conserved.

We now move on to considering dynamics that conserves the quadrupole moment
$Q_{xx}=\sum_{x}x^{2}\rho_{x}$ (related to the standard derivation as
$\sigma=\sqrt{Q_{xx}/N}$) in addition to the center of mass. Any quadrupole-
conserving processes must involve the coordinated motion of at least three
particles [31], with the simplest allowed three-particle process given by
$(\Delta\rho_{x},\Delta\rho_{x+1},\cdots,\Delta\rho_{x+4})=(\pm 1,\mp 2,0,\pm
2,\mp 1)$, as depicted in Fig. 2(a). By writing down a master equation via
manipulations similar to those that led to (Derivation of the diffusion
equation:), we obtain the quadrupole-conserving diffusion equation

$\partial_{t}\rho=\partial_{x}^{3}J,\qquad
J=D\rho^{3}\partial_{x}^{3}\ln(\rho).$ (15)

Maximizing the entropy subject to the conservation of both dipole and
quadrupole moments yields Gaussian steady states of the form

$\rho^{\rm eq}(x)=\frac{N}{\sqrt{{2\pi}\sigma^{2}}}e^{-(x-x_{\rm
cm})^{2}/2\sigma^{2}},$ (16)

which are also static zero-current solutions to (15). This thus predicts the
scenario whereby a generic initial distribution eventually “congeals” into a
Gaussian, and then remains that way for all further times. The simulation
results shown in Fig. 2(b) indeed confirm this expectation.

Due to the presence of the nonzero conserved length scale $\sigma$, no single-
parameter scaling ansatz $\rho(x^{z}/t)$ with finite $z$ is possible in the
present case. Small deviations about $\rho^{\rm eq}$ will however relax with
$z=6$, as seen by linearizing (15).

##### General multipole moments:

We further generalize the discussion to a $d$-dimensional system that
conserves $\sum_{\mathbf{r}}p_{n}(\mathbf{r})\rho_{\mathbf{r}}$ for all
degree-$n$ polynomials $p_{n}(\mathbf{r})$. Maximizing the entropy yields
equilibrium states formed as exponentials of degree-$n$ polynomials in the
coordinates, which leads to generalizations of (1), (3) and (15) as

$\partial_{t}\rho=-\partial_{A}J^{A},\qquad
J^{A}=(-1)^{n+1}D\rho^{m}\partial^{A}\ln\rho,$ (17)

where $A=\\{a_{1},\dots,a_{n+1}\\}$ is a composite index and
$\partial_{A}\equiv\partial_{a_{1}}\cdots\partial_{a_{n+1}}$. The minimum
value for $m$ in the definition of the current $J^{A}$ (which determines the
mobility $D\rho^{m}$) is equal to the minimum number of particles needed to
form a microscopic $n$-pole conserving hopping process. In [31] we show that
$m=n+1$ for $n\leq 10$ and $n=12$; determining whether or not $m=n+1$ for all
$n$ turns out to be equivalent to a longstanding open problem in number theory
[34]. As a further consistency check, the two relations in (17) imply that the
entropy $S=-\int\rho\ln\rho$ is monotonically increasing with time, as they
yield

$\displaystyle\frac{dS}{dt}=D\int
dx\rho^{m}\left(\partial^{A}[\ln\rho+1]\right)^{2}\geq 0.$ (18)

##### Quantum mutipolar diffusion:

So far, we have treated the objects undergoing diffusion as classical
(distinguishable) particles. We now generalize by letting the particles obey
quantum statistics. This is particularly important in light of the fact that a
natural experimental realization of multipole-conserving dynamics is found in
quantum particles in strongly tilted optical lattices [8, 35, 36].

We begin by generalizing the expression (3) obtained for dipole-conserving
dynamics. This can be done by suitably modifying the lattice master equation.
In the classical case, we took the rate for a pair at locations $x_{1},x_{2}$
to hop to $x_{3},x_{4}$ to be proportional to $\rho_{x_{1}}\rho_{x_{2}}$. To
generalize, we simply modify this to
$\rho_{x_{1}}\rho_{x_{2}}(1+\zeta\rho_{x_{3}})(1+\zeta\rho_{x_{4}})$, where
$\zeta=+1,-1$ for bosons and fermions, respectively. These factors account for
the “statistical interactions” (attractive or repulsive) experienced by
indistinguishable particles, and follow from demanding that the dynamics
satisfy detailed balance [31].

With these modifications to the hopping rates, the same approach used
previously produces [31]

$\displaystyle J$
$\displaystyle=D(\rho\partial_{x}^{2}\rho-(\partial_{x}\rho)^{2}+\zeta(\rho^{2}\partial_{x}^{2}\rho-2\rho(\partial_{x}\rho)^{2}))$
(19)
$\displaystyle=D\rho^{2}(1+\zeta\rho)^{2}\partial_{x}^{2}\ln\left(\frac{1}{1/\rho+\zeta}\right).$

From the second line, it is easy to see that the zero-current steady states
are real-space Bose-Einstein ($\zeta=1$) or Fermi-Dirac ($\zeta=-1$)
distributions:

$\rho^{\rm eq}(x)=[e^{(x-\mu)/T}-\zeta]^{-1},$ (20)

as could be expected on entropy-maximization grounds (here $T,\mu$ are
determined by $N,Q^{x}$). Note that (19) reduces to (3) in the (classical)
low-density limit ($\zeta=0$). However, for bosons at high densities, $J$
becomes cubic in $\rho$, giving a crossover from $z=5$ to $z=6$ dynamics as
the density is increased.

We now point out a few curious features of the steady-state distributions that
arise in this case. For illustration, first consider fermions in a 1d box
$x\in[-L/2,L/2]$. At half-filling $N=L/2$ one finds $\mu=0$, with $T$ fixed by
the center of mass: $T=\infty$ corresponds to $x_{\rm cm}=0$, while
$T=0^{\pm}$ gives $x_{\rm cm}=\mp L/4$ (here ${\rm sgn}(T)=-{\rm sgn}(x_{\rm
cm})$). The fermions thus form a Fermi surface in real space, with the Fermi
“radius” roughly fixed by $N$, and the “temperature” roughly fixed by $x_{\rm
cm}$ 333 Since the dynamics freezes out below a critical (hopping-range-
dependent) density, we expect that the density will follow a Fermi-Dirac
distribution only in regions where $\rho^{\rm eq}(x)\gtrsim 1/\ell_{max}$,
where $\ell_{max}$ is the range of the longest-ranged hopping process [30,
29]..

Now consider $N$ dipole-conserving bosons in the positive-coordinate space
$\mathbb{P}^{d}\equiv[0,\infty)^{d}$, and suppose for simplicity that all
components of the dipole moment $Q^{a}$ are equal. The new physics afforded by
bosonic statistics is the possibility of forming a Bose-Einstein condensate
(BEC). The BECs that form in the present case are distinguished by the fact
that the ordering occurs in real space, with a macroscopic number of bosons
condensing onto a single lattice site.

The transition “temperature” $T_{*}$ at which a real-space BEC forms is
determined by the point at which the chemical potential vanishes, calculated
by solving

$N=\int_{\mathbb{P}^{d}}\frac{d^{d}r}{e^{\sum_{a}r_{a}/T_{*}}-1}$ (21)

for $T_{*}$. For $d=1$ this integral has a short-distance divergence, which
precludes ordering by fixing $T_{*}\rightarrow 0$ (cf. the long-distance
divergences that destroy conventional 1d BECs). When $d>1$, a BEC can form at
$T_{*}\propto N^{1/d}$ [31]. By calculating the dipole moment $Q^{x}$ at
$T=T_{*}$, we find that a BEC forms when $N>N_{*}$, with

$N_{*}\propto(Q^{x}/a)^{\frac{d}{d+1}},$ (22)

where we have restored the lattice constant $a$. To understand this, imagine
increasing the number of particles in the system while keeping $Q^{x}$ fixed.
For $N<N_{*}$, the system will smoothly adjust its overall density profile as
particles are added. When $N>N_{*}$, all of the particles subsequently added
to the system will join the condensate at $\mathbf{r}={\mathbf{0}}$, where
they do not contribute to the dipole moment. In this regime the system will
thus possess a condensate at the origin, together with a surrounding “normal
fluid” localized nearby.

Bosons with higher moment conservation behave similarly. As an example,
quadrupole-conserving bosons in infinite space have $\rho^{\rm
eq}(r)=[e^{(r^{2}-\mu)/T}-1]^{-1}$ (assuming rotational symmetry), rendering
the problem equivalent to the nonrelativistic ideal Bose gas after the
interchange of coordinates and momenta. A short-distance divergence prevents a
BEC from forming when $d<3$ [31]. For $d\geq 3$, the critical boson number is
determined in terms of the quadrupole moment $Q^{xx}$ as
$N_{*}\propto(Q^{xx}/a^{2})^{\frac{d}{d+2}}$. Understanding other physical
properties of these unusual BECs, especially the nature of their condensation
transitions, are questions that we leave to future work.

##### Discussion:

We have demonstrated that multipole conservation laws can radically alter
diffusive processes. Even in the absence of interactions, these conservation
laws lead to exponentially localized steady-state distributions, with the
approach to equilibrium governed by unusual dynamical exponents (e.g. $z=4+d$
in the dipole-conserving case).

It would be of great interest to pursue experiments realizing multipole-
conserving dynamics in the lab. A natural platform in this regard is found in
strongly tilted optical lattices [2, 4, 8, 35, 36], where the lattice tilt
enforces emergent dipole conservation over a long prethermal time scale. The
breakdown of diffusion due to ergodicity breaking can be avoided by working at
sufficiently high densities, with the envisioned experiment rather similar to
the one performed in Ref. [36].

There are many natural extensions of our work. Instead of full multipole
moment conservation, one could examine non-maximal multipole groups [38, 21]
in $d>1$, where the dynamics conserves
$\sum_{\mathbf{r}}P(\mathbf{r})\rho_{\mathbf{r}}$ for only a subset of
degree-$n$ polynomials $P(\mathbf{r})$, e.g. a 2d system conserving
$P(\mathbf{r})=x^{2},y^{2}$, but not $P(\mathbf{r})=xy$. It may also be
interesting to examine the importance of the lattice geometry, with more
complicated lattices potentially leading to emergent subsystem symmetries á la
the mechanism of Ref. [39].

Another intriguing direction is to break the time-reversal symmetry of the
diffusion process. Such modification would make the system _active_ and is
known to result in nonequilibrium phenomena that cannot be captured by entropy
maximization [40, 41, 42].

##### Acknowledgments:

We are grateful to Soonwon Choi, Johannes Feldmeier, Byungmin Kang, Andy
Lucas, Francisco Machado, Seth Musser, Rahul Nandkishore, Marvin Qi, Brian
Skinner, Charles Stahl, Daniel Swartz, and Hongzheng Zhao for helpful
discussions and feedback. J.H.H. was supported by the National Research
Foundation of Korea(NRF) grant funded by the Korea government(MSIT) (No.
2023R1A2C1002644). He also acknowledges financial support from EPIQS Moore
theory centers at MIT and Harvard, where this work was initiated.

## References

  * Pai  _et al._ [2019] S. Pai, M. Pretko, and R. M. Nandkishore, Localization in fractonic random circuits, Physical Review X 9, 021003 (2019).
  * Khemani  _et al._ [2020] V. Khemani, M. Hermele, and R. Nandkishore, Localization from hilbert space shattering: From theory to physical realizations, Physical Review B 101, 174204 (2020).
  * Rakovszky  _et al._ [2020] T. Rakovszky, P. Sala, R. Verresen, M. Knap, and F. Pollmann, Statistical localization: From strong fragmentation to strong edge modes, Physical Review B 101, 125126 (2020).
  * 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, Physical Review X 10, 011047 (2020).
  * Iaconis  _et al._ [2021] J. Iaconis, A. Lucas, and R. Nandkishore, Multipole conservation laws and subdiffusion in any dimension, Physical Review E 103, 022142 (2021).
  * Feldmeier  _et al._ [2020] J. Feldmeier, P. Sala, G. De Tomasi, F. Pollmann, and M. Knap, Anomalous diffusion in dipole-and higher-moment-conserving systems, Physical Review Letters 125, 245303 (2020).
  * Schulz  _et al._ [2019] M. Schulz, C. Hooley, R. Moessner, and F. Pollmann, Stark many-body localization, Physical review letters 122, 040606 (2019).
  * Scherg  _et al._ [2021] S. Scherg, T. Kohlert, P. Sala, F. Pollmann, B. H. Madhusudhana, I. Bloch, and M. Aidelsburger, Observing non-ergodicity due to kinetic constraints in tilted fermi-hubbard chains, Nature Communications 12, 1 (2021).
  * Feng and Skinner [2022] X. Feng and B. Skinner, Hilbert space fragmentation produces an effective attraction between fractons, Physical Review Research 4, 013053 (2022).
  * Moudgalya  _et al._ [2021] S. Moudgalya, A. Prem, D. A. Huse, and A. Chan, Spectral statistics in constrained many-body quantum chaotic systems, Physical Review Research 3, 023176 (2021).
  * Pretko [2017] M. Pretko, Subdimensional particle structure of higher rank u (1) spin liquids, Physical Review B 95, 115139 (2017).
  * Prem  _et al._ [2018] A. Prem, M. Pretko, and R. M. Nandkishore, Emergent phases of fractonic matter, Physical Review B 97, 085116 (2018).
  * Pretko [2018] M. Pretko, The fracton gauge principle, Physical Review B 98, 115134 (2018).
  * Gromov  _et al._ [2020] A. Gromov, A. Lucas, and R. M. Nandkishore, Fracton hydrodynamics, Physical Review Research 2, 033124 (2020).
  * Glorioso  _et al._ [2022] P. Glorioso, J. Guo, J. F. Rodriguez-Nieva, and A. Lucas, Breakdown of hydrodynamics below four dimensions in a fracton fluid, Nature Physics 18, 912 (2022).
  * Gorantla  _et al._ [2022] P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, Global dipole symmetry, compact lifshitz theory, tensor gauge theory, and fractons, Phys. Rev. B 106, 045112 (2022).
  * Pielawa  _et al._ [2011] S. Pielawa, T. Kitagawa, E. Berg, and S. Sachdev, Correlated phases of bosons in tilted frustrated lattices, Physical Review B 83, 205135 (2011).
  * Sachdev  _et al._ [2002] S. Sachdev, K. Sengupta, and S. Girvin, Mott insulators in strong electric fields, Physical Review B 66, 075128 (2002).
  * Lake  _et al._ [2022a] E. Lake, M. Hermele, and T. Senthil, Dipolar bose-hubbard model, Phys. Rev. B 106, 064511 (2022a).
  * Lake  _et al._ [2022b] E. Lake, H.-Y. Lee, J. H. Han, and T. Senthil, Dipole condensates in tilted bose-hubbard chains, arXiv:2210.02470 (2022b).
  * Stahl  _et al._ [2022] C. Stahl, E. Lake, and R. Nandkishore, Spontaneous breaking of multipole symmetries, Physical Review B 105, 155107 (2022).
  * Zechmann  _et al._ [2022] P. Zechmann, E. Altman, M. Knap, and J. Feldmeier, Fractonic luttinger liquids and supersolids in a constrained bose-hubbard model, arXiv:2210.11072 (2022).
  * Lake and Senthil [2023] E. Lake and T. Senthil, Non-fermi liquids from kinetic constraints in tilted optical lattices, arXiv preprint arXiv:2302.08499 (2023).
  * Anakru and Bi [2023] A. Anakru and Z. Bi, Non-fermi liquids from dipolar symmetry breaking, arXiv preprint arXiv:2304.01181 (2023).
  * Chen  _et al._ [2021] S. A. Chen, J.-K. Yuan, and P. Ye, Fractonic superfluids. ii. condensing subdimensional particles, Physical Review Research 3, 013226 (2021).
  * Yuan  _et al._ [2020] J.-K. Yuan, S. A. Chen, and P. Ye, Fractonic superfluids, Physical Review Research 2, 023267 (2020).
  * Guo  _et al._ [2022] J. Guo, P. Glorioso, and A. Lucas, Fracton hydrodynamics without time-reversal symmetry, Physical Review Letters 129, 150603 (2022).
  * Note [1] Note that ordinary diffusion equation can be rewritten as $\partial_{t}\rho=\partial_{x}J$ with $J=D\rho\partial_{x}\ln(\rho)$.
  * 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, Physical Review B 101, 214205 (2020).
  * [31] See the supplemental material for details.
  * Lefèvre and Biroli [2007] A. Lefèvre and G. Biroli, Dynamics of interacting particle systems: stochastic process and field theory, Journal of Statistical Mechanics: Theory and Experiment 2007, P07024 (2007).
  * Hexner and Levine [2017] D. Hexner and D. Levine, Noise, diffusion, and hyperuniformity, Phys. Rev. Lett. 118, 020601 (2017).
  * Borwein [2002] P. Borwein,  _Computational excursions in analysis and number theory_ (Springer Science & Business Media, 2002).
  * Guardado-Sanchez  _et al._ [2020] E. Guardado-Sanchez, A. Morningstar, B. M. Spar, P. T. Brown, D. A. Huse, and W. S. Bakr, Subdiffusion and heat transport in a tilted two-dimensional fermi-hubbard system, Physical Review X 10, 011042 (2020).
  * Zahn  _et al._ [2022] H. Zahn, V. Singh, M. Kosch, L. Asteria, L. Freystatzky, K. Sengstock, L. Mathey, and C. Weitenberg, Formation of spontaneous density-wave patterns in dc driven lattices, Physical Review X 12, 021014 (2022).
  * Note [2] Since the dynamics freezes out below a critical (hopping-range-dependent) density, we expect that the density will follow a Fermi-Dirac distribution only in regions where $\rho^{\rm eq}(x)\gtrsim 1/\ell_{max}$, where $\ell_{max}$ is the range of the longest-ranged hopping process [30, 29].
  * Gromov [2019] A. Gromov, Towards classification of fracton phases: the multipole algebra, Physical Review X 9, 031035 (2019).
  * Bulmash  _et al._ [2023] D. Bulmash, O. Hart, and R. Nandkishore, Multipole groups and fracton phenomena on arbitrary crystalline lattices, arXiv preprint arXiv:2301.10782 (2023).
  * Soto and Golestanian [2014] R. Soto and R. Golestanian, Run-and-tumble dynamics in a crowded environment: Persistent exclusion process for swimmers, Physical Review E 89, 012706 (2014).
  * Thompson  _et al._ [2011] A. G. Thompson, J. Tailleur, M. E. Cates, and R. A. Blythe, Lattice models of nonequilibrium bacterial dynamics, Journal of Statistical Mechanics: Theory and Experiment 2011, P02029 (2011).
  * Kourbane-Houssene  _et al._ [2018] M. Kourbane-Houssene, C. Erignoux, T. Bodineau, and J. Tailleur, Exact hydrodynamic description of active lattice gases, Physical review letters 120, 268003 (2018).
  * 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, Physical Review B 104, 235116 (2021).
  * Note [3] The nonlinearity measured by $m$ discussed here refers to the largest power of the density appearing on the RHS of the master equation. When the hopping process is such that at most one particle hops from any given site, all terms on the RHS will be of the same order. In other cases (such as the 3-site dipole hopping process discussed in App. II) smaller powers may additionally appear. In this appendix we are interested in determining the minimal value that the largest power on the RHS can be.
  * Macdonald [1998] I. G. Macdonald,  _Symmetric functions and Hall polynomials_ (Oxford university press, 1998).

Scaling and localization in kinetically constrained diffusion:
Supplementary information

Jung Hoon Han1, Ethan Lake2, and Sunghan Ro2

1Department of Physics, Sungkyunkwan University, Suwon 16419, South Korea
2Department of Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts 02139, USA

(Dated: )

### Contents

  * •

Section I: details of calculations relevant for including quantum statistics.

  * •

Section II: discussion of how to pass from the lattice master equation to the
continuum, as well as alternate forms of the master equation which produce
modifications to the currents derived in the main text.

  * •

Section III: analysis of the minimal number of particles needed to participate
in multipole conserving hopping processes.

  * •

Section IV: derivation of the dipolar Langevin equation.

##  I Quantum multipolar diffusion

In this section we provide a few details relating to the incorporation of
quantum statistics. We will begin by deriving a master equation which
correctly encorporates indistinguishability, and see how this can be used to
derive a diffusion equation that holds for any type of statistics (classical,
bosonic, or fermionic). We then elaborate on a few phenomenological details of
the resulting equilibrium distributions.

###  I.1 Master equation and diffusion equation

Our goal in this subsection is to write down an extension of our classical
multipole-conserving diffusion process to the setting where the microscopic
particles are either bosons or fermions. This can be done with a simple
modification of the rates appearing in the master equation, as we now explain.
We will mostly focus on the dipole-conserving case, as the extension to higher
multipole moments is straightforward.

####  I.1.1 detailed balance

First let us review the case where the statistics is classical. Let
$A=\\{\rho_{x}^{A}\\}$ and $B=\\{\rho^{B}_{x}\\}$ be two configurations of
particles related by a dipole hopping process that hops particles at sites
$x_{1},x_{2}$ to sites $x_{3},x_{4}$. Let $P_{A\rightarrow B}$ denote the
probability for a transition between states $A,B$ to occur (and for simplicity
assume that $x_{1}\neq x_{2},x_{3}\neq x_{4}$). The most natural
generalization of conventional diffusion is to let $P_{A\rightarrow B}$ to be
proportional to the number of pairs available to hop, so that

$P_{A\rightarrow B}\propto\rho^{A}_{x_{1}}\rho^{A}_{x_{2}},\qquad
P_{B\rightarrow
A}\propto\rho^{B}_{x_{3}}\rho^{B}_{x_{4}}=(\rho^{A}_{x_{3}}+1)(\rho^{A}_{x_{4}}+1).$
(S1)

In the main text we stated that this rule satisfies detailed balance. To check
this, we need to specify the distribution with respect to which detailed
balance is satisfied. This distribution is simply the usual one, which assigns
equal weights to each microstate, i.e. each configuration of particles (which
for now are treated as being distinguishable). Thus the equilibrium
probability of being in state $A$ is

$P_{A}\propto\frac{N!}{\prod_{x}\rho^{A}_{x}!},$ (S2)

where $N$ is the total number of particles. We then see that detailed balance
is obeyed, as

$\frac{P_{A}}{P_{B}}=\frac{\prod_{x=1,2,3,4}\rho^{B}_{x}!}{\prod_{x=1,2,3,4}\rho^{A}_{x}!}=\frac{(\rho^{A}_{x_{1}}-1)!(\rho^{A}_{x_{2}}-1)!(\rho^{A}_{x_{3}}+1)!(\rho^{A}_{x_{4}}+1)!}{\rho^{A}_{x_{1}}!\rho^{A}_{x_{2}}!\rho^{A}_{x_{3}}!\rho^{A}_{x_{4}}!}=\frac{(\rho^{A}_{x_{3}}+1)(\rho^{A}_{x_{4}}+1)}{\rho^{A}_{x_{1}}\rho^{A}_{x_{2}}}=\frac{P_{B\rightarrow
A}}{P_{A\rightarrow B}}.$ (S3)

We now generalize to the setting of indistinguishable particles. We claim that
the natural generalization yields transition rates

$\displaystyle P^{\zeta}_{A\rightarrow B}$
$\displaystyle\propto\rho^{A}_{x_{1}}\rho^{A}_{x_{2}}(1+\zeta\rho^{A}_{x_{3}})(1+\zeta\rho^{A}_{x_{4}}),$
(S4) $\displaystyle P^{\zeta}_{B\rightarrow A}$
$\displaystyle\propto\rho^{B}_{x_{3}}\rho^{B}_{x_{4}}(1+\zeta\rho^{B}_{x_{1}})(1+\zeta\rho^{B}_{x_{2}})=(\rho^{A}_{x_{3}}+1)(\rho^{A}_{x_{4}}+1)(1+\zeta\rho^{A}_{x_{1}}-\zeta)(1+\zeta\rho^{A}_{x_{2}}-\zeta),$

where $\zeta=+1,-1,0$ for bosons, fermions, and classical particles,
respectively. When $\zeta=0$ we clearly recover the classical case treated
above. When $\zeta=-1$, the additional factors can be understood simply as a
consequence of Pauli blocking; when $\zeta=+1$ these factors are a consequence
of the statistical attraction experienced by bosons. The necessity of
including these factors can be seen by checking detailed balance. For
$\zeta=\pm 1$, $P^{\zeta}_{A\rightarrow B}/P^{\zeta}_{B\rightarrow A}$ is
equal to unity for all allowed (non-Pauli-blocked) processes. This is
consistent with detailed balance because for indistinguishable particles
$P_{A}$ is independent of $A$: the configurations $\\{\rho_{x}\\}$ are now
themselves microstates, and hence are all assigned equal probabilities in
equilibrium.

####  I.1.2 master equation

With the above modifications to the hopping rates, it is now easy to write
down a master equation valid for any value of $\zeta$. For the 4-site
microscopic hopping processes considered in the main text, the full expression
is

$\displaystyle\partial_{t}\rho=r(1+\zeta\rho_{x})\Big{(}\rho_{x-1}\rho_{x+2}(1+\zeta\rho_{x+1})+\rho_{x-2}\rho_{x+1}(1+\zeta\rho_{x-1})$
(S5)
$\displaystyle\qquad+\rho_{x+1}\rho_{x+2}(1+\zeta\rho_{x+3})+\rho_{x-2}\rho_{x-1}(1+\zeta_{x+3})\Big{)}$
$\displaystyle-r\rho_{x}\Big{(}\rho_{x-1}(1+\zeta\rho_{x+1})(1+\zeta\rho_{x-2})+\rho_{x+1}(1+\zeta\rho_{x-1})(1+\zeta\rho_{x+2})$
$\displaystyle\qquad+\rho_{x-3}(1+\zeta\rho_{x-2})(1+\zeta\rho_{x-1})+\rho_{x+3}(1+\zeta\rho_{x+1})(1+\zeta\rho_{x+2})\Big{)}.$

Master equations for processes conserving higher multipole moments are derived
similarly.

####  I.1.3 diffusion equation

We now perform a derivative expansion on the master equation
(LABEL:quantum_master), normalizing the lattice constant to $a=1$.
Unilluminating algebra then yields $\partial_{t}\rho=-\partial_{x}^{2}J$,
where the current is now

$\displaystyle J$
$\displaystyle=D(\rho\partial_{x}^{2}\rho-(\partial_{x}\rho)^{2}+\zeta(\rho^{2}\partial_{x}^{2}\rho-2\rho(\partial_{x}\rho)^{2})$
(S6)
$\displaystyle=D\rho^{2}(1+\zeta\rho)^{2}\partial_{x}^{2}\ln\left(\frac{1}{1/\rho+\zeta}\right)$

where $D=4r$ for the 4-site hopping process considered above. From the
rewriting on the second line, it is easy to derive the steady states

$\rho^{\rm eq}(x)=\frac{1}{e^{\beta(x-\mu)}-\zeta},$ (S7)

which are the Bose-Einstein ($\zeta=1$) and Fermi-Dirac ($\zeta=-1$)
distributions that we expect on entropy-maximization grounds. This also
suggests that the generalization to $n$-pole conserving dynamics can be done
by taking $\partial_{t}=-\partial_{x}^{n+1}J$, where now

$J=(-1)^{n+1}D(\rho+\zeta\rho^{2})^{n+1}\partial_{x}^{n+1}\ln\left(\frac{1}{1/\rho+\zeta}\right).$
(S8)

As a consistency check, we see that regardless of $\zeta$, our diffusion
equation agrees with the $\zeta=0$ result in the low-density (classical)
limit, where the effects of indistinguishability become unimportant. Thus for
e.g. dipole-conserving dynamics, the density relaxes with dynamical exponent
$z=5$ in all cases. For high-density bosons however, things are different:
here the term proportional to $\zeta$ in (S6) dominates. Since this term
contains three powers of $\rho$, the dynamical exponent is instead $z=6$.

###  I.2 Equilibrium physics

In this section we elaborate on the phenomenology of the steady states in the
bosonic and fermionic cases.

####  I.2.1 bosons

As a simple illustrative example, consider $N$ dipole-conserving bosons living
in the space $\mathbb{P}^{d}\equiv[0,\infty)^{d}$. For simplicity, we will
assume that all $d$ components of the dipole moment
$Q^{a}\equiv\int_{\mathbf{r}}r^{a}\rho^{\rm eq}(\mathbf{r})$ are equal. A
point to note is that, in this section, we will be discussing how the
equilibrium distribution changes when the total particle number $N$ is varied.
For this reason it will be more convenient to work at fixed dipole moment,
rather than at fixed center of mass $r_{\rm cm}^{a}=Q^{a}/N$. In this case we
may write the equilibrium distribution as

$\rho^{\rm eq}(\mathbf{r})=\frac{1}{e^{\beta(\sum_{a}r_{a}-\mu)}-1}$ (S9)

with $\beta>0$. The constant $\mu$ — which by positivity of $\rho^{\rm
eq}(\mathbf{r})$ and $\beta$ cannot be positive — is determined by fixing the
overall number of particles to be $N$:

$N=\int_{\mathbb{P}}d^{d}r\frac{1}{e^{\beta(\sum_{a}r_{a}-\mu)}-1}.$ (S10)

$\beta$ on the other hand is determined by fixing the total dipole moment as

$Q^{a}=\int_{\mathbb{P}}d^{d}r\frac{r^{a}}{e^{\beta(\sum_{a}r_{a}-\mu)}-1}.$
(S11)

Since $\mu$ is monotonically increasing with $N$, there is a
($\beta$-dependent) value of $N$ beyond which (S10) cannot be satisfied. For
larger values of $N$, real-space Bose condensation will occur, with the
remaining bosons occupying the “zero-energy” state at $\mathbf{r}=0$. The
critical value $N_{*}$ beyond which a real-space BEC forms is determined by
setting $\mu=0$ at fixed $\beta$ in (S10), yielding

$N_{*}=\beta_{*}^{-d}\int_{\mathbb{P}}d^{d}r\,\frac{1}{e^{\sum_{a}r^{a}}-1}\equiv\beta_{*}^{-d}C(d)$
(S12)

where $\beta_{*}$ will be determined by fixing the dipole moment, and

$C(d)=\begin{dcases}\infty&d=1\\\ \pi^{2}/6&d=2\\\ \zeta(3)&d=3\\\
\frac{\pi^{4}}{90}&d=4\end{dcases}\,\,.$ (S13)

The divergence for $d=1$ — which in the present setting is a short-distance
divergence — prevents a real-space BEC from forming in one dimension.

We can relate $N_{*}$ to the dipole moment through the requirement that

$Q^{a}=\beta_{*}^{-(d+1)}\int_{\mathbb{P}}d^{d}r\,\frac{r^{a}}{e^{\sum_{a}r^{a}}-1}=\beta_{*}^{-(d+1)}C(d+1).$
(S14)

Solving this equation for $\beta_{*}$ and substituting into the result for
$N_{*}$, we find

$N_{*}=\frac{C(d)}{C(d+1)^{\frac{d}{d+1}}}(Q^{a})^{\frac{d}{d+1}}.$ (S15)

Imagine keeping the dipole moment $Q^{a}$ fixed while increasing $N$.
Initially the “temperature” of the Bose-Einstein distribution will decrease
with $N$, with the particles becoming progressively more localized around the
origin. Once $N>N_{*}$, all of the remaining particles will be added at the
origin, where they join the real-space condensate. From the point of view of
the particles away from the origin, particle number conservation is thus
spontaneously broken, leading to a phenomenon one might refer to as “off-
diagonal short-range order”.

Note that since we are fixing the dipole moment — as opposed to the center of
mass — the number of particles at the origin will increase without bound as
$N$ is increased. This can be done at fixed $Q^{a}$ since adding particles to
the origin does not modify the dipole moment (note that this discussion would
be modified if we fixed the total center of mass, rather than the total dipole
moment).

Bosons conserving higher multipole moments can be analyzed in much the same
way. As an example, consider $N$ quadrupole-conserving bosons living in
$\mathbb{R}^{d}$. We will again assume for simplicity that all $d$ diagonal
components of the quadrupole moment
$Q^{aa}=\int_{\mathbf{r}}(r^{a})^{2}\rho^{\rm eq}(\mathbf{r})$ are equal, and
that the off-diagonal components vanish; $Q^{ab}\propto\delta^{a,b}$.
Following the same analysis as in the dipole case, a real-space BEC forms when
the number of particles exceeds the value

$N_{*}=\beta_{*}^{-d/2}\int_{\mathbb{R}^{d}}d^{d}r\,\frac{1}{e^{\sum_{a}(r^{a})^{2}}-1}\equiv\beta_{*}^{-d/2}A(d),$
(S16)

where

$A(d)=\begin{dcases}\infty&d=1\\\ \infty&d=2\\\
\pi^{3/2}\zeta(3/2)&d=3\end{dcases}\,\,.$ (S17)

Here the short-distance divergence which preempts a real-space BEC is stronger
than in the dipole-conserving case on account of the “energy”
$\sum_{a}(r^{a})^{2}$ vanishing more quickly near the origin (note that the
present quadupole-conserving case is formally equivalent to the non-
relativistic Bose gas, except with the roles of position space and momentum
space interchanged). We may again relate $N_{*}$ to the (fixed) quadrupole
moment $Q^{aa}$ by calculating

$Q^{aa}=\beta_{*}^{-(d/2+1)}\int_{\mathbb{R}^{d}}d^{d}r\,\frac{(r^{a})^{2}}{e^{r^{2}}-1}\equiv\beta_{*}^{-(d/2+1)}B(d)$
(S18)

where

$B(d)=\begin{dcases}\frac{\sqrt{\pi}}{2}\zeta(3/2)&d=1\\\
\frac{\pi^{3}}{12}&d=2\\\ \frac{\pi^{3/2}}{2}\zeta(5/2)&d=3\end{dcases}\,\,.$
(S19)

Therefore we may write the critical particle number in terms of $Q^{aa}$ as

$N_{*}=\frac{A(d)}{B(d)^{\frac{d}{d+2}}}(Q^{aa})^{\frac{d}{d+2}}.$ (S20)

The interpretation of this equation is similar to the dipole-conserving case:
as we increase $N$ past $N_{*}$ at fixed $Q^{aa}$, the additional particles
are all added to the real-space condensate at $\mathbf{r}={\mathbf{0}}$.

####  I.2.2 Fermions

Figure S1:  Equilibrium distributions for $N=50$ dipole-conserving fermions
living on the half-space $[0,\infty)$, plotted for various values of the
dipole moment $Q_{x}=\int dx\,x\,\rho^{\rm eq}(x)$. Larger dipole moments
yield larger effective temperatures.

The analysis for fermions is similar, although of course condensation does not
occur. For this reason the physics is similar in all spatial dimensions, and
accordingly we will specify to $d=1$ for simplicity. Consider first $N$
dipole-conserving fermions on $\mathbb{P}=[0,\infty)$. The parameters
$\beta,\mu$ appearing in the equilibrium distribution are determined through

$N=\int_{\mathbb{P}}dr\,\frac{1}{e^{\beta(r-\mu)}+1},\qquad
Q^{a}=\int_{\mathbb{P}}dr\,\frac{r}{e^{\beta(r-\mu)}+1}.$ (S21)

The first equation determines $\mu$ as a function of $\beta$ via

$\beta\mu=\ln(e^{\beta N}-1),$ (S22)

while the second equation implicitly determines $\beta$ through

$Q^{a}=-\beta^{-2}{\rm Li}_{2}(1-e^{\beta N}).$ (S23)

The equilibrium distribution for $N=50$ fermions on $\mathbb{P}$ is shown in
Fig. S1 for different values of the dipole moment $Q_{x}$. The minimal value
of $Q_{x}$ is obtained when all of the fermions pile up around the origin;
this setup gives $Q_{x}=N^{2}/2$. When $Q_{x}$ is increased the effective
“temperature” of $\rho^{\rm eq}(x)$ also increases, as seen in the figure.
Higher multipole moments are analyzed similarly.

##  II continuum limits, scaling, and modified master equations

In this appendix we discuss different types of master equations and the
scaling relations needed to obtain continuum differential equations from them.
We will restrict our attention to one-dimensional dipole-conserving systems
for concreteness; more general settings can be analyzed similarly.

###  II.1 Continuum limits

As a review, first consider how one would derive the conventional diffusion
equation from a microscopic master equation. The simplest setting involves
non-interacting classical particles hopping in 1d. Here one considers a
process such that in each time interval $dt$, any given particle has a
probability $r\,dt$ of hopping one lattice site to the left, and the same
probability of hopping one site to the right. On a site $x$ with $\rho_{x}$
particles, the probabilities that some particle hops to the left or the right
are consequently equal to $r\,dt\rho_{x}$. Summing up the contributions from
particles which hop into and out of a given site $x$, we then have

$\rho_{x}(t+dt)=\rho_{x}(t)+r\,dt(\rho_{x+1}(t)+\rho_{x-1}(t)-2\rho_{x}(t)),$
(S24)

which in the limit $dt\rightarrow 0$ gives the master equation

$\partial_{t}\rho_{x}=r(\rho_{x+1}+\rho_{x-1}-2\rho_{x}).$ (S25)

One scaling limit which produces a continuum differential equation is as
follows. Consider a system containing $L$ lattice sites with lattice spacing
$a$. We then take $a\rightarrow 0$ while sending
$r\rightarrow\infty,L\rightarrow\infty$ in such a way that $ra^{2}$ and the
dimensionful system size $aL$ remain fixed. This then produces the usual
diffusion equation

$\partial_{t}\rho=D\left(\partial_{x}^{2}\rho+\frac{a^{2}}{12}\partial_{x}^{4}\rho+\cdots\right),\qquad
D\equiv ra^{2},$ (S26)

where the $\cdots$ contain higher derivatives and higher powers of $a$. In the
usual continuum limit we imagine that all of the continuum derivatives of
$\rho_{x}$ are of the same order, and hence after sending $a\rightarrow 0$
only the first term survives.

In order to have a well-defined continuum limit, we require that $\rho_{x}$ be
finite, and not parametrically divergent or vanishing as $a\rightarrow
0,r,L\rightarrow\infty$. In terms of the original lattice model, this means
that we must consider the dilute limit where $n_{i}\rightarrow 0$. In other
words, this continuum limit corresponds to keeping the number of particles
present in the system fixed while simultaneously fine-graining the lattice.

This type of continuum limit — one which corresponds to the dilute limit of a
lattice model with lattice spacing $a\rightarrow 0$ — is not the correct one
to take for the kinematically constrained systems studied in this paper.
Indeed, if we let the number of lattice sites over which any given elementary
hopping process acts be finite, the kinetic constraints then completely
prevent all motion in this limit, and the particles all become trivially
localized. We can get around this by simultaneously taking the range of the
elementary particle hopping processes to diverge as $1/a$, but this produces a
system (a dilute multipole-conserving continuum gas) whose dynamics is
potentially very different from the non-dilute lattice gases studied in the
present work.

Instead, the limit we are interested in is one where the average number of
particles on each lattice site is nonzero, and for typical configurations is
larger than the critical value below which the dynamics begins to freeze out
(as mentioned above, this critical value is set by the range of the elementary
hopping processes [29]). For $\rho_{x}$ to be well-defined we thus need to
keep the lattice spacing fixed, with a derivative expansion of $\rho_{x}$
being performed simply under the assumption that near equilibrium, high-order
derivatives of $\rho_{x}$ are numerically smaller than low-order ones, but not
parametrically smaller in any scaling parameter. For example, in the case of
particles living on the half-line $x>0$ we have (4)

$\rho^{\rm eq}(x)=\frac{N}{x_{\rm cm}}e^{-x/x_{\rm cm}},$ (S27)

and higher derivatives are thus suppressed in powers of $a/x_{\rm cm}$, which
is finite (but still usually rather small) in the limit we are interested in.

We note that these sorts of subtleties, namely the existence of multiple
inequivalent ways of taking the continuum limit, have also recently been
discussed from a field theoretic point of view in e.g. Ref. [43].

###  II.2 Different forms of the master equation

We now take a critical look at the 1d dipolar master equation, and show that
it can take on qualitatively different forms depending on the type of
elementary hopping process considered. This analysis will also help further
motivate the preceeding discussion of continuum limits and will help sharpen
up the meaning of “high density regime”.

As discussed in the main text, the 4-site elementary hopping process
“$0110\leftrightarrow 1001$” captured by (5), (6) gives rise to the master
equation

$\partial_{t}\rho_{x}=r(\rho_{x-1}\rho_{x+2}+\rho_{x-2}\rho_{x+1}+\rho_{x+1}\rho_{x+2}+\rho_{x-2}\rho_{x-1})-r\rho_{x}(\rho_{x-1}+\rho_{x+1}+\rho_{x-3}+\rho_{x+3}),$
(S28)

where $r$ is again the hopping rate. Upon performing a derivative expansion of
$\rho$ as in the main text, we obtain the dipolar diffusion equation
$\partial_{t}\rho=-\partial_{x}^{2}J_{4-{\rm site}}$ with current

$J_{4-{\rm
site}}=4ra^{5}((\partial_{x}\rho)^{2}-\rho\partial_{x}^{2}\rho+a^{2}(\cdots)),$
(S29)

where the $\cdots$ are suppressed by higher powers of derivatives.

It is instructive to consider modifying the form of the elementary dipole-
conserving hopping processes that go into the derivation of the master
equation. The most interesting modifications are those where two particles
leave from or arrive at the same site (other modifications simply lead to
different values for $D$). For example, consider the 3-site
“$020\leftrightarrow 101$” hopping process, which is a more compact version of
the 4-site process considered above. The change in particle occupancies for a
hopping process occuring on sites $x,x+1,x+2$ is

$(\Delta\rho_{x},\Delta\rho_{x+1},\Delta\rho_{x+2})=(\pm 1,\mp 2,\pm 1)$ (S30)

with the top sign corresponding to an “out” hopping process and the bottom to
an “in” process. Within the same type of mean-field approximation used in the
main text, the rate for the “in” process is $rdt\rho_{x}\rho_{x+2}$, while the
rate for the “out” process is proportional to the number of pairs on site $i$,
viz. $rdt\rho_{x+1}(\rho_{x+1}-1)$ (it is easy to verify that these rates
satisfy detailed balance).

The important part here is that the rate for the “out” process contains a term
linear in $n$. This produces a master equation

$\partial_{t}\rho_{x}=r(2\rho_{x-1}\rho_{x+1}+\rho_{x+1}(\rho_{x+1}-1)+\rho_{x-1}(\rho_{x-1}-1)-\rho_{x}(2(\rho_{x}-1)+\rho_{x+2}+\rho_{x-2}))$
(S31)

which contains terms linear in $n$. Performing a derivative expansion then
yields $\partial_{t}\rho=-\partial_{x}^{2}J_{3-{\rm site}}$, with the current
now given by

$J_{3-{\rm
site}}=ra^{5}((\partial_{x}\rho)^{2}-\rho\partial_{x}^{2}\rho+a^{2}(\cdots)-a^{-3}\rho),$
(S32)

where the $\cdots$ are terms quadratic in $\rho$ that contain four
derivatives. This expression for the current differs from the 4-site version
by the extra $-a^{-3}\rho$ term on the RHS, which represents a diffusive term
with a negative effective diffusion constant, favoring charge localization
(interestingly, general hydrodynamic arguments can be given which would have
appeared to rule out the existence of such a term [27]).

Dropping the terms in $(\cdots)$, the zero-current solutions obtained from
(S32) are no longer simple exponentials, and instead take the form

$\rho^{\rm eq}_{3-{\rm site}}(x)=\frac{a^{-3}}{2\alpha^{2}}\sinh^{2}(\alpha
x+\beta),$ (S33)

where $\alpha,\beta,$ are fixed by the total number of particles and total
dipole moment.

Note that if we were to take the standard continuum limit (in which we send
$a\rightarrow 0$ while keeping $\rho$ finite), the last $a^{-3}\rho$ term in
$J_{3-{\rm site}}$ would dominate, and the physics would look very different
to the 4-site variant. Despite this, we emperically find that simulations run
with only 3-site hoppings give results essentially identical to those run with
only 4-site hoppings, at least for the large-density regimes studied in the
main text, where the dynamics is more universal (in both cases the steady-
state distributions are exponentially localized, and the dynamic exponent is
$z=5$). Accordingly, these simulations must be run in a regime where the
$(\partial_{x}\rho)^{2},\rho\partial_{x}^{2}\rho$ terms are not subleading
compared to $a^{-3}\rho$. This necessitates that the average density be larger
than the amount by which higher derivatives of $\rho$ are suppressed. For the
distribution with characteristic length scale $\ell$ (e.g. $\ell=x_{\rm cm}$
for (S27), or $\ell=\sigma$ for the Gaussian
$\rho(x)=\frac{N}{\sqrt{2\pi}\sigma}e^{-x^{2}/2\sigma^{2}}$, the
$(\partial_{x}\rho)^{2},\rho\partial_{x}^{2}\rho$ terms dominate when

$N\gtrsim(\ell/a)^{3},$ (S34)

which is satisfied in the “high-density” situations considered in the main
text. For now we will content ourselves with this characterization, but in the
future it will be valuable to precisely characterize the universal features
(or lack thereof) exhibited by general multipolar dynamics in different
parameter regimes.

##  III minimal nonlinearities of the consitutive relation

In this section we lower-bound the minimal number $m$ of powers of $\rho$ that
must enter into the constitutive relation for the current in $n$-pole
conserving diffusive dynamics, restricting to 1d lattice systems for
simplicity. In the main text we saw that the current was quadratic in $\rho$
($m=2$) in the dipole-conserving case ($n=1$), and cubic in $\rho$ ($m=3$) in
the quadrupole-conserving case ($n=2$) 444The nonlinearity measured by $m$
discussed here refers to the largest power of the density appearing on the RHS
of the master equation. When the hopping process is such that at most one
particle hops from any given site, all terms on the RHS will be of the same
order. In other cases (such as the 3-site dipole hopping process discussed in
App. II) smaller powers may additionally appear. In this appendix we are
interested in determining the minimal value that the largest power on the RHS
can be.. By physicist’s induction, it is then tempting to conclude that
$m=n+1$. While this (likely) is actually not true in complete generality, we
are able to prove the following proposition, which show that $m\geq n+1$:

Proposition.  Consider particles moving on a 1d lattice under translation-
invariant dynamic that conserves the $n$-pole moment of the charge. Then the
constitutive relation expressing the current in terms of the density $\rho$
must involve at least $m\geq n+1$ powers of $\rho$, and at least $n+1$ spatial
derivatives. This lower bound can be saturated for at least all $n\leq 10$, as
well for $n=12$.

Proof: The number $m$ is determined by finding the smallest number of
particles that can participate in an $n$-pole preserving particle hopping
process. For a process involving $m$ particles, let
$\mathbf{a},\mathbf{b}\in\mathbb{Z}^{m}$ denote the initial and final
positions of the particles, respectively. Note that if $a_{i}=b_{j}$ for any
$i,j\in 1,\dots,m$, then at least one particle does not actually move during
the hopping process; thus wolog we may restrict to
$\\{a_{i}\\}\cap\\{b_{i}\\}=\emptyset$.

For translation-invariant dynamics, the conservation of the $n$th multipole
moment implies the conservation of the $p$-pole moments for all $p<m$. Thus
since the process $\mathbf{a}\mapsto\mathbf{b}$ must preserve all multipole
moments up to $n$, we require that

$\sum_{i}a_{i}^{p}=\sum_{i}b_{i}^{p}\,\,\forall\,\,p\in\\{0,\dots,n\\}.$ (S35)

The goal is thus to find the smallest-length pair of disjoint integer-valued
sets such that the above equation is satisfied.

We now show that (S35) can be satisfied only when $m>n$. Indeed, suppose that
$m=n$. It is a well-known fact [45] that the set of “power-sum polynomials”

$P_{p}(\mathbf{a})\equiv\sum_{i=1}^{m}a_{i}^{p}$ (S36)

for $p\in\\{0,\dots,m\\}$ generate all symmetric polynomials in $m$ variables.
This means that the symmetric polynomials
$\prod_{i=1}^{m}a_{i},\prod_{i=1}^{m}b_{i}$ can be expressed as linear
combinations of the $P_{p}(\mathbf{a}),P_{p}(\mathbf{b})$, respectively. But
since we have assumed that all $o\leq m$ multipole moments are conserved, we
have

$P_{p}(\mathbf{a})=P_{p}(\mathbf{b})\,\,\forall\,\,p\in\\{0,\dots,m\\}\,\,\implies\,\,\prod_{i=1}^{m}a_{i}=\prod_{i=1}^{m}b_{i}.$
(S37)

Now we use the fact that (S35) is invariant under translations, in that (S35)
implies

$P_{p}(\mathbf{a}-x)=P_{p}(\mathbf{b}-x)$ (S38)

for any constant vector $x$ (this follows from the fact that
$P_{p}(\mathbf{a}-x)$ is also a symmetric polynomial). In particular, we may
choose $x=a_{1}$ to set $a_{1}=0$ wolog. But then from the above we see that
since $\prod_{i}a_{i}=0$, we must also have $\prod_{i}b_{i}=0$. In particular,
one of the $b_{i}$ must equal zero, which contradicts our assumption about the
disjointness of $\\{a_{i}\\},\\{b_{i}\\}$.

This shows that any $n$-pole conserving hopping process must involve at least
$n+1$ particles, so that the constitutive relation must involve at least $n+1$
powers of the density. This also implies that the current must involve at
least $n+1$ derivatives; if it involved fewer than $n+1$ derivatives then it
would be proportional to an overall power of the density $J=\rho^{k}(\cdots)$
(where $\cdots$ involves derivatives of $\rho$), and the $\rho^{k}$ in front
could be stripped away without affecting the nature of the zero-current
solutions.

We have shown that $n+1$-body terms are necessary, but not that they are
sufficient. Showing this would amount to proving that one can always find two
mutually disjoint size-$n+1$ sets of integers satisfying (S35), which is in
fact a famous open problem in number theory known as the Prouhet-Tarry-Escott
problem [34]. We will of course not attempt to make progress on this problem
here, and simply note that solutions are known to exist for $n\leq 10$ and
$n=12$, but that for other values of $n$ the answer is unknown (with the
existence of solutions at large $n$ seeming unlikely) [34]. ∎

The explicit $m=n+1$ solutions of (S35) (when they exist) quickly get rather
complicated for large $n$, meaning that they involve very long-range hopping
processes. For example, when $n=7$ the simplest solution is an eight-body
hopping process with

$\displaystyle a$ $\displaystyle=[1,5,10,24,28,42,47,51]$ (S39) $\displaystyle
b$ $\displaystyle=[2,3,12,21,31,40,49,50]$

meaning that the minimal-body hopping process extends over 50 lattice sites.

Depending on the physical context it may be unreasonable to consider such
long-range hopping terms, instead considering shorter-range but higher-body
hopping processes. The minimal-range hopping processes involve a number of
particles that scales exponentially in $n$ but has range equal to $n+2$, which
can be constructed with the aid of Pascal’s triangle by letting

$\mathbf{a}=\bigoplus_{i=0}^{\lfloor(n+1)/2\rfloor}(2i)^{\oplus{n+1\choose
2i}},\qquad\mathbf{b}=\bigoplus_{i=1}^{\lfloor(n+2)/2\rfloor}(2i-1)^{\oplus{n+1\choose
2i-1}},$ (S40)

where $\oplus$ denotes concatenation, so that e.g. $x^{\oplus k}$ is the
length-$k$ vector with all entries equal to $x$. For dipoles this gives the
2-body process $\mathbf{a}=(1,3),\,\mathbf{b}=(2,2)$, for quadrupoles the
4-body process $\mathbf{a}=(1,3,3,3),\,\mathbf{b}=(2,2,2,4)$, and so on.
Processes involving a non-minimal number of particles can be checked to lead
to difussion equations with additional powers of $\rho$ appearing in the
constitutive relation for the current. For example, the (non-minimal) 4-body
quadrupole-conserving process just described gives rise to the relation
$J=D\rho^{4}\partial_{x}^{3}\ln(\rho)$, which has one more power of $\rho$ as
compared to the minimal solution $J=D\rho^{3}\partial_{x}^{3}\ln(\rho)$.

##  IV Derivation of the dipolar Langevin equation

The dipolar diffusion equation given in the main text is deterministic, and
misses the stochastic nature of the pairwise hopping. Properly including
fluctuation effects can be done by statistically analyzing the full
microscopic motion and coarse-graining the result, as laid out e.g. in Ref.
[32]. We will use notation appropriate for a 1d chain; generalizations to
higher dimensions are straightforward. We will work in units where the lattice
spacing $a=1$ throughout.

In each time step $t\to t+dt$, each pair of particles in the system undergoes
pairwise diffusion with a rate $r$, causing the number of particles in each
lattice to change by $\Delta\rho_{x}(t)\equiv\rho_{x}(t+dt)-\rho_{x}(t)$. Then
$\Delta\rho_{x}(t)$ for all possible $x$ and $t$ together determines the
trajectory of the system.

Next, we introduce the moment-generating function of the particle dynamics

$W[\hat{\rho}_{x}(t)]=\Big{\langle}\exp\Big{[}\sum_{x,t}\hat{\rho}_{x}(t)\Delta\rho_{x}(t)\Big{]}\Big{\rangle},$
(S41)

where $\hat{\rho}_{x}(t)$ is the conjugate field of $\Delta\rho_{x}(t)$ and
the average is taken over all possible trajectories of $\Delta\rho_{x}(t)$. In
the continuum time limit, functional derivatives of $\ln W$ with respect to
$\hat{\rho}$ can be used to calculate connected correlation functions of
$\partial_{t}\rho$.

Since the pairwise diffusive jump processes are independent, we can rewrite
the generating function as

$W[\hat{\rho}_{x}(t)]=\prod_{x,t}\Big{\langle}\exp\Big{[}\hat{\rho}_{x}(t)\Delta\rho_{x}(t)\Big{]}\Big{\rangle}.$
(S42)

For a lattice with $L$ sites and periodic boundary conditions, there are $2L$
possible jump processes; the ingoing ($1001\to 0110$) and outgoing ($0110\to
1001$) process for any given adjacent quartet of sites $(x,x+1,x+2,x+3)$.

To calculate Eq. (S42) over these $2L$ processes, we refer to
$(x,x+1,x+2,x+3)$ sites with the index $i$ and take an average over ingoing
and outgoing processes. To make the subsequent expressions concise, we omit
$t$ from the arguments of $\rho_{x}(t)$ and $\hat{\rho}_{x}(t)$ in the
following. We obtain

$\displaystyle W[\hat{\rho}_{x}(t)]$
$\displaystyle=\prod_{x,t}\Big{[}\rho_{x+1}\rho_{x+2}rdte^{\hat{\rho}_{x}-\hat{\rho}_{x+1}-\hat{\rho}_{x+2}+\hat{\rho}_{x+3}}+\rho_{x}\rho_{x+3}rdte^{-\hat{\rho}_{x}+\hat{\rho}_{x+1}+\hat{\rho}_{x+2}-\hat{\rho}_{x+3}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+1-[\rho_{x+1}\rho_{x+2}+\rho_{x}\rho_{x+3}]rdt\Big{]}.$
(S43)

In the square brackets of the expression above, the first term is the product
between the probability for the outgoing motion to occur,
$\rho_{x}\rho_{x+3}rdt$, and the conjugate field factor associated with the
changes of the particle numbers $\Delta\rho_{x}=\Delta\rho_{x+3}=+1$ and
$\Delta\rho_{x+1}=\Delta\rho_{x+2}=-1$. Similarly, the second term consists of
the probability of the ingoing motion and the conjugate field factor with
$\Delta\rho_{x}=\Delta\rho_{x+3}=-1$ and
$\Delta\rho_{x+1}=\Delta\rho_{x+2}=+1$. The last term is the probability for
no motion to occur, with the corresponding conjugate factor with
$\Delta\rho_{x}=\cdots=\Delta\rho_{x+3}=0$. Now, we introduce

$\displaystyle R_{i}^{+}$ $\displaystyle\equiv r\rho_{x+1}\rho_{x+2}\,,$ (S44)
$\displaystyle R_{i}^{-}$ $\displaystyle\equiv r\rho_{x}\rho_{x+3}\,,$
$\displaystyle\Sigma_{i}$
$\displaystyle\equiv\hat{\rho}_{x}-\hat{\rho}_{x+1}-\hat{\rho}_{x+2}+\hat{\rho}_{x+3}\,$

and write Eq. (S43) as

$\displaystyle W[\hat{\rho}_{x}(t)]$
$\displaystyle=\prod_{i,t}\Big{[}1+dt\big{\\{}R_{i}^{+}(e^{\Sigma_{i}}-1)+R_{i}^{-}(e^{-\Sigma_{i}}-1)\big{\\}}\Big{]},$
(S45)
$\displaystyle\simeq\exp\Big{[}\sum_{i,s}dt\\{R_{i}^{+}(e^{\Sigma_{i}}-1)+R_{i}^{-}(e^{-\Sigma_{i}}-1)\\}\Big{]}.$

To proceed, we take the continuum time limit and expand $\hat{\rho}$ in
gradients (which is allowed as long as we are only interested in long-distance
correlation functions of $\Delta\rho$), allowing us to write
$\Sigma_{x}\approx 2\partial_{x}^{2}\rho$. This gives

$\displaystyle\ln W[\hat{\rho}_{x}(t)]$
$\displaystyle=\sum_{x}\int_{t}(R^{+}_{x}(e^{2\partial_{x}^{2}\hat{\rho}}-1)+R^{-}_{x}(e^{-2\partial_{x}^{2}\hat{\rho}}-1)).$
(S46)

Taking the variation with respect to $\hat{\rho}$ and then sending
$\hat{\rho}\rightarrow 0$ gives the average $\langle\partial_{t}\rho\rangle$,
yielding

$\displaystyle\langle\partial_{t}\rho\rangle$
$\displaystyle=2\partial_{x}^{2}(R^{+}_{x}-R^{-}_{x})$ (S47)
$\displaystyle\rightarrow
D\partial_{x}^{2}\left[(\partial_{x}\rho)^{2}-\rho\partial_{x}^{2}\rho\right],$

where in the second line we have defined $D\equiv 4r$ and expanded $\rho_{x}$
up to quadratic order in derivatives. We thus recover the dipolar diffusion
equation, but now with $\partial_{t}\rho$ replaced by the statistical average
$\langle\partial_{t}\rho\rangle$.

Taking two functional derivatives of (S46) gives

$\displaystyle\langle\partial_{t}\rho(x,t)\partial_{t^{\prime}}\rho(x^{\prime},t^{\prime})\rangle_{c}$
$\displaystyle=4\partial_{x}^{2}\partial_{x^{\prime}}^{2}[\delta(x-x^{\prime})\delta(t-t^{\prime})(R^{+}_{x}+R^{-}_{x})]$
(S48) $\displaystyle\rightarrow
2D\partial_{x}^{2}\partial_{x^{\prime}}^{2}\left[\rho^{2}(x,t)\delta(x-x^{\prime})\delta(t-t^{\prime})\right].$

In order to reproduce this result, one must supplement the dipolar diffusion
equation with the fluctuation term, which results in the dipolar Langevin
equation (10) quoted in the main text,

$\partial_{t}\rho(x,t)=-D\partial_{x}^{2}(\rho^{2}\partial_{x}^{2}\ln(\rho))+\sqrt{2D}\partial_{x}^{2}\eta(x,t),$
(S49)

where the noise field $\eta$ obeys
$\langle\eta(x,t)\eta(x^{\prime},t^{\prime})\rangle=\delta(t-t^{\prime})\delta(x-x^{\prime})\rho^{2}(x,t)$.

One application of this is a derivation of the structure factor governing the
behavior of small density fluctuations about equilibirum. To this end, write
$\rho=\rho^{\rm eq}+\delta\rho$, and linearize the Langevin equation in
$\delta\rho$ while dropping all gradients of $\rho^{\rm eq}$ (this amounts to
considering fluctuations on length scales short enough so that $\rho^{\rm eq}$
can be treated as being approximately constant). Then

$\partial_{t}\delta\rho=-\widetilde{D}\partial_{x}^{4}\delta\rho+\sqrt{2D}\partial_{x}^{2}\eta,$
(S50)

where $\widetilde{D}\equiv D(\rho^{\rm eq})^{2}$ is an effective subdiffusion
constant, and
$\langle\eta(x,t)\eta(x^{\prime},t^{\prime})\rangle_{c}=(\rho^{\rm
eq})^{2}\delta(t-t^{\prime})\delta(x-x^{\prime})$. Going to Fourier space
gives

$\langle\delta\rho(k,t)\delta\rho(-k,t)\rangle_{c}=\int\frac{d\omega}{{2\pi}}\frac{2\widetilde{D}\rho^{\rm
eq}k^{4}}{\omega^{2}+(\widetilde{D}k^{4})^{2}}=\rho^{\rm eq},$ (S51)

which is the exact same result one would obtain by applying this formalism to
the regular diffusion equation.