# Einstein relation for subdiffusive relaxation in Stark chains

P. Prelovšek  Jožef Stefan Institute, SI-1000 Ljubljana, Slovenia  S. Nandy
Jožef Stefan Institute, SI-1000 Ljubljana, Slovenia  M. Mierzejewski
Institute of Theoretical Physics, Faculty of Fundamental Problems of
Technology, Wrocław University of Science and Technology, 50-370 Wrocław,
Poland

###### Abstract

We investigate chains of interacting spinless fermions subject to a finite
external field $F$ (also called Stark chains) and focus on the regime where
the charge thermalization follows the subdiffusive hydrodynamics. First, we
study reduced models conserving the dipole moment and derive an explicit
Einstein relation which links the subdiffusive transport coefficient with the
correlations of the dipolar current. This relation explains why the decay
rate, $\Gamma_{q}$, of the density modulation with wave-vector $q$ shows
$q^{4}$-dependence. In the case of the Stark model, a similar Einstein
relation is also derived and tested using various numerical methods. They
confirm an exponential reduction of the transport coefficient with increasing
$F$. On the other hand, our study of the Stark model indicates that upon
increasing $q$ there is a crossover from subdiffusive behavior,
$\Gamma_{q}\propto q^{4}$, to the normal diffusive relaxation,
$\Gamma_{q}\propto q^{2}$, at the wave vector $q^{*}$ which vanishes for $F\to
0$.

Introduction. Macroscopic systems driven by finite external forces/fields are
usually described within the extended concept of local equilibrium via local
thermodynamic parameters as, e.g., the temperature and chemical potential.
Such a description can fail when we are dealing with systems isolated from the
environment. Prominent recent examples are those of the Stark systems of
interacting fermions (which is fermion systems in the presence of a finite
field), as realized in tilted cold-atom lattices [1, 2, 3]. Such systems
reveal several novel theoretical insights and challenges, discussed mostly
within one-dimensional (1D) models. Since noninteracting particles subject to
a finite external field, $F$, exhibit Stark localization, the problem shares
some similarities with the many-body localization (MBL) (for an overview see
[4, 5, 6, 7]) and it is known as the Stark MBL [8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18]. It is also well established that effective models, derived at
large $F$ from the Stark model (usually by invoking the Schrieffer–Wolff
transformations), can exhibit Hilbert space fragmentation [19, 20, 21, 22, 23,
24] that violates the eigenstate thermalization hypothesis (ETH) [25, 26, 27].
In systems with strongly fragmented Hilbert space, the latter can be linked
with the emergence of additional integrals of motion [29, 28].

Still at moderate $F$, the cold-atom experiments [1] as well as numerical
simulations of 1D models [30] reveal the hydrodynamic relaxation of the
inhomogeneous particle distributions towards the steady state, which
corresponds to the infinite-temperature ($T\to\infty$) equilibrium. For small
wave-vectors $q\ll 1$ the relaxation rates follow a particular subdiffusive
(SD) law $\Gamma_{q}\propto q^{4}$, rather than the normal diffusive behavior,
$\Gamma_{q}\propto q^{2}$. This is well consistent with the fact that at
$F\neq 0$, the dipole moment emerges in macroscopic systems as an additional
conserved quantity [19, 31, 2, 3, 30] and the phenomenological description can
be given in terms of the fracton hydrodynamics [32, 31, 33, 34, 35, 36]. In
spite of a wide consensus that such a description is appropriate for isolated
driven systems, so far the theoretical studies are mostly based on
phenomenological hydrodynamic approaches. There are so far very few
quantitative results for the subdiffusion in microscopic models [30] as well
as theoretical attempts to express the SD transport coefficient $D_{S}$ [1,
37] in terms of the response functions.

In this Letter, we present the analysis, focused on the particle/density
relaxation and anomalous diffusion in the hydrodynamic regime. We consider the
standard 1D Stark model, i.e., the chain of interacting spinless fermions
subject to a finite external field $F$, as well as the effective models which
involve extended pair-hopping (EPH) interactions and strictly conserve the
dipole moment, $P$. On the one hand, from numerical results for the dynamical
density structure factor $S_{q}(\omega)$, we extract the relaxation rates
$\Gamma_{q}$ of the density modulation. These rates reveal subdiffusive
transport for $q\to 0$ as well as the value of the corresponding SD
coefficient $D_{S}$. On the other hand, employing the memory-function (MF)
formalism, we derive the Einstein relation that expresses $\Gamma_{q\to 0}$ in
terms of the uniform ($q=0$) correlations of the normal current, $J_{N}$, and
the dipolar current, $J_{D}$. If $P$ is conserved the response function is
determined solely by $J_{D}$ and the relaxation follows the SD relation
$\Gamma_{q}=D_{S}q^{4}$. In the EPH model, $J_{D}$ is a translationally-
invariant operator that governs the relaxation. In the case of the full Stark
model, $J_{D}$ dominates the response for large $L$ (for $F>0$) when one
observes the emergent conservation of the dipole moment. The Einstein relation
as well as the numerical results for $D_{S}$ in both models are tested with
alternative numerical approaches. It should be emphasized that the derived
Einstein relations remain valid beyond the considered models and even beyond
1D, which we mostly discuss below. Moreover, our numerical results in Stark
chains for larger $q>0$ reveal the crossover from SD $\Gamma_{q}\sim
D_{S}q^{4}$ to normal diffusion $\Gamma_{q}\sim D_{N}q^{2}$ at $q\sim
q^{*}(F)$ with vanishing $q^{*}(F\to 0)$, consistent with some
phenomenological theories [1, 37].

Dynamical density-modulation relaxation. In the following, we study two 1D
lattice models of interacting spinless fermions, as they emerge in the
presence of the finite external field $F$, whereby the chain has $L$ sites and
open boundary conditions (OBC). The coupling to the field enters the
Hamiltonian via $H^{\prime}=FP$ where $P$ is the dipole moment,
$P=\sum_{l}(l-L/2)n_{l}$, and $n_{l}$ is the particle number operator at site
$l$.

Isolated macroscopic Stark systems at finite $F>0$ develop (heat up) towards a
homogeneous steady state $\langle n_{l}\rangle=\bar{n}=N/L$ ($N$ representing
the total particle number), corresponding to $T\to\infty$ equilibrium. Further
on, we analyze dynamics of the periodic density modulation
$n_{q}=\sum_{l}e^{iql}\tilde{n}_{l}/\sqrt{L}$, $\tilde{n}_{l}=n_{l}-\bar{n}$,
and related correlation function $\phi_{q}(\omega)$,

$\displaystyle\phi_{q}(\omega)$ $\displaystyle=$
$\displaystyle\frac{\chi_{q}(\omega)-\chi^{0}_{q}}{\omega}=\frac{-\chi^{0}_{q}}{\omega+M_{q}(\omega)},$
$\displaystyle\chi_{q}(\omega)$ $\displaystyle=$
$\displaystyle\frac{i}{\beta}\int_{0}^{\infty}dte^{i\omega
t}\langle[n_{-q}(t),n_{q}]\rangle,$ (1)

whereby we define dynamical susceptibilities $\chi_{q}(\omega)$ and
$\chi^{0}_{q}=\chi_{q}(\omega=0)$ that remain nonzero even at $\beta=1/T\to
0$. In this limit, $\phi_{q}(\omega)$ is related with the standard dynamical
structure factor $S_{q}(\omega)=\mathrm{Im}\phi_{q}(\omega)/\pi$. In general,
$\phi_{q}(\omega)$ can be represented in terms of the memory function (MF),
$M_{q}(\omega)$, that determines the profile relaxation rate
$\Gamma_{q}=\mathrm{Im}M_{q}(\omega=0)$, which we later on extract also from
numerical results for $\phi_{q}(\omega)$.

Einstein relation. An analytical step towards $M_{q}(\omega)$ can be made
using the MF formalism [38, 39, 40, 41], with the introduction of the general
scalar product of two operators, $(A|B)$, and the Liouville operator ${\cal
L}A=[H,A]$. In the case of $\beta\to 0$, the latter scalar product reduces to
thermodynamic average, i.e., $(A|B)\sim\langle A^{\dagger}B\rangle$. Within
this formalism one can express the correlation function as
$\phi_{q}(\omega)=(n_{q}|({\cal L}-\omega)^{-1}|n_{q})$. The memory function
can be written in the hydrodynamic regime $q\to 0$ [39] (in analogy to the
perturbation theory [40]) as

$M_{q}(\omega)=({\cal L}n_{q}|({\cal L}-\omega)^{-1}|{\cal
L}n_{q})/\chi^{0}_{q}.$ (2)

Expanding $n_{q}$ in powers of $q$ one obtains,

${\cal L}n_{q}\simeq\frac{1}{\sqrt{L}}\left(iq{\cal L}P+iq^{2}\frac{i}{2}{\cal
L}Q\right),\quad Q=\sum_{l}l^{2}\tilde{n}_{l},$ (3)

where we assumed conservation of the particle number, ${\cal L}N=0$ with
$N=\sum_{l}n_{l}$. The first term in Eq.(3) represents the normal (uniform)
current, $J_{N}=i{\cal L}P$. It determines the hydrodynamic relaxation ($q\to
0$) in generic systems that do not conserve the dipole moment, ${\cal L}P\neq
0$. Namely, one obtains from Eq.(2) the standard Einstein relation [42, 39,
43, 44], $\Gamma_{q}=q^{2}\mathrm{Im}\phi_{N}(0)/\chi^{0}_{0}=D_{N}q^{2}$,
which links the diffusion constant, $D_{N}$, with the current correlation
function $\phi_{N}(\omega)=(J_{N}|({\cal L}-\omega)^{-1}|J_{N})/L$.

If the dipole moment is conserved, ${\cal L}P=0$, then the hydrodynamic
relaxation is determined by the second term in the expansion in Eq. (3), which
can be interpreted as the dipolar current $J_{D}=\frac{i}{2}{\cal L}Q$.
Similarly, one then obtains then the Einstein relation $\Gamma_{q}=D_{S}q^{4}$
with $D_{S}=\mathrm{Im}\phi_{D}(0)/\chi^{0}_{0}$, however, involving the
dipolar currents $\phi_{D}(\omega)=(J_{D}|({\cal L}-\omega)^{-1}|J_{D})/L$. To
conclude this part, we note that the MF formalism straightforwardly explains
the origin of SD in systems with dipole-moment conservation,
$\Gamma_{q}=D_{S}q^{4}$. Similarly to the phenomenological fracton
hydrodynamics [32, 31, 33, 34, 35, 36], the only assumption we made is the
conservation of the dipole moment.

Extended pair-hopping (EPH) model. As a first example we analyze the model
where $P$ is strictly a conserved. Starting from the full Stark models, one
can derive at $F\gtrsim 1$ the EPH model, either via the Schrieffer-Wolff
transformation [19, 20, 31], or expanding the interaction in the Stark basis
[28],

$H_{EPH}=\sum_{l}\zeta_{dr}[c_{l-r}^{\dagger}c_{l}c_{l+d+r}^{\dagger}c_{l+d}+\mathrm{H.c.}]+H_{d}.$
(4)

Here, $c^{\dagger}_{l},c_{l}$ refer to localized Stark states and $H_{d}$
represents the Hartree-Fock diagonal term, while $\zeta_{dr}$ can be derived
for given $F$, at least to lowest order in the interaction [28]. Here, by
construction ${\cal L}P=0$. One can derive explicit expression for the dipolar
current $J_{D}=\frac{i}{2}{\cal L}Q=\sum_{dr}\zeta_{dr}J_{D}(d,r)$, where

$J_{D}(d,r)=-r(r+d)\sum_{j}(ic_{j-r}^{\dagger}c_{j}c_{j+d+r}^{\dagger}c_{j+d}+\mathrm{H.c.}).$
(5)

It is remarkable that explicit $l$-dependence cancels out and $J_{D}$ emerges
as a translationally-invariant operator.

Figure 1: Dynamical structure factor $q^{4}S_{q}(\omega)$ for the EPH model
with $L=32$ sites, obtained via MCLM for the lowest $q=2m\pi/L,m=1,2$. The
inset shows the extracted
$\tilde{D}_{S}(\omega)=\mathrm{Im}M_{q_{1}}(\omega)/q_{1}^{4}$, compared to
the Einstein-relation result
$D_{S}(\omega)=\mathrm{Im}\phi_{D}(\omega)/\chi^{0}_{0}$.

Here, we do not aim to investigate closer the EPH models with realistic
parameters, but rather test the existence of SD and the direct expression for
the coefficient $D_{S}$ for a simplified case with $H_{d}=0$ and (rather
arbitrary) assuming $r=1$ and values $\zeta_{d1}=1,0.75,0.5,0.25$ for
$d=1,...,4$. The motivation for including longer-range $d>1$ terms is that the
basic pair-hopping model with only $\zeta_{11}$ is known to exhibit strong
Hilbert space fragmentation [19, 20, 22, 23] which invalidates basic ETH
concept, while additional terms with $\zeta_{d>1,r>1}$ are expected to supress
this effect [28]. In the following we calculate numerically
$\phi_{q}(\omega)$, Eq. (1), using the microcanonical Lanczos method (MCLM)
for finite $L$ systems [45, 46, 47], employing large number of Lanczos steps
up to $N_{L}\sim 10^{5}$, to achieve the frequency resolution
$\delta\omega\lesssim 10^{-4}$ which allows for reliable extraction of
$M_{q}(\omega)$ [47] even for small $q\ll 1$. The advantage of the EPH model
(relative to the Stark model) is that in addition to $N$ we use also
conservation of $P$ (choosing $P=0$) to reduce the Hilbert space and to reach
$L=32$ with $N_{st}\sim 10^{7}$ basis states.

In Fig. 1 we present results for the density structure factor
$q^{4}S_{q}(\omega)$ as calculated via MCLM for two lowest
$q=q_{m}=2m\pi/L,m=1,2$. Since we consider the half-filled system
$\bar{n}=1/2$ with effective $T\to\infty$, we know analytically
$\chi^{0}_{q}\sim\chi_{0}^{0}\sim\bar{n}(1-\bar{n})=1/4$. Results confirm very
sharp peak at $\omega\sim 0$, being consistent with SD hydrodynamics, i.e.,
$\pi S_{q}(\omega\sim
0)\sim\chi_{q}^{0}/\Gamma_{q}\sim\chi_{0}^{0}/(D_{S}q^{4})$. Moreover, we
extract also the corresponding (dynamical) SD coefficient
$\tilde{D}_{S}(\omega)=\mathrm{Im}M_{q}(\omega)/q^{4}$ for smallest $q=q_{1}$
and present results in the inset, together with numerically evaluated Einstein
relation $D_{S}(\omega)=\mathrm{Im}\phi_{D}(\omega)/\chi^{0}_{0}$, using
$J_{D}$ from Eq. (5). The agreement is reasonable given that both numerical
approaches can suffer from finite-size effects. Moreover, the considered EPH
model can still exhibit some features of the Hilbert-space fragmentation [19,
20, 22, 23], which could influence the presumed ETH.

Stark model. We turn further to the properties of the prototype Stark model,
i.e., 1D chain of interacting spinless fermions in the presence of a finite
external field $F$,

$\displaystyle H$ $\displaystyle=$ $\displaystyle
t\sum_{i}(c^{\dagger}_{l+1}c_{l}+c^{\dagger}_{l}c_{l+1})+V\sum_{l}\tilde{n}_{l+1}\tilde{n}_{l}+$
(6) $\displaystyle+$ $\displaystyle
V^{\prime}\sum_{l}\tilde{n}_{l+2}\tilde{n}_{l}+FP,$

with $\tilde{n}_{l}=n_{l}-\bar{n},n_{l}=c^{\dagger}_{l}c_{l}$. Fermions
interact via the nearest-neighbor ($V$) and next-nearest-neighbor
($V^{\prime}$) repulsion. We consider half-filling, i.e., $\bar{n}=N/L=1/2$,
and set $t=1$ as the unit of energy. We introduce $V^{\prime}\neq 0$ in order
to suppress the integrability (and dissipationless transport) at $F\to 0$,
although the latter effect appears not to be important for $F\gg 0$. In the
main text, we restrict the numerical results to the case $V=V^{\prime}=1$,
while in [48] we discuss also results for $V=2,V^{\prime}=0$.

In the case of the Stark model, Eq. (6), we cannot apply the same analysis as
for EPH model, since ${\cal L}P\neq 0$ and the conservation of $P$ emerges
only in the thermodynamic limit, $L\to\infty$ [30]. Still, one can derive from
Hamiltonian. (6) both contributions to ${\cal L}n_{q}$ in Eq. (3):
$J_{N}=i{\cal L}P=\sum_{l}J_{l}$ and $J_{D}=\frac{i}{2}{\cal
L}Q=J_{N}/2+\sum_{l}lJ_{l}$, where $J_{l}=itc^{\dagger}_{l}c_{l+1}+{\rm
H.c.}$. Neglecting possible offdiagonal correlations, we can then express the
corresponding Einstein relation from Eq. (2),

$M_{q}(\omega)\simeq[q^{2}\phi_{N}(\omega)+q^{4}\phi_{D}(\omega)]/\chi^{0}_{0}.$
(7)

Since in general $\phi_{N}(\omega)\neq 0$, we can proceed by showing that for
finite $F>0$ and $L\to\infty$ it follows $\phi_{N}(\omega\to 0)=0$, related to
emergent conservation of $P$ [30], as verified also numerically in [48].
Consequently, we expect the at $F>0$ the hydrodynamic $q\to 0$ behavior will
be dominated by SD with $\Gamma_{q}=D_{S}q^{4}$ with $D_{S}={\rm
Im}\phi_{D}(\omega\to 0)/\chi_{0}^{0}$.

Figure 2: (a) Dipolar-current correlations $\mathrm{Im}\phi_{D}(\omega)$ (in
log scale) as calculated with MLCM for the Stark model on $L=28$ sites for
different fields $F$. (b) Subdiffusion coefficient $D_{S}(F)$, calculated from
the rates $\Gamma_{q_{1}}$, obtained via: time-evolution (TE) of the density
profile and via MCLM on $L=24$ chain, and directly from the Einstein relation
(ER), i.e., from $\phi_{D}(\omega)$.

We further present numerical results for $\mathrm{Im}\phi_{D}(\omega)$ in Fig.
2(a). Since by construction $J_{D}$ requires OBC, one cannot guarantee
$L$-independent result for $\phi_{D}(\omega\to 0)$. Still, results in Fig.
2(a), obtained with MCLM on $L=28$ chain, indicate that beyond
$F>F_{*}(L)\gtrsim 0.4$, there is a well defined value
$D_{S}=\mathrm{Im}\phi_{D}(\omega\to 0)/\chi^{0}_{0}$, revealing already an
exponential-like dependence on $F$. We also note that
$\mathrm{Im}\phi_{D}(\omega)$ for larger $\omega$, as in Fig. 2(a), does not
have a direct relation to $M_{q}(\omega)$, since it neglects
$\phi_{N}(\omega)$ contribution in Eq. (7).

Results for the transport coefficient $D_{S}(F)$ are summarized in Fig. 2(b).
Besides the results from the Einstein relation evaluated via
$\phi_{D}(\omega\to 0)$ (see Fig. 2(a)), we include also results for
$D_{S}=\Gamma_{q}/q^{4}$ obtained from two alternative approaches applied for
$L=24$ chain and the smallest $q_{1}=2\pi/L$. Namely, we extract $\Gamma_{q}$
directly from the MCLM results for $S(q,\omega\to 0)$, as well as from the
decay of the inhomogeneous density profile where $F$ is introduced via the
time-dependent flux [49, 30]. The latter method evaluates the relaxation rate
$\Gamma_{q}$ (see [50] for the details), with the advantage of periodic
boundary conditions and consequently resolving also very small $D_{S}$, i.e.,
reaching larger $F\simeq 2$. Results in Fig. 2(b) are quantitatively
consistent in the broad regime of $F>F_{*}(L)\sim 0.4$, confirming the
validity of the Einstein relation for SD transport as well as the exponential
dependence of $D_{S}(F)$. In [48] we present results also for other
parameters, in particular for $V=2,V^{\prime}=0$, where we employ also the
boundary-driven open systems [30, 50], allowing for the analysis of
considerably larger $L\leq 50$ as well as the direct evaluation of $D_{S}$ via
the form of the nonequilibrium stationary density profile [48].

Crossover to normal diffusion. Let us finally examine closer $M_{q}(\omega)$
at larger $q$ as extracted numerically from $\phi_{q}(\omega)$, again for
parameters $V=V^{\prime}=1$. Fig. 3 shows MCLM results for ${\rm
Im}M_{q}(\omega)/q^{2}$ obtained for various $q=2m\pi/L$ with $1\leq m\leq 4$
and $L=24$ so that $q\leq\pi/3$. One can resolve three regimes. Despite
finite-size limitations, results at small $F=0.2<F_{*}(L)$, shown in Fig.
3(a), are approximately consistent with normal diffusion for all presented
$q_{m}$, i.e., ${\rm Im}M_{q}(\omega)\propto q^{2}$. At intermediate
$F_{*}(L)<F<F_{c}\sim 1$, as in Fig. 3(b), only the smallest $q=q_{1}$
evidently deviates, the latter being the signature of the SD transport ${\rm
Im}M_{q}(\omega\sim 0)\propto q^{4}$. Still, the relaxation functions ${\rm
Im}M_{q}(\omega)/q^{2}$ nearly overlap for larger $q$ which can be interpreted
as an effective normal diffusion ${\rm Im}M_{q}(\omega)\propto q^{2}$.
Finally, for large $F\gtrsim F_{c}$, as in Fig. 3(c), an anomalous SD-like
relaxation appears for all $q<1$.

Figure 3:  Dynamical relaxation rates $\Gamma_{q}(\omega)/q^{2}$ (in log
scale), as extracted from MCLM results for $\phi_{q}(\omega)$ on $L=24$ chain
and $q_{m}=2m\pi/L,m=1-4$, for three different regimes of $F$.

In Fig. 4(a) we collect results for an effective normal diffusion coefficient
$\tilde{D}_{N}={\rm Im}M_{q}(\omega=0)/q^{2}$ at largest $q=q_{4}=\pi/3$ as in
Fig 3. Taking results for $D_{S}$ from Fig. 2(b) and $\tilde{D}_{N}$ from Fig.
4(a) one may estimate that the crossover between SD and diffusive relaxations
is expected at the wavevector $q^{*}$ fulfilling the relation
$D_{S}(q^{*})^{4}=\tilde{D}_{N}(q^{*})^{2}$. Numerical results from such
estimate are shown in Fig. 4(b), representing a rough phase diagram of normal-
SD transport, relevant at least for moderate $F<F_{c}$. While we expect a
continuous vanishing of $q^{*}$ for $F\to 0$, it is hard to numerically
determine the dependence $q^{*}(F\ll 1)$ due to finite-size limitations.
Still, the general trend of $q^{*}(F)$ is well consistent with
phenomenological approaches [1, 37]. We should also note that the explicit MF
expression, Eq. (7), is restricted only to the hydrodynamic regime $q\to 0$
when $M_{q}(\omega)$ is small and $\phi_{N}(\omega\to 0)$ strictly vanishes.
Therefore, it cannot be extended to the discussion of the normal/SD crossover
at larger $q$.

Figure 4:  (a) Effective normal diffusion coefficient
$\tilde{D}_{N}=\Gamma_{q}/q^{2}$ for $q=\pi/3$ vs. $F$, evaluated at $q=\pi/3$
via MCLM on $L=24$ sites. (b) Effective phase diagram $q(F)$ with the
SD/normal diffusion crossover at $q^{*}(F)$, with the dotted line representing
the qualitative guess for small $F$.

Conclusions. In this Letter, we present the analysis of the Stark models,
based on the memory-function approach to hydrodynamics. It enables a direct
consideration of the subdiffusion and the corresponding transport coefficient
$D_{S}$. In models with strictly conserved dipole moment $P$, the derivation
yields a subdiffusive relaxation rates, $\Gamma_{q}=D_{S}q^{4}$ and an
explicit Einstein relation which links $D_{S}$ with correlations of the
uniform ($q=0$) dipolar current $J_{D}$. In the full Stark problem, where $P$
is conserved only in the thermodynamic limit, the analogous treatment is valid
only in more restricted hydrodynamic regime $q\to 0$. However, the
correlations of the dipolar current are still linked with $D_{S}$ via the
Einstein relation. Moreover, the obtained numerical results agree with
alternative numerical approaches, at least in the range of $F$ where the
finite-size limitations allow for reliable numerical studies. As a general
observation, we note that the coefficient $D_{S}$ reveals a strong
exponential-like reduction with $F$, which for large $F\gg 1$ can effectively
appear as Stark MBL concept [8, 9, 10, 11, 12, 13, 15, 14, 16, 17, 18], as
well as the Hilbert-space fragmentation [19, 20, 21, 22, 23, 24]. It should be
stressed that the presented Einstein relations are not specific to considered
model and even not to one-dimensional systems. They can be easily generalized
to models which are more relevant for experiments as, e.g., the tilted Fermi-
Hubbard model.

###### Acknowledgements.

We acknowledge the support by Slovenian Research and Innovation Agency (ARIS)
via the program P1-0044 (P.P., S.N.) and the project J1-2463 (S.N.) as well as
the National Science Centre, Poland via project 2020/37/B/ST3/00020 (M.M.). S.
N. acknowledges also the EU support via QuantERA grant T-NiSQ, and also
computing time at the supercomputer Vega at IZUM, Slovenia. The numerical
calculation were partly carried out at the facilities of the Wroclaw Centre
for Networking and Supercomputing. Our TEBD code was written using ITensors
Library in Julia [51].

## References

  * 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, Phys. Rev. X 10, 011042 (2020).
  * Scherg  _et al._ [2021] S. Scherg, T. Kohlert, P. Sala, F. Pollmann, B. Hebbe Madhusudhana, I. Bloch, and M. Aidelsburger, Observing non-ergodicity due to kinetic constraints in tilted Fermi-Hubbard chains, Nature Communications 12, 4490 (2021).
  * Kohlert  _et al._ [2023] T. Kohlert, S. Scherg, P. Sala, F. Pollmann, B. Hebbe Madhusudhana, I. Bloch, and M. Aidelsburger, Exploring the regime of fragmentation in strongly tilted Fermi-Hubbard chains, Phys. Rev. Lett. 130, 010201 (2023).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * Alet and Laflorencie [2018] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics, Comptes Rendus. Physique 19, 498 (2018).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Sierant  _et al._ [2024] P. Sierant, M. Lewenstein, A. Scardicchio, L. Vidmar, and J. Zakrzewski, Many-body localization in the age of classical computing (2024), arXiv:2403.07111 [cond-mat.dis-nn] .
  * Schulz  _et al._ [2019] M. Schulz, C. A. Hooley, R. Moessner, and F. Pollmann, Stark many-body localization, Phys. Rev. Lett. 122, 040606 (2019).
  * Nieuwenburg  _et al._ [2019] E. Nieuwenburg, Y. Baum, and G. Rafael, From Bloch oscillations to many-body localization in clean interacting systems, Proceedings of the National Academy of Sciences of the United States of America 116, 9269 (2019).
  * Taylor  _et al._ [2020] S. R. Taylor, M. Schulz, F. Pollmann, and R. Moessner, Experimental probes of Stark many-body localization, Phys. Rev. B 102, 054206 (2020).
  * Bhakuni  _et al._ [2020] D. S. Bhakuni, R. Nehra, and A. Sharma, Drive-induced many-body localization and coherent destruction of Stark many-body localization, Phys. Rev. B 102, 024201 (2020).
  * Zhang  _et al._ [2021] L. Zhang, Y. Ke, W. Liu, and C. Lee, Mobility edge of Stark many-body localization, Phys. Rev. A 103, 023323 (2021).
  * Morong  _et al._ [2021] W. Morong, F. Liu, P. Becker, K. S. Collins, L. Feng, A. Kyprianidis, G. Pagano, T. You, A. V. Gorshkov, and C. Monroe, Observation of Stark many-body localization without disorder, Nature 599, 393 (2021).
  * Wang  _et al._ [2021] Y.-Y. Wang, Z.-H. Sun, and H. Fan, Stark many-body localization transitions in superconducting circuits, Phys. Rev. B 104, 205122 (2021).
  * Doggen  _et al._ [2021] E. V. H. Doggen, I. V. Gornyi, and D. G. Polyakov, Stark many-body localization: Evidence for Hilbert-space shattering, Phys. Rev. B 103, L100202 (2021).
  * Gunawardana and Buča [2022] T. M. Gunawardana and B. Buča, Dynamical l-bits and persistent oscillations in Stark many-body localization, Phys. Rev. B 106, L161111 (2022).
  * Zisling  _et al._ [2022] G. Zisling, D. M. Kennes, and Y. Bar Lev, Transport in Stark many-body localized systems, Phys. Rev. B 105, L140201 (2022).
  * Wei  _et al._ [2022] X. Wei, X. Gao, and W. Zhu, Static and dynamical Stark many-body localization transition in a linear potential, Phys. Rev. B 106, 134207 (2022).
  * Sala  _et al._ [2020] P. Sala, T. Rakovszky, R. Verresen, M. Knap, and F. Pollmann, Ergodicity breaking arising from Hilbert space fragmentation in dipole-conserving Hamiltonians, Phys. Rev. X 10, 011047 (2020).
  * Khemani  _et al._ [2020] V. Khemani, M. Hermele, and R. Nandkishore, Localization from Hilbert space shattering: From theory to physical realizations, Phys. Rev. B 101, 174204 (2020).
  * Yang  _et al._ [2020] Z.-C. Yang, F. Liu, A. V. Gorshkov, and T. Iadecola, Hilbert-space fragmentation from strict confinement, Phys. Rev. Lett. 124, 207602 (2020).
  * Moudgalya and Motrunich [2022] S. Moudgalya and O. I. Motrunich, Hilbert space fragmentation and commutant algebras, Phys. Rev. X 12, 011050 (2022).
  * Moudgalya  _et al._ [2022] S. Moudgalya, B. A. Bernevig, and N. Regnault, Quantum many-body scars and Hilbert space fragmentation: a review of exact results, Reports on Progress in Physics 85, 086501 (2022).
  * Feng and Skinner [2022] X. Feng and B. Skinner, Hilbert space fragmentation produces an effective attraction between fractons, Phys. Rev. Res. 4, 013053 (2022).
  * Deutsch [1991] J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016), 1509.06411 .
  * Łydżba  _et al._ [2024] P. Łydżba, P. Prelovšek, and M. Mierzejewski, Local integrals of motion in dipole-conserving models with Hilbert space fragmentation, (2024), arXiv:2401.17097 [cond-mat.str-el] .
  * Rakovszky  _et al._ [2020] T. Rakovszky, P. Sala, R. Verresen, M. Knap, and F. Pollmann, Statistical localization: From strong fragmentation to strong edge modes, Phys. Rev. B 101, 125126 (2020).
  * Nandy  _et al._ [2023] S. Nandy, Z. Lenarčič, E. Ilievski, M. Mierzejewski, J. Herbrych, and P. Prelovšek, Spin diffusion in a perturbed isotropic Heisenberg spin chain, Phys. Rev. B 108, L081115 (2023).
  * 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, Phys. Rev. Lett. 125, 245303 (2020).
  * Nandkishore and Hermele [2019] R. M. Nandkishore and M. Hermele, Fractons, Annual Review of Condensed Matter Physics 10, 295 (2019).
  * Gromov  _et al._ [2020] A. Gromov, A. Lucas, and R. M. Nandkishore, Fracton hydrodynamics, Phys. Rev. Res. 2, 033124 (2020).
  * Burchards  _et al._ [2022] A. G. Burchards, J. Feldmeier, A. Schuckert, and M. Knap, Coupled hydrodynamics in dipole-conserving quantum systems, Phys. Rev. B 105, 205127 (2022).
  * Guo  _et al._ [2022] J. Guo, P. Glorioso, and A. Lucas, Fracton hydrodynamics without time-reversal symmetry, Phys. Rev. Lett. 129, 150603 (2022).
  * Głódkowski  _et al._ [2023] A. Głódkowski, F. Peña Benítez, and P. Surówka, Hydrodynamics of dipole-conserving fluids, Phys. Rev. E 107, 034142 (2023).
  * Zechmann  _et al._ [2022] P. Zechmann, A. Bastianello, and M. Knap, Tunable transport in the mass-imbalanced Fermi-Hubbard model, Phys. Rev. B 106, 1 (2022), 2205.12970 .
  * Mori [1965] H. Mori, Transport, collective motion, and Brownian motion, Progress in Theoretical Physics 33, 423 (1965).
  * Forster [1975] D. Forster, Hydrodynamic Fluctuations, Broken Symmetry and Correlation Functions, Taylor - Francis Group, CRC Press (1975).
  * Götze and Wölfle [1972] W. Götze and P. Wölfle, Homogeneous dynamical conductivity of simple metals, Phys. Rev. B 6, 1226 (1972).
  * Jung and Rosch [2007] P. Jung and A. Rosch, Lower bounds for the conductivities of correlated quantum systems, Phys. Rev. B 75, 245104 (2007).
  * Kadanoff and Martin [1963] L. P. Kadanoff and P. C. Martin, Hydrodynamic equations and correlation functions, Annals of Physics 24, 419 (1963).
  * Bonča and Jaklič [1995] J. Bonča and J. Jaklič, Spin diffusion of the t-J model, Phys. Rev. B 51, 16083 (1995).
  * Steinigeweg and Gemmer [2009] R. Steinigeweg and J. Gemmer, Density dynamics in translationally invariant spin- 1 2 chains at high temperatures: A current-autocorrelation approach to finite time and length scales, Phys. Rev. B 80, 40 (2009).
  * Long  _et al._ [2003] M. W. Long, P. Prelovšek, S. El Shawish, J. Karadamoglou, and X. Zotos, Finite-temperature dynamical correlations using the microcanonical ensemble and the Lanczos algorithm, Phys. Rev. B 68, 235106 (2003).
  * Prelovšek and Bonča [2013] P. Prelovšek and J. Bonča, Ground state and finite temperature Lanczos methods, in  _Strongly Correlated Systems - Numerical Methods_, edited by A. Avella and F. Mancini (Springer, Berlin, 2013).
  * Herbrych  _et al._ [2012] J. Herbrych, R. Steinigeweg, and P. Prelovšek, Spin hydrodynamics in the $s=\frac{1}{2}$ anisotropic Heisenberg chain, Phys. Rev. B 86, 115106 (2012).
  * [48] See Supplemental Material for the anlysis the relaxation of dipole moment and related current correlations, as well as for results for other model parameters, involving also the additional numerical approach using the boundary-driven open systems.
  * Mierzejewski and Prelovšek [2010] M. Mierzejewski and P. Prelovšek, Nonlinear current response of an isolated system of interacting fermions, Phys. Rev. Lett. 105, 1 (2010).
  * Nandy  _et al._ [2024] S. Nandy, J. Herbrych, Z. Lenarčič, A. Głódkowski, P. Prelovšek, and M. Mierzejewski, Emergent dipole moment conservation and subdiffusion in tilted chains, Phys. Rev. B 109, 115120 (2024).
  * Fishman  _et al._ [2022] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022).
  * Prelovšek  _et al._ [2022] P. Prelovšek, S. Nandy, Z. Lenarčič, M. Mierzejewski, and J. Herbrych, From dissipationless to normal diffusion in the easy-axis Heisenberg spin chain, Phys. Rev. B 106, 245104 (2022).

a

Supplemental Material:
Einstein relation for subdiffusive relaxation in Stark chains

P. Prelovšek1, S. Nandy1, M. Mierzejewski2

1Department of Theoretical Physics, J. Stefan Institute, SI-1000 Ljubljana,
Slovenia

2Institute of Theoretical Physics, Faculty of Fundamental Problems of
Technology,
Wrocław University of Science and Technology, 50-370 Wrocław, Poland

In the Supplemental Material we analyze the relaxation of dipole moment and
related current correlations, We present also results for other model
parameters, involving the additional numerical approach using the boundary-
driven open systems.

##  S1 Relaxation of the dipole moment and current correlations

Studying the Stark chain in the main text we argue that for large $L$, nonzero
$F$ and small $\omega$, the dominating contribution to the memory function,
Eq. (7), comes from the dipolar currents, $\phi_{D}(\omega)$, instead of the
normal current correlations, $\phi_{N}(\omega)$. The latter expectation
originates from the emergent conservation of the dipole moment in macroscopic
systems, $L\to\infty$ [50]. Below we provide numerical results which support
this claim. However, we note that in a finite system one gets
$\phi_{N}(\omega=0)=0$ just due to open boundary conditions (OBC) simply
because a finite system with OBC cannot host a steady current. In order to
disentangle the influence of OBC in a finite system from the conservation of
the dipole moment we study the relaxation of the dipole moment
$P=\sum_{l}(l-L/2)n_{l}$. To this end we calculate the correlation function

$\phi_{P}(\omega)=(P|({\cal
L}-\omega)^{-1}|P)=\frac{-\chi^{0}_{P}}{\omega+N(\omega)}.$ (S1)

In equation (S1) we introduced the corresponding MF

$N(\omega)=\frac{L}{\chi^{0}_{P}}\tilde{\phi}_{N}(\omega),$ (S2)

where $\chi^{0}_{P}=\langle P^{2}\rangle$ and formally MF can be expressed as
[39],

$\tilde{\phi}_{N}(\omega)=({\cal L}P|({\cal\tilde{L}}-\omega)^{-1}|{\cal
L}P)/L=(J_{N}|({\cal\tilde{L}}-\omega)^{-1}|J_{N})/L,$ (S3)

where we used the identity ${\cal L}P=-iJ_{N}$. In contrast to Eq. (2) in the
main text that is valid in the hydrodynamic (perturbative) regime, there is no
small parameter in Eq. (S3) thus $\cal\tilde{L}$ contains additional
projectors [39],

${\cal\tilde{L}}={\cal QLQ},\quad{\cal Q}=1-{\cal P},\quad{\cal
P}=\frac{|P)(P|}{(P|P)}.$ (S4)

For this reason in a finite system $\tilde{\phi}_{N}(\omega)$ in Eq. (S3) is
not equal to $\phi_{N}(\omega)=(J_{N}|({\cal L}-\omega)^{-1}|J_{N})/L$ but
both quantities should merge for $L\to\infty$ due to the conservation of $P$.

Figure S1: (a) The current correlation function
$\mathrm{Im}\tilde{\phi}_{N}(\omega)$, as extracted form dipole-moment
correlations $\phi_{P}(\omega)$, calculated for Stark model using MCLM on
$L=28$ chain, for different fields $F=0-1.0$, (b) directly evaluated current
correlations $\mathrm{Im}\phi_{N}(\omega)$ in the same system.

Due to the presence of projectors, it is hard to determine
$\tilde{\phi}_{N}(\omega)$ from Eq. (S3). However, this quantity can be
obtained also from the correlations of the dipole moment via Eqs. (S1) and
(S2). In order to confirm the relation between $\tilde{\phi}_{N}(\omega)$ and
$\phi_{N}(\omega)$ and related conjectures, we numerically investigate the
Stark model with the same parameters $V=V^{\prime}=1$ for different fields
$F\leq 1$. On one hand, using MCLM on the $L=28$ system with OBC, we evaluate
$\phi_{P}(\omega)$ and extract the corresponding
$\mathrm{Im}\tilde{\phi}_{N}(\omega)$. On the other hand, we calculate within
the same system directly $\mathrm{Im}\phi_{N}(\omega)$. Numerical results
presented in Figs. S1(a) and S1(b) confirm that at $F>F^{*}\sim 0.4$ (for
considered $L=28$) both approaches indeed yield (at least qualitatively)
$\tilde{\phi}_{N}(\omega)\sim\phi_{N}(\omega)$, up to finite-size effects. As
expected, the disagreement occurs at $F\sim 0$ since $\phi_{N}(\omega\to 0)$
vanishes due to OBC, while $\tilde{\phi}_{N}(\omega\to 0)$ does not. Due to
the same argument the qualitative discrepancy in the range $F<F^{*}(L)$ can be
attributed to finite-size limitations.

The main conclusion is that $P$ behaves as a conserved quantity and, at the
same time, $\mathrm{Im}\phi_{N}(\omega\to 0)\to 0$, so that in Eq. (7) in the
main text one can neglect the latter contribution in the hydrodynamic regime
$q\to 0$. Consequently, the relaxation is then dominated by the SD
contribution, i.e., $M_{q}(\omega\sim 0)\sim q^{4}\phi_{D}(\omega\to
0)/\chi^{0}_{0}$. Based on results in Fig. S1 this conclusion holds true at
least for $F>F^{*}\sim 0.4$, whereby we expect that $F^{*}$ vanishes in the
thermodynamic limit.

##  S2 Subdiffusion coefficient from boundary driven open systems

As an additional technique to compute the subdiffusion constant, we study an
open chain described by Eq. (6), driven via boundary Lindblads operators
inducing a weak particle current by creating a bias $\mu$ at the edges of the
system. The master equation governing the time-evolution of the density matrix
is given by

$\displaystyle\partial_{t}\rho=-i[H,\rho]+\hat{\mathcal{D}}\rho$ (S5)

where $\hat{\mathcal{D}}$ denotes the dissipator operator written in terms of
the Lindblads as $\hat{\mathcal{D}}\rho=\sum_{k}L_{k}\rho
L_{k}^{\dagger}-\frac{1}{2}\\{L_{k}^{\dagger}L_{k},\rho\\}$. The details of
the Lindblads employed here are discussed in [30, 50, 52]. To evolve the
density matrix towards the steady state $\rho_{ss}$, we use the time-evolving
block decimation (TEBD) for vectorized density matrices. In particular, we use
the fourth-order TEBD with a time step $dt=0.2$, bond dimension $\chi\sim
140$, and weak bias $\mu=0.01$.

The method allows us to establish a nonequilibrium steady state (NESS) for
which the normalized current $I/\mu$ and the spatial profile of particles
$\langle\tilde{n}_{l}\rangle$ can be easily computed. One remarkable advantage
of this technique is that the NESS profile carries the signature of the nature
of particle transport [50]. In particular, for system with both particle-
number and dipole-moment conservation, the hydrodynamic equation of motion
$\partial_{t}n+D\partial^{4}_{x}n=0$ implies the NESS particle profile of the
form $\langle\tilde{n}_{l}\rangle=ax+bx^{3},x=2\ell/L-1$.

Figure S2: Steady-state density profiles $\langle\tilde{n}_{l}\rangle/\mu$
obtained by the boundary-driven open-system approach, obtained for different
fields $F$ on $L=40$ system. Inset shows the ratio of profile coefficients
$|a/b|$ vs. $F$ for two system sizes $L$.

We further present and analyse the NESS results for the case
$V=2,V^{\prime}=0$, which is at $F=0$ equivalent to the isotropic Heisenberg
model with the corresponding particle superdiffusion [30], while for $F>0$ the
SD is expected with results similar as presented in the main text for
$V=V^{\prime}=1$. Fig. S2 clearly demonstrates how the normalized particle
profile changes with increasing $F$. The profile at $F=0$ is indeed a
signature of particle superdiffusion where
$\langle\tilde{n}_{l}\rangle\sim\arcsin(x)$ (see [30] for details). As evident
in Fig. S2, increasing $F>0$ leads to qualitative change of the profile.
Whereas it is close to linear (corresponding to normal diffusion) for small
finite $F\lesssim 0.4$, eventually the profile becomes cubic for larger
$F\gtrsim 0.6$, being a direct signature of the crossover to SD transport.
Here, we refer to reference [50], for further details regarding the fitting of
coefficients $a,b$. Finally, for large $F$ with the SD transport we establish
the ratio $|a/b|\sim 3$, as the consequence of vanishing derivative
$\delta_{x}(ax+bx^{3})|_{x=\pm 1}=0$ at the edges, being evident in the main
panel of Fig. S3. The inset shows that indeed the ratio $|a/b|\to 3$ with
increasing $F$ and also with increasing $L$. We argue that the edge flattness
carries the signature of the fact that in SD systems the normal current
becomes zero (or negligibly small) for $F>0$ in the $L\to\infty$ limit. Using
the relation $|a/b|=3$, the SD coefficient can be then expressed directly with
the obtained NESS current as $D_{S}=IL^{3}/(12\mu)$.

##  S3 Results for other model parameters

Figure S3: Subdiffusion coefficient $D_{S}(F)$, as calculated: from rates
$\Gamma_{q_{1}}$, obtained via time-evolution (TE) and via MCLM on $L=24$
chain, and directly from Einstein relation (ER) and $\phi_{D}(\omega)$ on
$L=28$ chain, as well via the boundary-driven open systems (OS) on $L\leq 50$
chains.

We present in Fig. S3 the summary of results for the SD coefficient $D_{S}$,
as obtained by different numerical approaches for parameters
$V=2,V^{\prime}=0$, corresponding to the isotropic Heisenberg model subject to
constant field (magnetic field gradient). Besides the methods presented and
used in the main text, i.e., the time-evolution (TE) and MCLM where $D_{S}$ is
extracted from relaxation rates $D_{S}=\Gamma_{q_{1}}/q_{1}^{4}$ for the
smallest $q_{1}=2\pi/L$ at $L=24$, and the direct evaluation using the
Einstein relation (ER) on $L=28$ chain, we include now also the NESS results
for the boundary-driven systems with $L\lesssim 50$. Results are presented for
fields $F>F^{*}(L)\sim 0.4$ where the SD transport is dominant for numerically
accessible $L$. Employed methods have also different numerical limitations at
large $F>1$. While the NESS approach is limited by smallness of the current
$I$, the MCLM and Einstein relation are limited by the frequency resolution
and the influence of OBC, so the time-evolution approach has the advantage of
smallest reachable $\Gamma_{q}$ and consequently of largest reachable
$F\lesssim 2$ due to PBC in this case. The overall conclusion following from
Fig. S3 is that results by different methods are even quantitatively well
consistent, revealing again the strong/exponential dependence of $D_{S}$ on
$F$. On the other hand, by comparing with Fig. 2 (b) in the main text, we also
notice that $D_{S}$ only weakly depends on the interaction parameters $V$ and
$V^{\prime}$. For stronger interactions we observe an increase $D_{S}$ at
larger $F$ since the interaction suppresses localization, while at weaker $F$
the trend might be opposite, but further studies are needed to establish such
relations.