# Anomalous diffusion in quantum system driven by heavy-tailed stochastic
processes

Chenyue Guo<EMAIL_ADDRESS>Wilczek Quantum Center and Key Laboratory
of Artificial Structures and Quantum Control, School of Physics and Astronomy,
Shanghai Jiao Tong University, Shanghai 200240, China

###### Abstract

In this paper, we study a stochastically driven non-equilibrium quantum system
where the driving protocols consist of hopping and waiting processes. The
waiting times between two hopping processes satisfy a heavy-tailed
distribution. By calculating the squared width of the wavepackets, our
findings demonstrate the emergence of various anomalous transport phenomenons
when the system remains unchanged within the heavy-tailed regime, including
superdiffusive, subdiffusive, and standard diffusive motion. Only subdiffusion
occurs when the system has evolved during the waiting process. All these
transport behaviors are accompanied by a breakdown of ergodicity, highlighting
the complex dynamics induced by the stochastic driving mechanism.

##  I introduction

The continuous-time random walk (CTRW) model [Montroll and Weiss, 1965; Weiss,
1994; Klafter and Sokolov, 2011] is a mathematical framework used to describe
the dynamic behavior of systems, particularly the diffusion process of
particles in complex systems and non-uniform media [Scher and Montroll, 1975;
Bouchaud and Georges, 1990; Barkai  _et al._ , 2000; Weiss, 2013]. This model
applies to multiple fields such as physics [Shlesinger, 1974; Pfister and
Scher, 1977; Geisel  _et al._ , 1985; Solomon  _et al._ , 1993; Ohtsuki and
Kawarabayashi, 1997], finance [Mandelbrot, 1963; Masoliver  _et al._ , 2003;
Scalas, 2006], and many others [Helmstetter and Sornette, 2002; Mega  _et al._
, 2003; H$\ddot{\rm{o}}$fling and Franosch, 2013], because it can handle
discontinuities in time and space as well as long-range dependencies. The CTRW
model incorporates two independent stochastic variables [Scher and Montroll,
1975; Metzler and Klafter, 2000]: jump length $x$ and waiting time $t$, as
described by the probability distribution functions (PDFs) $\lambda(x)$ and
$\psi(t)$, respectively. They characterize a particle’s waiting time $t$ at
its initial position, then followed by a jump of length $x$, after which the
process restarts. When the PDFs $\lambda(x)$ and $\psi(t)$ exhibit finite
first and second moments, such as $\lambda(x)$ and $\psi(t)$ follow a
Poissonian and Gaussian distribution, the motion of system corresponds to
normal diffusion [Klafter and Sokolov, 2011]. Conversely, if the waiting time
PDF $\psi(t)$ adheres to a power-law distribution with a long tail, such that
$\psi(t)\propto 1/t^{(1+\alpha)}$ with $0<\alpha<1$, the first moment,
$\langle\tau\rangle=\int_{0}^{\infty}t\psi(t)dt$, diverges. The ensemble-
averaged mean squared displacements satisfy $\langle x^{2}(t)\rangle_{ens}\sim
t^{\alpha}$, indicating a subdiffusive motion [Metzler and Klafter, 2000;
Sokolov  _et al._ , 2002]. Additionally, the ensemble-averaged mean squared
displacements do not equal the time-averaged mean squared displacements,
indicating ergodicity breaking [Bouchaud, 1992; Barkai and Margolin, 2004;
Margolin and Barkai, 2005; Lubelski  _et al._ , 2008; He  _et al._ , 2008].

A profusion of equilibrium quantum systems, including generic one-dimensional
integrable models [Zotos  _et al._ , 1997; Batchelor, 2005; Miller  _et al._ ,
2020], chaotic systems [Karrasch  _et al._ , 2014;
$\rm{\check{Z}nidari\check{c}}$  _et al._ , 2016; Blake  _et al._ , 2017;
Friedman  _et al._ , 2020] or disordered models [Anderson, 1958; Agarwal  _et
al._ , 2015; Bar Lev  _et al._ , 2015; Luitz  _et al._ , 2016] typically
exhibit ballistic, diffusive or subdiffusive (localization) motion. The
dynamics of non-equilibrium have garnered significant interest due to
substantial experimental progress [Bloch  _et al._ , 2008; Chang  _et al._ ,
2018]. By employing methods such as quenching [Polkovnikov  _et al._ , 2011;
Luca D’Alessio and Rigol, 2016], driving [Struck  _et al._ , 2012; Cai  _et
al._ , 2017; Harper  _et al._ , 2020; Guo  _et al._ , 2023] or coupling to
environment [Sieberer  _et al._ , 2013; Cai and Barthel, 2013; Marino and
Diehl, 2016; Lee  _et al._ , 2022; Erbanni  _et al._ , 2023], one can drive
the system far from equilibrium and then analyze the transport phenomena
emerged in the non-equilibrium system. Notably, recent studies [Wang  _et al._
, ; $\rm{\check{Z}nidari\check{c}}$, 2024] have unveiled an innovative
approach that demonstrates the emergence of superdiffusion by coupling quantum
lattice systems with multi-site dephasing dissipation. Ref [Wang  _et al._ , ]
provided an insightful theoretical arguement using the framework of
L$\rm{\acute{e}}$vy walk. The L$\rm{\acute{e}}$vy walk (L$\rm{\acute{e}}$vy
flight) model, with a finite (infinite) velocity of a random walker, is known
for generating anomalous diffusion due to jump determined by the heavy-tailed
distribution [Zaburdaev  _et al._ , 2015]. And this model is often regarded as
the theoretical background for explaining superdiffusion [Wang  _et al._ , ;
$\rm{\check{Z}nidari\check{c}}$, 2024; Zaburdaev  _et al._ , 2015; Schuckert
_et al._ , 2020; Joshi  _et al._ , 2022].

In this paper, we study a stochastically driven non-equilibrium quantum
system. Unlike previous studies on stochastically driven systems [Nandy  _et
al._ , 2017; Halimeh  _et al._ , 2018; Pan  _et al._ , 2020; Kos  _et al._ ,
2021; Cai, 2022], the driving time between two hopping processes satisfies the
power-law PDF, $\psi(t)\propto 1/t^{(1+\alpha)}$ with $0<\alpha<1$, a typical
heavy-tailed distribution. It is noteworthy that the case is the quantum
analog of the CTRW model. The toy model is composed of two parts, one depicted
by Hamiltonian. (3), which exhibits ballistic transport behavior distinct from
the diffusion seen in classical random walks (Brownian motion) due to phase
coherence and quantum interference effects. We are particularly interested in
the transport behavior of the quantum system after introducing an evolution
time $t$ distributed as a heavy-tailed distribution. This process is governed
by another Hamiltonian. (2). Here, we use the width of the wavepacket to
characterize the system’s dynamics. For the case $V=0$ in Hamiltonian. (2), we
observe that the ensemble-averaged squared width of the wavepackets grows with
time as $t^{2\alpha}$, corresponding to subdiffusive, superdiffusive and
normal diffusive motion, with $\alpha$ varying in different intervals.
Furthermore, we calculate both the time-averaged squared width of the
wavepackets and its ensemble average. It grows with time as $t^{1+\alpha}$.
These two types of averages do not converge, leading to nonergodicity. For the
case where $V\neq 0$, we observe that the ensemble-averaged (time-averaged)
squared width of the wavepackets grows with time as $t^{\alpha}$ ($t^{1}$),
corresponding to subdiffusive (diffusive) motion, indicating the nonergodicity
reoccurs.

The structure of the paper is organized as follows: Sec. II presents the
single-particle models defined in a one-dimensional (1D) lattice and methods
used in this paper. In Sec. III, we explore the transport behavior of
ensemble-averaged and time-averaged squared width of the wavepackets. Sec. IV
provides the conclusion and outlook.

Figure 1: (Color online). Example of a typical trajectory of the stochastic
driven quantum system. Control parameter $f(t)$ in the Hamiltonian. (1) equals
to $1$ or $0$ corresponds to the duration of $\hat{H}_{0}$ and $\hat{H}_{1}$,
respectively. The time for waiting process follows the power law distribution.
Figure 2: (Color online). Illustration to ensemble-averaged and time-averaged
method in this paper. (a) $\langle\cdots\rangle_{ens}$ represents the quantity
over different waiting time PDF trajectories. (b) $\langle\cdots\rangle_{T}$
represents the quantity over different time parts within a single waiting time
PDF trajectory.

##  II Model and method

We consider a 1D single-particle model with Hamiltonian defined as:

$\hat{H}(t)=f(t)\hat{H}_{0}+[1-f(t)]\hat{H}_{1},$ (1)

where

$\displaystyle\hat{H}_{0}$
$\displaystyle=\sum_{i}{(-1)^{i}V\hat{c}_{i}^{\dagger}\hat{c}_{i}},$ (2)
$\displaystyle\hat{H}_{1}$
$\displaystyle=\sum_{i}{-J\hat{c}_{i}^{\dagger}\hat{c}_{i+1}+h.c.}.$ (3)

Where $\hat{H}_{1}$ and $\hat{H}_{0}$, including the free fermion hopping term
and on-site potential term, correspond to the process of hopping and waiting
governed by $f(t)$, respectively. $\hat{c}_{i}$ ($\hat{c}_{i}^{\dagger}$)
denotes the annihilation (creation) operator of a spinless fermion on site
$i$. The parameter $J$ represents the amplitude for nearest-neighbor (NN)
single-particle hopping and $V$ is the on-site chemical potential. Open
boundary condition is adapted.

Starting from an initial state $\phi(t=0)=\delta(i=L/2)$ where a particle is
located at site $i=L/2$, the wave function $\phi(t)$ evolves according to the
Schr$\ddot{\rm{o}}$dinger eqation $i\hbar\partial\hat{H}_{0}(t)/\partial
t=\hat{H}_{0}(t)\phi(t)$, with the time $t$ for the waiting process following
a power-law distribution $\psi(t)\propto 1/t^{1+\alpha}$, with $\alpha=0.65$
in the main text [it is applicable for other $\alpha$ values in Appendix A].
Then the particle starts hopping governed by the Hamiltonian $\hat{H}_{1}$
with $i\hbar\partial\hat{H}_{1}(t)/\partial t=\hat{H}_{1}(t)\phi(t)$, and the
process is renewed. For simplification, the hopping time in the quantum model
is quantified as lasting one second.

To quantify the transport behavior emerged in the driven quantum model, we
introduce the squared width of the wavepacket $W^{2}(t)$, defined as:

$W^{2}(t)=\sum_{i}\rho_{i}(t)[i-\bar{x}(t)]^{2},$ (4)

where the density distribution
$\rho_{i}(t)=\langle\phi(t)|\hat{c}_{i}^{\dagger}\hat{c}_{i}|\phi(t)\rangle$
and the center of matter (COM) of the wave packet
$\bar{x}(t)=\sum_{i}i\rho_{i}(t)$. We can calculate the ensemble-averaged
squared width of the wavepacket $\langle W^{2}(t)\rangle_{ens}$ over different
waiting time PDF trajectories, as illustrated in Fig. 2(a).

In order to calculate the time-averaged squared width of the wavepackets
$\langle W^{2}(t)\rangle_{T}$ for a given waiting time PDF trajectory, we
consider the time evolution of the wave function over time $\tau$ at different
initial position $t_{i}^{\prime}$, i.e.,
$|\phi(t_{i}^{\prime}+\tau)\rangle=\mathcal{\hat{T}}exp[-i\int_{t_{i}^{\prime}}^{\tau+t_{i}^{\prime}}\hat{H}(t^{\prime\prime})dt^{\prime\prime}]|\phi(t_{i}^{\prime})\rangle$.
Where $\mathcal{\hat{T}}$ is the time-ordering operator.
$\hat{H}(t^{\prime\prime})=\hat{H}_{0}$ or $\hat{H}_{1}$ is determined by
whether the time $t^{\prime\prime}$ lies in the waiting or hopping process.
For every finite total evolution time $T$, $t_{i}^{\prime}\in[0,T-t]$ with
$i=1,2,..,\frac{T-t}{dt_{i}^{\prime}}$, as illustrated in Fig. 2(b). It is
noteworthy that the width of the wavepacket has spread over many lattice sites
after the evolution time interval $\tau$, making it necessary to reposition
the particle at $\delta(i=L/2)$ at each time $t_{i}^{\prime}$, i.e.,
$|\phi(t_{i}^{\prime})\rangle=\delta(i=L/2)$. we can calculate a series of the
squared width of the wavepackets defined in eq. (4) for different
$t_{i}^{\prime}$. The time-averaged squared width of the wavepackets $\langle
W^{2}(t)\rangle_{T}$ is then obtained by averaging these results.
Additionally, for every finite total evolution time $T$, this quantity is also
a random variable. Thus we consider its ensemble average $\langle\langle
W^{2}(t)\rangle_{T}\rangle_{ens}$ over different waiting time PDF
trajectories.

Figure 3: (Color online). (a) Width of the wavepacket under three different
waiting-time PDFs. (b) Ensemble-averaged square width of the wavepackets. (c)
Ensemble time-averaged square of width of the wavepackets. $L=8000$ for
(a)-(b). $L=1200$ for (c). The solid (dashed) black (red) line with circular
markers represents the real (fitted) data. The parameter $\Delta t=0.1J^{-1}$,
$V=0$ and $\psi(t)\propto 1/t^{(1+\alpha)}$ with $\alpha=0.65$. Figure 4:
(Color online). (a) Width of the wavepacket under three different waiting-time
PDFs. (b) Ensemble-averaged square width of the wavepackets. (c) Ensemble
time-averaged square of width of the wavepackets. $L=2500$ for (a)-(b).
$L=600$ for (c). The solid (dashed) black (red) line with circular markers
represents the real (fitted) data. The parameter $\Delta t=0.1J^{-1}$, $V=4J$
and $\psi(t)\propto 1/t^{(1+\alpha)}$ with $\alpha=0.65$.

##  III Results

###  III.1 The case for free fermion

In this section, we first consider the wave function remaining unchanged
during the waiting process, i.e., $V$ is always set to zero (zero-potential).
Specifically, we start from an initial state $\phi(t=0)=\delta(i=L/2)$, the
wave function $\phi(t)$ keeps unchanged [$\hat{H}_{0}=0$] for waiting time $t$
following a power-law distribution. Then the particle starts hopping governed
by the Hamiltonian $\hat{H}_{1}$, and the process is renewed.

Fig. 3(a) intuitively exhibits three trajectories of the squared width of the
wavepackets. We choose the discrete time step $\Delta t=0.1J^{-1}$ and system
size $L=8000$, both shown to be sufficiently accurate in Appendix B.
Considering $W(t)$ as a random variable dependent on different waiting-time
PDFs $\psi_{i}(t)$, we calculate the ensemble-averaged squared width of the
wavepackets $\langle W(t)^{2}\rangle_{ens}$ over $5\times 10^{4}$ waiting time
PDF trajectories, as shown as a log-log plot in Fig. 3(b). We observe that it
grows with time $t$ as a function of $t^{1.3155}$ ($\sim t^{2\alpha}$). Given
that the jumping process corresponds to a free fermion hopping model within
$W(t)\sim t$, ballistic transport behavior, we only restrict the particle
waiting for time $t_{i}$ based on the ballistic transport behavior. This
allows us to explain the result from classical point of view. Specifically,
the PDF of displacement $x$ satisfies $\lambda(x)=\delta(x-1)$ for ballistic
transport and $\psi(t)\propto 1/t^{(1+\alpha)}$ with $0<\alpha<1$. Their
Fourier and Laplace transforms read as

$\displaystyle\lambda(k)$ $\displaystyle\sim 1+ik-\frac{1}{2}k^{2},$ (5)
$\displaystyle\psi(u)$ $\displaystyle\sim 1-u^{\alpha}.$ (6)

In fact, any pair of PDFs leading to finite $\langle x^{2}\rangle$ and
$\langle t\rangle$ results in the same outcomes [Klafter  _et al._ , 1987].
Considering the variable $x$ and $t$ are independent of each other, so the PDF
to find a walker at position $x$ and time $t$ is determined by

$p(k,u)=\psi(u)\lambda(k)=1+ik-\frac{1}{2}k^{2}-u^{\alpha}.$ (7)

The Fourier-Laplace transform [Metzler and Klafter, 2000; Klafter  _et al._ ,
1987] is given by

$\displaystyle P(k,u)$
$\displaystyle=\frac{1-\psi(u)}{u}\frac{W_{0}(k)}{1-p(k,u)}$
$\displaystyle=\frac{2u^{\alpha-1}W_{0}(k)}{k(-2i+k)+2u^{\alpha}}.$ (8)

Where $W_{0}(k)$ denotes the Fourier transform of the initial condition
$W_{0}(x)=\delta(x-1)$. The time dependence of the moments of the
displacement,
$M_{n}(t)=\mathscr{L}^{-1}[(-i)^{n}\frac{\mathrm{d}^{n}P(k,u)}{\mathrm{d}k^{n}}]$
where $\mathscr{L}^{-1}$ is the inverse Laplace transform. Thus, we can derive
the expressions for the first and second moments,

$\displaystyle\langle x(t)\rangle$
$\displaystyle=M_{1}(t)=\frac{W_{0}(k)}{\Gamma(1+\alpha)}t^{\alpha},$
$\displaystyle\langle x^{2}(t)\rangle$
$\displaystyle=M_{2}(t)=\frac{W_{0}(k)}{\Gamma(1+\alpha)}t^{\alpha}+\frac{2W_{0}(k)}{\Gamma(1+2\alpha)}t^{2\alpha}.$
(9)

Then the mean squared displacement of $x$ reads

$\displaystyle\sigma^{2}[x(t)]$ $\displaystyle=\langle x^{2}(t)\rangle-\langle
x(t)\rangle^{2},$
$\displaystyle=\frac{W_{0}(k)}{\Gamma(1+\alpha)}t^{\alpha}+\left[\frac{2W_{0}(k)}{\Gamma(1+2\alpha)}-\frac{W_{0}^{2}(k)}{\Gamma^{2}(1+\alpha)}\right]t^{2\alpha}.$
(10)

In the large $t$ limit where $t^{2\alpha}\gg t^{\alpha}$, eq. (III.1) reduced
to $\sigma^{2}[x(t)]\sim t^{2\alpha}$. Therefore, the squared width of the
wavepacket $W^{2}(t)$ correspond to the $\sigma^{2}[x(t)]$, $\sim
t^{2\alpha}$.

On a parallel front, we calculate the time-averaged squared width of the
wavepacket $\langle W^{2}(t)\rangle_{T}$ and its ensemble average
$\langle\langle W^{2}(t)\rangle_{T}\rangle_{ens}$ over $5\times 10^{4}$
waiting time PDF trajectories, as shown in Fig. 3(c). Here, we choose the
system size $L=1200$. We observe the quantity $\langle\langle
W^{2}(t)\rangle_{T}\rangle_{ens}$ grows with time as a function of
$t^{1.6387}$ $(\sim t^{1+\alpha})$, is not consistent with $\langle
W^{2}(t)\rangle_{ens}$, indicating nonergodicity. This directly corresponds to
the classical case where $\langle\delta(x)\rangle\neq 0$ and $\langle\langle
W^{2}(t)\rangle_{T}\rangle_{ens}\sim t^{\alpha}t=t^{1+\alpha}$ as discussed in
Ref [He  _et al._ , 2008].

###  III.2 The effect of potential

Figure 5: (Color online). Ensemble-averaged square of width of the wavepackets
over $10^{4}$ samples. The solid (dashed) black (red) line with circular
markers represents the real (fitted) data. The parameter $L=2500$, $\Delta
t=0.01J^{-1}$ and $V(t)\propto 1/V^{(1+\alpha)}$ with $\alpha=0.65$.

In this section, we consider the effect of the on-site potential $V$. For
simplicity, we set $V$ as a time-independent constant (constant-potential).
Specifically, we start from an initial state $\phi(t=0)=\delta(i=L/2)$, the
wave function $\phi(t)$ evolves according to
$i\hbar\partial\hat{H}_{0}(t)/\partial t=\hat{H}_{0}(t)\phi(t)$ for waiting
time $t$ following a power-law distribution. Then the particle starts hopping
governed by the Hamiltonian $\hat{H}_{1}$ with
$i\hbar\partial\hat{H}_{1}(t)/\partial t=\hat{H}_{1}(t)\phi(t)$, and the
process is renewed.

Fig. 4(a) illustrates three particle trajectories with the potential amplitude
$V$ set to $4J$. The ensemble-averaged squared width of the wavepackets
$\langle W(t)^{2}\rangle_{ens}$ is shown as a log-log plot in Fig. 4(b). We
observe that it grows with time $t$ as a function of $t^{0.6508}$ ($\sim
t^{\alpha}$) which is different from the case mentioned in Sec. III.1. In this
case, the unitary time evolution operator during the waiting process is given
by $e^{-iH_{0}t}=e^{\pm iVt}$ [$H_{0}$ is diagonal], where the absolute value
of global phase factor $Vt$ can become significantly large due to the PDF
$\psi(t)$ satisfied a power law distribution. To determine which key element,
$V$ or $t$, is more significant, we make $V$ follow the same power law
distribution $\psi(V)\propto 1/V^{(1+\alpha)}$ and set the waiting time to
four seconds to ensure the same global phase factor. The results are shown in
Fig. 5. It reveals that the ensemble-averaged squared width of the wavepackets
$\langle W(t)^{2}\rangle_{ens}\sim t$ and $\langle W(t)\rangle_{ens}\sim
t^{0.5}$ [normal diffusion], indicating that the chemical potential $V$, with
significant changes, is treated as white noise. Therefore, the presence of
anomalous diffusive transport is a consequence of the long-tailed distribution
of $t$ for the waiting process, which diminishes the amplitude and frequency
of the wavepacket changes. Hence we can interpret the results shown in Fig.
4(b) from a classical perspective. The dynamics of the wavepacket adhere to
normal diffusive behavior under the large amplitude of the phase factor,
analogous to classical diffusion motion (Brownian motion) where
$\lambda(x)=\frac{1}{2}\left[\delta(x-1)+\delta(x+1)\right]$. Thus, the eq.
(5) is rewritten as $\lambda(k)\sim 1-k^{2}/2$, and the first and second
moment are given by:

$\displaystyle\langle x(t)\rangle$ $\displaystyle=0,$ $\displaystyle\langle
x^{2}(t)\rangle$
$\displaystyle=M_{2}(t)=\frac{W_{0}(k)}{\Gamma(1+\alpha)}t^{\alpha}.$ (11)

Then the mean squared displacement of $x$, $\sigma^{2}[x(t)]=\langle
x^{2}(t)\rangle\sim t^{\alpha}$, which is consistent with our result.

Similarly, we calculate the time-averaged squared width of the wavepacket
$\langle W^{2}(t)\rangle_{T}$ and its ensemble average $\langle\langle
W^{2}(t)\rangle_{T}\rangle_{ens}$, as shown in in Fig. 4(c). Here, we choose
the system size $L=600$. It reveals that the quantity grows with time as
$t^{1.0102}$ $(\sim t)$, is not consistent with $\langle
W^{2}(t)\rangle_{ens}$, indicating nonergodicity once again. It directly
corresponds to the classical case where $\langle\delta(x)\rangle=0$ and
$\langle\langle W^{2}(t)\rangle_{T}\rangle_{ens}\sim t$ as discussed in Ref
[He  _et al._ , 2008].

##  IV Conclusion and outlook

In conclusion, our work explores the quantum analog of the CTRW model,
revealing various anomalous diffusive transport behaviors accompanied by the
breakdown of ergodicity. For the zero-potential case, we observe that the
ensemble-averaged squared width of the wavepackets grow with time as
$t^{2\alpha}$, encompassing subdiffusive, superdiffusive and the standard
diffusive motion. Additionally, the time-averaged square of width of
wavepackets and its ensemble average grow with time as $t^{1+\alpha}$. The
lack of convergence between these two averages indicates ergodicity breaking.
For the constant-potential case, the ensemble-averaged squared width of the
wavepackets grow with time as $t^{\alpha}$, corresponding to subdiffusive
motion. Meanwhile, $\langle\langle W^{2}(t)\rangle_{T}\rangle_{ens}$ grows
with time as $t^{1}$, signifying the recurrence of nonergodicity. The results
highlight the abundant dynamics induced by the stochastic driving mechanism.

The power-law distribution is known for its ability to describe complex
systems with a wide range of temporal scales. Up to this point, we have
somewhat artificially set the long-tailed waiting time in our models. This
raises an intriguing question: can we develop a model in which the time
correlation of physical quantities naturally conforms to a power-law
probability distribution? If successful, such a model would likely exhibit
richer dynamics compared to those governed by exponential decay interactions.

##  V Acknowledgments

We thank Igor M. Sokolov for helpful and insightful suggestions on the
continuous-time random walk. We also thank Zi Cai for the valuable discussions
and guidance throughout this work. This work is supported by the National Key
Research and Development Program of China (Grant No. 2020YFA0309000), NSFC of
China (Grant No.12174251), the Natural Science Foundation of Shanghai (Grant
No.22ZR142830), and the Shanghai Municipal Science and Technology Major
Project (Grant No. 2019SHZDZX01).

Figure 6: (Color online). (a)(d) Width of the wavepacket under three different
waiting-time PDFs. (b)(e) [(c)(f)] Ensemble-averaged [time-averaged] square
width of the wavepackets over $5\times 10^{4}$ samples. $L=6000$ for (a)(b)
and $L=1600$ for (c). $L=2500$ for (d)(e) and $L=600$ for (f). The solid
(dashed) black (red) line with circular markers represents the real (fitted)
data for the right panel. The parameter $\Delta t=0.1J^{-1}$, $V=0$ ($V=4J$)
for the upper (lower) panel and $\alpha=0.5$. Figure 7: (Color online).
Comparison between the width of the wavepacket for a waiting-time PDF with
different $\Delta t$ with $V=0$ (a) and $V=4J$ (b)-(c). (c) Comparison between
the ensemble-averaged square of width of the wavepackets over $5\times 10^{4}$
samples with different $\Delta t$. $L=8000$ for (a) and $L=2500$ for (b)-(c).
$\alpha=0.65$. Figure 8: (Color online). Comparison between the width of the
wavepacket for a waiting-time PDF with different system size $L$ with $V=0$
(a) and $V=4J$ (b). The parameter $\Delta t=0.1J^{-1}$ and $\alpha=0.65$.

##  Appendix A universal of numerical results

In the main text, we chose $\alpha=0.65$. The results are consistent for other
$\alpha$ values between 0 and 1. As shown in the upper (lower) panel of Fig.
6, where $\alpha=0.5$, we observe the same behavior as with $\alpha=0.65$ for
the zero-potential (constant-potential) case, indicating the universality of
the anomalous transport behavior.

##  Appendix B Convergence of numerical results

###  B.1 Discrete time step dependence

Throughout the main text, we selected a discrete time step of $\Delta
t=0.1J^{-1}$. The choice of $\Delta t$ is subtle because part of the system’s
evolution time is composed of waiting time, which directly depends on $\Delta
t$. To verify the convergence of our results with respect to $\Delta t$, we
selected different $\Delta t$ values of $0.2J^{-1}$, $0.1J^{-1}$, and
$0.05J^{-1}$, and compared their results. For the zero-potential shown in Fig.
7(a), our simulation results are consistent with the smaller $\Delta
t=0.05J^{-1}$, indicating that the $\Delta t$ chosen is small enough to
neglect numerical errors caused by the time discretization. For the constant-
potential shown in Fig. 7(b), the results differs for different $\Delta t$, we
infer that this is because the system continuously evolves due to the
existence of $V$, so the state is slightly distinct for different $\Delta t$
when the hopping process starts, and the distinction can be amplified during
time evolution. However, the non-convergence can be revised by the large
number of samples because the small difference between these $\Delta t$ values
almost does not affect the distribution of time $\psi(t)$. As shown in Fig.
7(c), the ensemble averaged square of width of wavepacket $\langle
W^{2}(t)\rangle_{ens}$ is converged to the smaller $\Delta t=0.05J^{-1}$,
indicating that the $\Delta t$ chosen in our simulation is sufficient.

###  B.2 System size dependence

The main text outlines a single-particle simulation on a 1D lattice with
finite system sizes of $L=8000$ for the zero-potential case and $L=2500$ for
the constant-potential case, which limits the maximum simulation time. After
this period, the wave packet will reflect off the boundaries of the $1D$
lattice under open boundary conditions. It is necessary to check the system
boundary dependence of our results. As shown in Fig. 8(a) and (b), the results
converge well for different system sizes, indicating we have ruled out
boundary effects.

## References

  * Montroll and Weiss (1965) E. W. Montroll and G. H. Weiss, Journal of Mathematical Physics 6, 167 (1965).
  * Weiss (1994) G. H. Weiss, Journal of Statistical Physics 79, 497 (1994).
  * Klafter and Sokolov (2011) J. Klafter and I. M. Sokolov, Oxford University Press (2011).
  * Scher and Montroll (1975) H. Scher and E. W. Montroll, Phys. Rev. B 12, 2455 (1975).
  * Bouchaud and Georges (1990) J.-P. Bouchaud and A. Georges, Physics Reports 195, 127 (1990).
  * Barkai  _et al._ (2000) E. Barkai, R. Metzler, and J. Klafter, Phys. Rev. E 61, 132 (2000).
  * Weiss (2013) M. Weiss, Phys. Rev. E 88, 010101 (2013).
  * Shlesinger (1974) M. F. Shlesinger, Journal of Statistical Physics 10, 421 (1974).
  * Pfister and Scher (1977) G. Pfister and H. Scher, Phys. Rev. B 15, 2062 (1977).
  * Geisel  _et al._ (1985) T. Geisel, J. Nierwetberg, and A. Zacherl, Phys. Rev. Lett. 54, 616 (1985).
  * Solomon  _et al._ (1993) T. H. Solomon, E. R. Weeks, and H. L. Swinney, Phys. Rev. Lett. 71, 3975 (1993).
  * Ohtsuki and Kawarabayashi (1997) T. Ohtsuki and T. Kawarabayashi, Journal of the Physical Society of Japan 66, 314 (1997).
  * Mandelbrot (1963) B. Mandelbrot, The Journal of Business 36, 394 (1963).
  * Masoliver  _et al._ (2003) J. Masoliver, M. Montero, and G. H. Weiss, Phys. Rev. E 67, 021112 (2003).
  * Scalas (2006) E. Scalas, Physica A 362, 225 (2006).
  * Helmstetter and Sornette (2002) A. Helmstetter and D. Sornette, Phys. Rev. E 66, 061104 (2002).
  * Mega  _et al._ (2003) M. S. Mega, P. Allegrini, P. Grigolini, V. Latora, L. Palatella, A. Rapisarda, and S. Vinciguerra, Phys. Rev. Lett. 90, 188501 (2003).
  * H$\ddot{\rm{o}}$fling and Franosch (2013) F. H$\ddot{\rm{o}}$fling and T. Franosch, Reports on Progress in Physics 76, 046602 (2013).
  * Metzler and Klafter (2000) R. Metzler and J. Klafter, Physics Reports 339, 1 (2000).
  * Sokolov  _et al._ (2002) I. Sokolov, J. Klafter, and A. Blumen, Physics Today - PHYS TODAY 55, 48 (2002).
  * Bouchaud (1992) J. P. Bouchaud, J. Phys. I France 2, 1705 (1992).
  * Barkai and Margolin (2004) E. Barkai and G. Margolin, Israel Journal of Chemistry 44, 353 (2004).
  * Margolin and Barkai (2005) G. Margolin and E. Barkai, Phys. Rev. Lett. 94, 080601 (2005).
  * Lubelski  _et al._ (2008) A. Lubelski, I. M. Sokolov, and J. Klafter, Phys. Rev. Lett. 100, 250602 (2008).
  * He  _et al._ (2008) Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev. Lett. 101, 058101 (2008).
  * Zotos  _et al._ (1997) X. Zotos, F. Naef, and P. Prelovsek, Phys. Rev. B 55, 11029 (1997).
  * Batchelor (2005) M. T. Batchelor, Journal of Physics A: Mathematical and General 38, 3245 (2005).
  * Miller  _et al._ (2020) S. Miller, J. Tennyson, T. R. Geballe, and T. Stallard, Rev. Mod. Phys. 92, 035003 (2020).
  * Karrasch  _et al._ (2014) C. Karrasch, J. E. Moore, and F. Heidrich-Meisner, Phys. Rev. B 89, 075139 (2014).
  * $\rm{\check{Z}nidari\check{c}}$  _et al._ (2016) M. $\rm{\check{Z}nidari\check{c}}$, A. Scardicchio, and V. K. Varma, Phys. Rev. Lett. 117, 040601 (2016).
  * Blake  _et al._ (2017) M. Blake, R. A. Davison, and S. Sachdev, Phys. Rev. D 96, 106008 (2017).
  * Friedman  _et al._ (2020) A. J. Friedman, S. Gopalakrishnan, and R. Vasseur, Phys. Rev. B 101, 180302 (2020).
  * Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Agarwal  _et al._ (2015) K. Agarwal, S. Gopalakrishnan, M. Knap, M. M$\rm{\ddot{u}}$ller, and E. Demler, Phys. Rev. Lett. 114, 160401 (2015).
  * Bar Lev  _et al._ (2015) Y. Bar Lev, G. Cohen, and D. R. Reichman, Phys. Rev. Lett. 114, 100601 (2015).
  * Luitz  _et al._ (2016) D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 93, 060201 (2016).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Chang  _et al._ (2018) D. E. Chang, J. S. Douglas, A. Gonz$\rm{\acute{a}}$lez-Tudela, C.-L. Hung, and H. J. Kimble, Rev. Mod. Phys. 90, 031002 (2018).
  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Luca D’Alessio and Rigol (2016) A. P. Luca D’Alessio, Yariv Kafri and M. Rigol, Advances in Physics 65, 239 (2016).
  * Struck  _et al._ (2012) J. Struck, C. $\rm{\ddot{O}}$lschl$\rm{\ddot{a}}$ger, M. Weinberg, P. Hauke, J. Simonet, A. Eckardt, M. Lewenstein, K. Sengstock, and P. Windpassinger, Phys. Rev. Lett. 108, 225304 (2012).
  * Cai  _et al._ (2017) Z. Cai, C. Hubig, and U. Schollwöck, Phys. Rev. B 96, 054303 (2017).
  * Harper  _et al._ (2020) F. Harper, R. Roy, M. S. Rudner, and S. L. Sondhi, Annual Review of Condensed Matter Physics 11, 345 (2020).
  * Guo  _et al._ (2023) C. Guo, W. Cui, and Z. Cai, Phys. Rev. A 107, 033330 (2023).
  * Sieberer  _et al._ (2013) L. M. Sieberer, S. D. Huber, E. Altman, and S. Diehl, Phys. Rev. Lett. 110, 195301 (2013).
  * Cai and Barthel (2013) Z. Cai and T. Barthel, Phys. Rev. Lett. 111, 150403 (2013).
  * Marino and Diehl (2016) J. Marino and S. Diehl, Phys. Rev. Lett. 116, 070407 (2016).
  * Lee  _et al._ (2022) K. H. Lee, V. Balachandran, C. Guo, and D. Poletti, Phys. Rev. E 105, 024120 (2022).
  * Erbanni  _et al._ (2023) R. Erbanni, X. Xu, T. F. Demarie, and D. Poletti, Phys. Rev. A 108, 032619 (2023).
  * (50) Y.-P. Wang, C. Fang, and J. Ren,  arXiv:2310.03069 .
  * $\rm{\check{Z}nidari\check{c}}$ (2024) M. $\rm{\check{Z}nidari\check{c}}$, Phys. Rev. B 109, 075105 (2024).
  * Zaburdaev  _et al._ (2015) V. Zaburdaev, S. Denisov, and J. Klafter, Rev. Mod. Phys. 87, 483 (2015).
  * Schuckert  _et al._ (2020) A. Schuckert, I. Lovas, and M. Knap, Phys. Rev. B 101, 020416 (2020).
  * Joshi  _et al._ (2022) M. K. Joshi, F. Kranzl, A. Schuckert, I. Lovas, C. Maier, R. Blatt, M. Knap, and C. F. Roos, Science 376, 720 (2022).
  * Nandy  _et al._ (2017) S. Nandy, A. Sen, and D. Sen, Phys. Rev. X 7, 031034 (2017).
  * Halimeh  _et al._ (2018) J. C. Halimeh, M. Punk, and F. Piazza, Phys. Rev. B 98, 045111 (2018).
  * Pan  _et al._ (2020) L. Pan, X. Chen, Y. Chen, and H. Zhai, Nature Physics 16, 767 (2020).
  * Kos  _et al._ (2021) P. Kos, B. Bertini, and T. c. v. Prosen, Phys. Rev. Lett. 126, 190601 (2021).
  * Cai (2022) Z. Cai, Phys. Rev. Lett. 128, 050601 (2022).
  * Klafter  _et al._ (1987) J. Klafter, A. Blumen, and M. F. Shlesinger, Phys. Rev. A 35, 3081 (1987).


In fact, (<ref>), combined with the method in [15] will be sufficient for validating (i)$\Leftrightarrow$(ii) in Rademacher type Theorem <ref>. Unfortunately, since we are coping with Hörmander vector field, a barrier arises when modifying their proofs. Indeed, their method uses a result by [17] that every $x$-measurable Hamiltonian $H$ can be approximated by a sequence of smooth Hamiltonians $\{H_n\}_n$ such that two intrinsic distances $\wz d_\lz^{H_n}$ and $d_\lz^{H_n}$ induced by $H_n$ by means of (<ref>) and (<ref>) satisfy $ \lim_{n \to \fz}d_\lz^{H_n} = d_\lz$ and $ \limsup_{n \to \fz}\wz d_\lz^{H_n} \le \wz d_\lz$ uniformly on $ U\times U$ respectively. The process of the proof of [17] is based on a $C^1$ Lusin approximation property for curves. Namely, given a Lipschitz curve $\gz: \ [0,1] \to U$ joining $x$ and $y$, for any $\ez >0$, there exists a $C^1$ curve $\widetilde{\gz}: \ [0,1] \to U$ with the same endpoints such that

$$ \mathcal{L}^1 ( \{t \in [0,1]  \ | \ \widetilde{\gz}(t) \neq \gz(t) \quad \text{or} \quad \widetilde{\gz}'(t)\neq \gz'(t)\})<\ez$$

where $\mathcal{L}^1$ denotes the one dimensional Lebesgue measure. Besides,

$$ \|\widetilde{\gz}\|_{L^\fz} \le c\|\gz\|_{L^\fz}, $$

for some constant $c$ depending only on $n$.
Although this version of $C^1$ Lusin approximation property holds for horizontal curves in Heisenberg groups ([47]) and step 2 Carnot groups ([43]),
it fails for some horizontal curve in Engel group ([47]).

In summary, it is difficult to generalize the properties of the pseudo metric $\wz d_\lz$ not only from Euclidean space to the case of Hörmander vector fields but also from lower-semicontinuous $H(x,p)$ to measurable $H(x,p)$. Hence we would like to pose the following open problem.

Under the assumptions (H0)-(H3), does (<ref>) holds?

thmapp[lemapp] Theorem
corapp[lemapp] Corollary
remapp[lemapp]  Remark
defnapp[lemapp]  Definition

§ RADEMACHER'S THEOREM  IN EUCLIDEAN DOMAINS—REVISIT

In this appendix, we state some consequence  of Rademacher's theorem (Theorem <ref>) for   Sobolev and Lipschitz classes,  see Lemma <ref> and Lemma <ref> below.
They were well-known in the literature and also partially  motivated our Theorem <ref> and Corollary <ref>.
For reader's convenience, we give the details.

Recall that $\Omega\subset\rn$ is  always a domain.
The homogeneous Sobolev space $\dot W^{1,\fz}(\Omega)$ is the collection of all functions $u\in L^1_\loc(\Omega)$ with its distributional derivative $\nabla u=(\frac{\partial u}{\partial x_i})_{1\le i\le n}\in L^\fz(\Omega)$.   We equip $\dot W^{1,\fz}(\Omega)$
with the semi-norm
$$\|u\|_{\dot W^{1,\fz}(\Omega)}=\|\nabla u\|_{L^\fz(\Omega)}.$$

$\dot W^{1,\fz}_\loc(\Omega)$ as the collection of all functions $u$ in $\Omega$ so that $u\in \dot W^{1,\fz}(V)$ whenever $V\Subset\Omega$. Here and below,  $V\Subset\Omega$ means that
$V$ is bounded domain with $\overline V\subset\Omega$.
On the other hand, denote by $\lip(\Omega)$
the collection of all Lipschitz functions $u$ in $\Omega$,
that is,   all functions $u$ satisfying (<ref>). We equip $\lip(\Omega)$ with the semi-norm

$$\lip(u,\Omega)=\sup_{x,y\in\Omega,x\ne y}\frac{|u(y)-u(x)|}{|x-y|}=\inf\{ \mbox{$\lz\ge 0$ satisfiying \eqref{lip}}\} .$$

Denote by $\lip_\loc(\Omega)$ the collection of all   functions $u$ in $\Omega$  so that $u\in\lip(V)$ for any $V\Subset\Omega$.
denote by $\lip^\ast(\Omega)$ the  collection of all functions  $u$ in $\Omega$

\begin{equation}\label{suplip}\lip^\ast(u,\Omega):=\sup_{x \in\Omega }\lip u(x)<\fz.
\end{equation}

Obviously, $ \lip(\Omega)\subset\lip^\ast(\Omega)$ with the seminorm bound
$\lip^\ast(u,\Omega)\le \lip(u,\Omega).$

Next, we have the following relation.

We have $  \lip^\ast(\Omega)\subset\lip_\loc(\Omega)$.
For any convex subdomain $V\subset\Omega$ and $u\in  \lip^\ast(\Omega)$, we have

\begin{equation}\label{ast}
|  u(x)-  u(y)|\le \|u\|_{\lip^\ast(V ) } |x-y|\quad\forall x,y\in V.
\end{equation}

Let $u\in  \lip^\ast(\Omega)$. To see $ u\in \lip_\loc(\Omega)$, it suffices to
prove $u\in\lip(B)$ for any ball $B\Subset\Omega$.  Given any $x,y\in B$, denote by $\gz(t)=x+t(y-x)$. Since $A_{x,y}=\sup_{t\in[0,1]} \lip u(\gz(t))<\fz$,
for each $t\in[0,1]$ we can find $r_t>0$ such that
$|u(\gz(s))-u(\gz(t))|\le A_{x,y}|\gz(s)-\gz(t)|=A_{x,y}|s-t||x-y| $ whenever $|s-t|\le r_t$ and $s \in [0,1]$.
Since $[0,1]\subset\cup_{t\in[0,1]}(t-r_t,t+r_t)$ we can find an increasing sequence
$t_i\in[0,1]$ with $t_0 = 0$ and $t_N=1$ such that $[0,1]\subset \cup_{i=1}^N(t_i-\frac12r_{t_i},t_i+\frac12 r_{t_i})$. Write $x_i=\gz(t_i)$ for $i =0 , \cdots , N$. We have

\sum_{i=0}^{N-1}|u(x_i)-u(x_{i+1})|\le A_{x,y}\sum_{i=0}^{N-1}|x_i-x_{i+1}|=A_{x,y}|x-y|.$$

Noticing that $A_{x,y}  \le \lip^\ast(u,B) \le \lip^\ast(u,\Omega) < \fz$ for all $x,y \in B$, we deduce that $u \in \lip(B)$ and hence $u \in \lip_\loc(\Omega)$.

If $V \Subset \Omega$ is convex, then for any $x,y \in \Omega$, the line-segment joining $x$ and $y$ lies in $V$. Hence similar to the above discussion, we have

$$ |u(x)-u(y)|\le A_{x,y}|x-y|, \ \forall x,y \in V $$

$A_{x,y}  \le \|u\|_{\lip^\ast(V)}< \fz$ for all $x,y \in V$. Therefore, (<ref>) holds and the proof is complete.

On the other hand, functions $\dot W^{1,\fz}_\loc(\Omega)$ admit  continuous representatives.

(i) Each $u\in  \dot W^{1,\fz}_\loc(\Omega)$
admits a unique continuous representative $\wz u$, that is,
$\wz u\in \dot W^{1,\fz}_\loc(\Omega)$    with $\wz u=u$ almost everywhere in $\Omega$. Moreover, $\wz u\in \lip_\loc(\Omega)$, and
for any convex subdomain $V\subset\Omega$, we have

$$|\wz u(x)-\wz u(y)|\le \|u\|_{\dot W^{1,\fz}(\Omega)}|x-y|\quad\forall x,y\in V.$$

(ii) Each   $u\in  \dot W^{1,\fz} (\Omega)$
admits a unique continuous representative $\wz u$, that is,  $\wz u\in  \dot W^{1,\fz} (\Omega)$   with $\wz u=u$ almost everywhere in $\Omega$. Moreover,
$\wz u\in \lip^\ast(\Omega)$ with $\lip^\ast(\wz u,\Omega)\le \|u\|_{\dot W^{1,\fz}(\Omega)}$.

Since (ii) can be shown in a similar way as (i), we only prove (i).
Given any  convex domain $V\Subset\Omega$, for any pair $x,y$ of  Lebesgue points  of $ u$, we have

u(y)-u(x)=\lim_{\dz\to0} [u\ast\eta_\dz(y)-u\ast\eta_\dz(x)]=\lim_{\dz\to 0}\nabla (u\ast\eta_\dz)(x+t_\dz(y-x))\cdot (y-x)$$

where $\eta_\dz$ is the standard mollifier in $\rr^n$ with its support  ${\rm spt}\eta_\dz \subset B(0,\dz)$ and $t_\dz \in [0,1]$.
Also, since for any $z\in V$,

$$|\nabla (u\ast\eta_\dz)(z)|= |(\nabla  u)\ast\eta_\dz)(z)|\le \|\nabla u\|_{L^\fz(B(z,\dz))}, $$

we deduce that for any pair $x,y$ of Lebesgue points of $ u$,

$$ |u(y)-u(x)|\le \|\nabla u\|_{L^\fz(V)} |y-x|.$$

If $z \in V$ is not a Lebesgue point of $u$, let $\{z_i\}_{i \in \nn} \subset V$ be a sequence of Lebesgue points of $u$ converging to $z$. We have

$$  \lim_{i \to \fz}  |u(z_i)-u(z_{i+1})|\le  \lim_{i \to \fz} \|\nabla u\|_{L^\fz(V)} |z_i-z_{i+1}| = 0,  $$

which implies $\{u(z_i)\}_{i \in \nn} \subset V$ is a Cauchy sequence.
Since $\|u\|_{L^\fz(V)} < \fz$, we know $\{u(z_i)\}_{i \in \nn}$ has a limit in $\rr$ independent of the choice of the sequence $\{u(z_i)\}_{i \in \nn}$. Define
$\wz u(z):= u(z) $ if $z \in V$ is a Lebesgue point of $u$ and $\wz u(z)= \lim_{i \to \fz} u(z_i)$ if $z \in V$ is not a Lebesgue point of $u$ where $\{z_i\}_{i \in \nn} \subset V$ is a sequence of Lebesgue points of $u$ converging to $z$.
We know $\wz u : V \to \rr$ is well-defined and moreover,

$$ | \wz u(y)-\wz u(x)|\le \|\nabla u\|_{L^\fz(V)} |y-x|\quad\forall x,y\in V.$$

Thus $\wz u\in\lip(V)$ with   $\sup_{x\in V}\lip  \wz u(x)\le \lip(\wz u,V)\le \|\nabla u\|_{L^\fz(V)}$. In particular, $\wz u$ is continuous, which shows (i).

Thanks to  lemma <ref>, below
for any function $u\in \dot W^{1,\fz}_\loc(\Omega)$ or $u\in \dot W^{1,\fz}(\Omega)$,
up to considering its continuous representative $\wz u$, we may assume that $u$ is continuous.
Under this assumption, Lemma <ref> further gives  $\dot W^{1,\fz}_\loc(\Omega)\subset \lip_\loc(\Omega)$,
and  $\dot W^{1,\fz} (\Omega)\subset \lip^\ast(\Omega)$ with a norm bound $\lip^\ast(u,\Omega)\le \|u\|_{\dot W^{1,\fz}(\Omega)}$.
Rademacher's theorem (Theorem <ref>) tells that their  converse are also true. Indeed, we have the following.

We have
$ \dot W^{1,\fz}_\loc(\Omega)=\lip_\loc(\Omega)$  and
$\lip( \Omega)\subset \dot W^{1,\fz} (\Omega)= \lip^\ast(\Omega)$ with
$\lip^\ast(u,\Omega)= \|u\|_{\dot W^{1,\fz}(\Omega)}\le \lip(u,\Omega)$.

(ii)  If $\Omega$ is convex, then  $\lip( \Omega)= \dot
W^{1,\fz} (\Omega)= \lip^\ast(\Omega)$ with
$\lip^\ast(u,\Omega)= \|u\|_{\dot W^{1,\fz}(\Omega)}= \lip(u,\Omega)$.

If $u\in \lip_\loc(\Omega)$, applying Rademacher's theorem (Theorem <ref>) to all subdomains  $V\Subset\Omega$, one has $u\in \dot W^{1,\fz}_\loc(\Omega)$
and $ |\nabla u(x)|=\lip u(x)$ for almost all $x\in \Omega$
(whenever $u$ is differentiable at $x$). Hence $\lip_\loc(\Omega) \subset W^{1,\fz}_\loc(\Omega)$. Combining Lemma <ref>(i), we know
$\lip_\loc(\Omega) = W^{1,\fz}_\loc(\Omega)$.

If $u\in \lip^\ast(\Omega)$,   that is, $\lip^\ast(u,\Omega) = \sup_{x\in\Omega}\lip u(x)<\fz$.
We have $u\in \lip_{\loc}(\Omega)$ and hence $u\in \dot W^{1,\fz}_\loc(\Omega)$
and   $ |\nabla u(x)|= \lip u(x) \le  \lip^\ast(u,\Omega)<\fz$ for almost all $x\in \Omega$.
Thus $u\in \dot W^{1,\fz}(\Omega)$.

By definition, it is obvious that $\lip( \Omega)\subset  \lip^\ast(\Omega)$. Hence Lemma <ref> (i) holds.

(ii) By Lemma <ref> (i), we only need to show $\lip^\ast(\Omega) \subset \lip( \Omega)$. Applying Lemma <ref> with $V=\Omega$ therein, (<ref>) becomes

$$   \frac{ |  u(x)-  u(y)|}{|x-y|} \le \|u\|_{\lip^\ast(\Omega ) } \quad\forall x,y\in \Omega.  $$

Taking supremum among all $x,y \in \Omega$ in the left hand side of the above inequality, we arrive at

$$   \|u\|_{\lip(\Omega ) }  \le \|u\|_{\lip^\ast(\Omega ) },  $$

which gives the desired result.

(i) Lemma <ref> and Lemma <ref> fail if we   relax $ \sup_{x\in\Omega}\lip u(x)$ in the definition  (<ref>)  to be $\|\lip u\|_{L^\fz(\Omega)}=\esup_{x\in\Omega}\lip u(x)$.
This is witted by  the standard Cantor function $w$ in $[0,1]$.
Denote by  $E$ the standard Cantor set.
It is well-known that $w$ is continuous but not absolute continuous in $[0,1]$. Since  Lipschitz functions are always absolutely continuous, we know that
$w$ is neither Lipschitz nor locally Lipschitz in $\Omega=(0,1)$,
and hence $w\notin \lip_\loc(\Omega)$.
On the other hand,  observe that $\Omega\setminus E$ consists of a sequence of open intervals which mutually disjoint, and $w$ is a constant in each such intervals and hence
$\lip w(x)=0 $  therein. So we know that $\lip w(x)=0 $
in $\Omega\setminus E$.
Since $|E|=0$, we have $\|\lip u\|_{L^\fz(\Omega)}=0$.

(ii) In general,  if $\Omega$ is not convex, one cannot expect that  $ \dot W^{1,\fz} (\Omega)\subset \lip (\Omega)$ with a norm bound.  Indeed, consider the planar domain

\begin{equation}\label{domainu}U:=\{x=(x_1,x_2) \in \rr^2 \ | \ |x|<1 \} \setminus [0,1) \times \{0\}.
\end{equation}

in the polar coordinate $(r,\tz)$, let $w: U \to \rr$ be

$$w(r,\tz) := r\tz \mbox{ for all $0< r <1$ and $0< \tz< 2\pi$}.$$

One can show that $w\in \dot W^{1,\fz}(U)$ so that $w(x_1,x_2)< \pi/3$ when $1/2\le x_1<1$ and $0<x_2<1/10$, and $w(x_1,x_2)> 5\pi/6$
when $1/2\le x_1<1$ and $-1/10<x_2<0$.  One has  $\lip (w,\Omega)=\fz$ and hence $w\notin \lip(\Omega)$.

The example in Remark <ref> (ii) also indicates that
the Euclidean distance does not match the geometry of domains and hence
$\lip (\Omega)$ defined via Euclidean distance is not the prefect one to understand  $ \dot W^{1,\fz} (\Omega)$.

Instead of Euclidean distance,  for any domain $\Omega$, we consider the intrinsic distance

\begin{equation}\label{dual}
d^\Omega_E(x,y)=\inf\{\ell(\gz) \ | \ \mbox{$\gz:[0,1]\to\Omega$ is absolute coninuous curve joining $x,y$}\},
\end{equation}

where $\ell(\gz):=\int_0^1|\dot\gz(t)|\,dt$ is the Euclidean length.
We have the dual formula.

(i) For any $x,y\in\Omega$,

\begin{equation}\label{intrinsic}
=\sup\{u(y)-u(x) \ | \ u\in\dot W^{1,\fz}(\Omega),\ \|\nabla u\|_{L^\fz(\Omega)}\le 1\}.
\end{equation}

(ii) If  $x,y\in\Omega$ with $|x-y|\le \dist(x,\partial\Omega)$, then  $d^\Omega_E(x,y)=|x-y|$.

(iii)  If $\Omega$ is convex, then  $d^\Omega_E(x,y)=|x-y|$ for all $x,y\in\Omega$.

(i) Set

\begin{equation}\label{wzd}
\wz d^\Omega_E(x,y)
=\sup\{u(y)-u(x) \ | \ u\in\dot W^{1,\fz}(\Omega),\ \|\nabla u\|_{L^\fz(\Omega)}\le 1\}.
\end{equation}

Notice that $d^\Omega_E(x, \cdot) \in \Lip^\ast(\Omega) = \dot W^{1,\fz}(\Omega)$ (Lemma <ref> (i)) and $\|\nabla d^\Omega_E(x, \cdot)\|_{L^\fz(\Omega)}\le 1$ for all $x \in \Omega$. Hence letting $d^\Omega_E(x, \cdot)$ be the test function in (<ref>), we see

$$  d^\Omega_E(x,y) \le \wz d^\Omega_E(x,y) \ \forall x,y \in \Omega.$$

To see the contrary, fix $x,y \in \Omega$. Let $\{u_i\}_{i \in \nn}$ be a sequence of test functions in (<ref>) such that

$$  \wz d^\Omega_E(x,y)
= \lim_{i \to \fz } (u_i(y)-u_i(x)).   $$

Let $\gz:[0,1] \to \Omega$ be an arbitrary absolute continuous curve joining $x$ and $y$. Then there exists a domain $U \Subset \Omega$ with $\gz \subset U$. Let $\{\eta_\dz\}_{\dz >0}$ be the standard mollifiers in $\rr^n$. For each $i \in \nn$, we know $u_i \ast \eta_\dz \in C^\fz(\Omega)$ and $\|\nabla (u_i \ast \eta_\dz)\|_{L^\fz(U)}\le \|\nabla u_i \|_{L^\fz(U)}  \le 1$.
Then we have

\begin{align*}
\wz d^\Omega_E(x,y) & = \lim_{i \to \fz} [u_i(y)-u_i(x)]  \\
& =\lim_{i \to \fz} \lim_{\dz \to 0}[u_i \ast \eta_\dz(y) - u_i \ast \eta_\dz(x)]  \\
& =\lim_{i \to \fz} \lim_{\dz \to 0} \int_{0}^1 \nabla (u_i \ast \eta_\dz) (\gz(t)) \cdot \dot \gz (t) \, dt \\
& \le \lim_{i \to \fz} \lim_{\dz \to 0} \int_{0}^1 |\nabla (u_i \ast \eta_\dz)(\gz(t))|  |\dot \gz (t)| \, dt \\
& \le  \int_{0}^1  |\dot \gz (t)| \, dt \\
& = \ell(\gz).
\end{align*}

Finally, taking infimum among all absolute continuous curves joining $x$ and $y$ in the above inequality, we conclude

$$  \wz d^\Omega_E(x,y) \le d^\Omega_E(x,y) \ \forall x,y \in \Omega.$$

(ii) If $|x-y|\le \dist(x,\partial\Omega)$, then the line-segment $\gz$ joining $x$ and $y$ is contained in $\Omega$. Letting $\gz$ be the absolute continuous curve in (<ref>),

$$  |x-y|  \le d^\Omega_E(x,y) \le  \ell(\gz) = |x-y| \ \forall x,y \in \Omega.$$

(iii) If $\Omega$ is convex, for all $x,y\in\Omega$, since the line-segment joining them is contained in $\Omega$,  similarly to (ii),  we have $d^\Omega_E(x,y)=|x-y|$.
The proof is complete.

Note that if $\Omega$ is not convex, one cannot expect   $d^\Omega_E(x,y)=|x-y|$ for all $x,y\in\Omega$. Indeed, if $\Omega$ is given by the domain $U$ as in (<ref>),
for points $(1/2,\ez)$ and $(1/2,-\ez)$ with $\ez\in(0,1/10)$, the Euclidean distance between them is $2\ez$.
However, note that any curve $\gz:[0,1]\to\Omega$ joining them must have intersection with $(-1,0)\times\{0\}$, which is call $z$. One then deduce
$$\ell(\gz)\ge |(1/2,\ez)-z|+|(-1/2,\ez)-z|\ge 1/2+1/2=1+2\ez.$$

Thus the intrinsic distance between $(1/2,\ez)$ and $(1/2,-\ez)$ is always larger than or
equals to $1+2\ez$.

With in Lemma <ref>, we show that the Lipschitz spaces defined via the intrinsic distance perfectly match
with the Sobolev space $ \dot W^{1,\fz}(\Omega)$, see Lemma <ref> below.
Denote by $\lip_{d^\Omega_E}(\Omega) $  the collection of all Lipschitz functions $u$ in $\Omega$ with respect to $d^\Omega_E$, that is,

$$ \lip_{d^\Omega_E}(u,\Omega):=\sup\frac{|u(x)-u(y)|}{d^\Omega_E(x,y)}<\fz.$$

We also denote by $\lip^\ast_{d^\Omega_E}(\Omega)$ the  collection of all functions  $u$ in $\Omega$

$$\lip^\ast_{d^\Omega_E}(u,\Omega):=\sup_{x \in\Omega }\lip_{d^\Omega_E} u(x)<\fz.$$

We have    $ \lip_{d^\Omega_E}(\Omega)= \dot W^{1,\fz} (\Omega)= \lip^\ast(\Omega)$ and

$$\|\nabla u\|_{L^\fz(\Omega)}=\lip_{d^\Omega_E} (u,\Omega)=\lip^\ast(u,\Omega)=\sup_{x\in\Omega}\lip_{d^\Omega_E} u(x).$$

Recall that Lemma <ref> gives $d^\Omega_E(x,y)=|x-y|$ whenever $|x-y|\le \dist(x,\partial\Omega)$.
One then has   $\lip_{d^\Omega_E}(\Omega) \subset\lip_\loc(\Omega)$,
and moreover,
$\lip u(x)= \lip_{d^\Omega_E} u(x)$ for all $x\in\Omega$, which gives $\lip^\ast_{d^\Omega_E}(\Omega)=\lip^\ast (\Omega)$.

Next, we show $\dot W^{1,\fz} (\Omega) \subset \lip_{d^\Omega_E}(\Omega)$ and $ \lip_{d^\Omega_E} (u,\Omega) \le \|\nabla u\|_{L^\fz(\Omega)}$. Let $u \in \dot W^{1,\fz} (\Omega)$.  Then $\|\nabla u\|_{L^\fz(\Omega)} =:\lz < \fz.$ If $ \lz >0$, then $\lz^{-1}u \in \dot W^{1,\fz} (\Omega)$ and $\|\nabla (\lz^{-1} u )\|_{L^\fz(\Omega)} =1$. Hence $\lz^{-1}u$ could be the test function in (<ref>), which implies

$$    \lz^{-1} u(y)- \lz^{-1}u(x) \le  d^\Omega_E (x,y) \ \forall x,y \in \Omega,  $$

or equivalently,

$$  \frac{|u(y)- u(x) |}{\|\nabla u\|_{L^\fz(\Omega)}} \le  d^\Omega_E (x,y) \ \forall x,y \in \Omega.  $$

Therefore, $u \in \lip_{d^\Omega_E}(\Omega)$ and $ \lip_{d^\Omega_E} (u,\Omega) \le \|\nabla u\|_{L^\fz(\Omega)}$. If $\lz=0$, then similar as the above discussion, we have for any $\lz'>0$

$$  \frac{|u(y)- u(x) |}{\lz'} \le  d^\Omega_E (x,y) \ \forall x,y \in \Omega.  $$

Therefore, $u \in \lip_{d^\Omega_E}(\Omega)$ and $ \lip_{d^\Omega_E} (u,\Omega) \le \lz'$ for any $\lz'>0$. Hence $ \lip_{d^\Omega_E} (u,\Omega) =0 = \|\nabla u\|_{L^\fz(\Omega)}$.

Moreover, we show $ \lip_{d^\Omega_E}(\Omega) \subset \lip^\ast(\Omega)$ and $ \lip^\ast(u,\Omega) \le \lip_{d^\Omega_E} (u,\Omega)$. Let $u \in \lip_{d^\Omega_E}(\Omega)$. Then $\lip_{d^\Omega_E} (u,\Omega)<\fz$. Since $\lip^\ast_{d^\Omega_E}(u,\Omega) \le \lip_{d^\Omega_E} (u,\Omega)$ and $\lip^\ast_{d^\Omega_E}(u,\Omega) = \lip (u,\Omega)$, we arrive at

$$   \lip (u,\Omega) \le \lip_{d^\Omega_E} (u,\Omega) < \fz. $$

Therefore, $u \in \lip^\ast(\Omega)$.

Finally, recalling that $\dot W^{1,\fz} (\Omega) = \lip^\ast(\Omega)$ and $\|\nabla u\|_{L^\fz(\Omega)} = \lip (u,\Omega) $ in Lemma <ref>, we finish the proof.

[1]
G. Aronsson,
Minimization problems for the functional $\sup_x F(x, f(x), f' (x))$.
Ark. Mat. 6 (1965), 33-53.

[2]
G. Aronsson, Minimization problems for the functional $\sup_x F(x, f(x), f' (x))$. II.
Ark. Mat. 6 (1966), 409-431.

[3]
G.  Aronsson, Extension of functions satisfying Lipschitz conditions.
Ark. Mat. 6 (1967), 551-561.

[4]
G. Aronsson, On the partial differential equation $u^2_x u_{xx} + 2u_x u_y u_{xy} + u^2_y u_{yy} = 0.$
Ark. Mat. 7, (1968), 395-425.

[5]
G. Aronsson,  Minimization problems for the functional $\sup_x F(x, f(x), f' (x))$. III.
Ark. Mat. 7 (1969),

[6]
S. N. Armstrong, M. G. Crandall, V. Julin, and C. K. Smart, Convexity criteria and uniqueness of absolutely minimizing functions.
Arch. Ration. Mech. Anal. 200 (2011), no. 2, 405-443.

[7]
G. Aronsson, M. G. Crandall, and P. Juutinen,   A tour of the theory of
absolutely minimizing functions.
Bull. Amer. Math. Soc. (N.S.) 41 (2004), 439-505.

[8]
N. Barron, Viscosity solutions and analysis in $L^\fz$. In: Nonlinear Analysis, Differential
Equations and Control (Montreal, QC, 1998). NATO Sci. Ser. C Math. Phys. Sci. 528.
Dordrecht: Kluwer Acad. Publ., pp. 1-60.(1999).

[10]  E. N. Barron, R. R. Jensen, and C. Y. Wang,
The Euler equation and absolute minimizers of $L^\fz$  functionals. Arch. Ration. Mech. Anal. 157 (2001),  255-283.

[11]
T. Bieske, On $\infty$-harmonic functions on the Heisenberg group. Comm. Partial Differential Equations 27 (2002), no. 3-4, 727-761.

[12]
A. Boutet de Monvel, D. Lenz and P. Stollmann,
Schnol's theorem for strongly local forms, Israel J. Math. 173 (2009), 189-211.

[13]
T. Champion,  L. De Pascale, and F.  Prinari, $\Gamma$-Convergence and absolute minimizers for supremal functionals. ESAIM Control Optim. Calc. Var. 10 (2004), 14-27.

[14]
V. M. Chernikov, S. K. Vodop'yanov, Sobolev Spaces and hypoelliptic
equations I,II. Siberian Advances in Mathematics. 6 (1996) no. 3, 27-67, no. 4, 64-96.
Translation from: Trudy In-ta matematiki RAN. Sib. otd-nie. 29 (1995), 7-62.

[15]
T. Champion  and L. De Pascale,
Principles of comparison with distance functions for absolute minimizers.
J. Convex Anal. 14 (2007),   515-541.

[16]
M. Crandall, An efficient derivation of the Aronsson equation. Arch. Ration. Mech. Anal. 167 (2003), 271-279.

[17]
A. Davini, Smooth approximation of weak Finsler metrics. Differential Integral Equations 18 (5) (2005) 509-530.

[18]
F. Dragoni, J. J. Manfredi and D. Vittone, Weak Fubini property and infinity harmonic functions in Riemannian and sub-Riemannian manifolds, Trans. Amer. Math.
Soc. Volume 365, Number 2, (2013), 837-859.

[21]
K. O. Friedrichs, The identity of weak and strong extensions of differential
operators. Trans. Amer. Math. Soc. 55 (1944), 132-151.

[22]
B. Franchi, P. Hajłasz, P. Koskela. Definitions of Sobolev classes on metric spaces. Annales de l'Institut Fourier, Volume 49 (1999) no. 6, pp. 1903-1924.

[23]
B. Franchi, R. Serapioni, F. Serra Cassano,
Meyers-Serrin type theorems and relaxation of variational integrals depending on vector fields. Houston J. Math. 22 (1996), no. 4, 859-890.

[24]
B. Franchi, R. Serapioni, F. Serra Cassano,
Approximation and imbedding
theorems for weighted Sobolev spaces associated with Lipschitz continuous vector
fields. Boll. Un. Mat. Ital. (7) 11-B (1997), 83-117.

[25]
R. L. Frank, D. Lenz, D. Wingert,
Intrinsic metrics for non-local symmetric Dirichlet forms and applications to spectral theory
Journal of Functional Analysis,
Volume 266, Issue 8,
Pages 4765-4808.

[26]
N. Garofalo and D. M. Nhieu, Isoperimetric and Sobolev inequalities for Carnot-Carathéodory
spaces and the existence of minimal surfaces, Comm. Pure Appl. Math. 49 (1996), 1081-1144.

[27]
N. Garofalo, D. Nhieu,
Lipschitz continuity, global smooth approximations and extension theorems for Sobolev functions in Carnot-Carathéodory spaces. J. Anal. Math. 74 (1998), 67-97.

[28]
R. Gariepy,  C. Y. Wang, and Y. Yu,
Generalized cone comparison principle for viscosity solutions of the Aronsson equation and absolute minimizers.
Comm. Partial Differential Equations 31 (2006),  1027-1046.

[29]
C. Y. Guo, C. Xiang, D. Yang,  $L^\infty$-variational problems associated to measurable Finsler structures.  Nonlinear Analysis. 132 (2015), 126-140.

[30]
P. Hajlasz, P. Koskela, Sobolev met Poincaré. Mem. Amer. Math. Soc. 145 (2000), no. 688.

[31]
J. Heinonen, P. Koskela, N. Shanmugalingam and J. T. Tyson, Sobolev spaces on metric measure spaces: an approach based on upper gradients, Cambridge Studies in Advanced Mathematics Series, Cambridge University Press, 2015.

[32]
R. Jensen,   Uniqueness of Lipschitz extensions: minimizing the sup norm of the gradient. Arch. Ration. Mech. Anal. 123 (1993), 51-74.

[33]
R. Jensen, C. Y. Wang, and Y.  Yu, Uniqueness and nonuniqueness of viscosity solutions to Aronsson's equation. Arch. Ration. Mech. Anal. 190 (2008), no. 2, 347-370.

[34]
D. Jerison, The Poincaré inequality for vector fields satisfying Hörmander's condition, Duke
Math. J. 53 (1986), 503-523.

[35]
D. Jerison, A. Sanchez-Calle, Subelliptic, second order differential operators.
In: Complex analysis, III (College Park, Md., 1985-86). pp. 46-77, Lecture Notes in
Math., 1277, Springer, 1987.

[37]
P. Juutinen, Minimization problems for Lipschitz functions via viscosity solutions, Ann.
Acad. Sci. Fenn. Math. Diss. No. 115 (1998), 53 pp.

[38]
P. Juutinen, Absolutely minimizing Lipschitz extensions on a metric space, Ann. Acad. Sci.
Fenn. Math. 27 (2002), no. 1, 57-67.

[39] P. Juutinen and N. Shanmugalingam,
Equivalence of AMLE, strong AMLE, and comparison with cones in metric measure space.   Math. Nachr. 279, (2006), 1083-1098.

[40] P. Koskela, N. Shanmugalingam, and Y. Zhou,
$L^\infty$-Variational problem associated to Dirichlet forms.  Math. Res. Lett.  19 (2012), 1263-1275.

[41]  P. Koskela, Y. Zhou,
Geometry and analysis of Dirichlet forms. Adv. Math. 231, (2012), 2755-2801.

[42]
P. Koskela, N. Shanmugalingam, and Y. Zhou,
Intrinsic geometry and analysis of diffusion process and $L^\infty$-variational problem.
Arch. Rational Mech. Anal.   214 (2014), no.1,  99-142.

[43]
E. Le Donne and G. Speight Lusin Approximation for Horizontal Curves in Step 2 Carnot Groups. Calculus of Variations and Partial Differential Equations 55(5) (2016), 1-22.

[44]
R. Monti , F S .Cassano, Surface measures in Carnot- Caratheodory spaces. Calc.var.partial Differential Equations. 13 (2001), 339-376.

[45]
A. Nagel, E. M. Stein and S, Wainger, Balls andmetrics defined by vectorfields L basic properties,
Acta Math. 155 (1985), 103-147.

[46]
P. Pansu, Metriques de Carnot-Carathéodory et quasiisométries des espaces symétriques
de rang un, Annals of Mathematics 129 (1989), 1-60.

[47]
G. Speight Lusin Approximation and Horizontal Curves in Carnot Groups. Revista Matematica Iberoamericana 32 (4) (2016), 1425-1446.

[48]
P. Stollmann,
A dual characterization of length spaces with application to Dirichlet metric spaces. Stud. Math. 198, (2010), 221-233.

[49]
K. T. Sturm,
Analysis on local Dirichlet spaces. I. Recurrence, conservativeness and $L^p$-Liouville properties. J. Rein. Angew. Math. 456, (1994), 173-196 .

[50]
K.T. Sturm, Is a diffusion process determined by its intrinsic metric? Chaos Solitons
Fractals 8 (1997), 1855-1860.

[51]
C. Y.  Wang,
The Aronsson equation for absolute minimizers of $L^\infty$ functionals associated with vector fields satisfying Hörmander's conditions. Trans. AMS. 359 (2007), 91-113.

Jiayin Liu

School of Mathematical Science, Beihang University, Changping District Shahe Higher Education Park
South Third Street No. 9, Beijing 102206, P. R. China

Department of Mathematics and Statistics, University of Jyväskylä, P.O. Box 35
(MaD), FI-40014, Jyväskylä, Finland

E-mail :

Yuan Zhou

School of Mathematical Sciences, Beijing Normal University, Haidian District Xinjiekou Waidajie No.19, Beijing 100875, P. R. China

E-mail :