††thanks: Corresponding author

# Spectral response of disorder-free localized lattice gauge theories

Nilotpal Chakraborty<EMAIL_ADDRESS>Max-Planck-Institut für Physik
komplexer Systeme, Nöthnitzer Straße 38, Dresden 01187, Germany  Markus Heyl
Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Straße 38,
Dresden 01187, Germany  Theoretical Physics III, Center for Electronic
Correlations and Magnetism, Institute of Physics, University of Augsburg,
D-86135 Augsburg, Germany  Petr Karpov  Max-Planck-Institut für Physik
komplexer Systeme, Nöthnitzer Straße 38, Dresden 01187, Germany  Roderich
Moessner  Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Straße
38, Dresden 01187, Germany

###### Abstract

We show that certain lattice gauge theories exhibiting disorder-free
localization have a characteristic response in spatially averaged spectral
functions: a few sharp peaks combined with vanishing response in the zero
frequency limit. This reflects the discrete spectra of small clusters of
kinetically active regions formed in such gauge theories when they fragment
into spatially finite clusters in the localized phase due to the presence of
static charges. We obtain the transverse component of the dynamic structure
factor, which is probed by neutron scattering experiments, deep in this phase
from a combination of analytical estimates and a numerical cluster expansion.
We also show that local spectral functions of large finite clusters host
discrete peaks whose positions agree with our analytical estimates. Further,
information spreading, diagnosed by an unequal time commutator, halts due to
real space fragmentation. Our results can be used to distinguish the disorder-
free localized phase from conventional paramagnetic counterparts in those
frustrated magnets which might realize such an emergent gauge theory.

Lattice gauge theories emerge as low-energy effective theories in a large
class of models relevant to quantum many-body physics, especially in
frustrated magnets Fradkin (2013); Moessner and Moore (2021). For example, in
quantum spin ice and hardcore dimer models, local constraints (such as the ice
rule of two spins pointing in and two pointing out locally) yield a 3+1
dimensional compact U(1) lattice gauge theory Moessner and Sondhi (2003);
Hermele  _et al._ (2004). The resulting algebraic spin correlations are,
strictly speaking, only present at zero temperature, in the absence of charges
(i.e., violations of the constraints). However, perhaps somewhat
counterintuitively, one can also find new phases by creating a high density of
such charges, provided that these are static. The concomitant sectors of the
lattice gauge theories have been shown to exhibit disorder-free localization
in one Smith  _et al._ (2017); Brenes  _et al._ (2018) and two dimensions
Karpov  _et al._ (2021); Chakraborty  _et al._ (2022). It thus becomes
necessary to distinguish disorder-free localization from a paramagnet, the
`standard' disordered phase, leading to the twin questions of to what extent
one can access novel phases at a finite density of static charges, and above
all, what are the distinguishing characteristics of the localized phases?

Figure 1: Accessing novel phases by increasing charge density (colored
circles: red, purple and black for $q_{r}=-1,1$ and $-2$,the arrows
$\rightarrow$ and $\downarrow$ represent spin $\ket{\uparrow}$, whereas
$\leftarrow$ and $\uparrow$ represent $\ket{\downarrow}$). In the presence of
a sufficient density of static charges, one can access the disorder-free
localized (DFL) phase ($\rm{C.D}>C_{\rm{dcr}}$). Rightmost panel: Sharp
spectral response in spatially averaged spectral functions (in this case, the
transverse component of the dynamical structure factor for $|\lambda/J|=0.2$)
with $S^{zz}(q,\omega)\rightarrow 0$ as $\omega\rightarrow 0$, in the
disorder-free localized phase.

Here, we present a spectral response dominated by discrete sharp peaks,
combined with vanishing weight in the zero frequency limit, as a
characteristic signature of the disorder-free localized phase. We show that
such a response arises due to formation of local motifs with discrete spectra
in such lattice gauge theories, a feature absent in conventional paramagnets
which behave diffusively at high temperatures Conlon and Chalker (2009). We
numerically obtain the spectral functions by calculating quantum dynamics of
the 2+1D U(1) quantum link model (QLM) via a controlled cluster expansion
method which allows us to access the dynamics deep in the localized phase. Our
conclusions are applicable to a large class of lattice gauge theories with
discrete degrees of freedom and static charges, of which the U(1) QLM is a
prototypical example.

Besides the aforementioned result we also gain the following insights: i)
Disorder-free localized lattice gauge theories provide an example of
information freezing due to real-space fragmentation in 2D. ii) In different
regimes of the U(1) QLM, there is a difference between local and global
response of spatially averaged spectral functions, such a difference could be
measured in experiments via local spectroscopy and serves as additional
evidence for disorder-free localization physics. There has also been a recent
surge in proposals for disorder-free localization Papaefstathiou  _et al._
(2020); Li  _et al._ (2021); Zhu and Heyl (2021); Hodson  _et al._ (2021), as
well as in proposals for realizing lattice gauge theories in various quantum
simulators Schweizer  _et al._ (2019); Yang  _et al._ (2020); Zhou  _et al._
(2022); Semeghini  _et al._ (2021) and stabilizing disorder-free localization
on such platforms Halimeh  _et al._ (2021). Hence, distinct experimental
signatures of disorder-free localization are of great current relevance.

Quantum Link model and percolation: The quantum link model (QLM) is a discrete
version of Wilsonian lattice gauge theories also occurring in gauge magnets
Chandrasekharan and Wiese (1997); Horn (1981); Orland and Rohrlich (1990);
Wiese (2013). Several models in condensed matter physics, for example the
toric code Kitaev (2003), quantum dimer and ice models Rokhsar and Kivelson
(1988); Moessner and Sondhi (2001); Shannon  _et al._ (2004) can be cast as
quantum link models. We use the 2+1 D U(1) QLM with

$\hat{H}=J\sum_{\square}(\hat{U}_{\square}+\hat{U}_{\square}^{\dagger})-\lambda\sum_{\square}(\hat{U}_{\square}+\hat{U}_{\square}^{\dagger})^{2},$
(1)

where
$\hat{U}_{\square}=\hat{S}^{+}_{\textbf{r},i}\hat{S}^{+}_{\textbf{r}+\textbf{i},j}\hat{S}^{-}_{\textbf{r}+\textbf{j},i}\hat{S}^{-}_{\textbf{r},j}$
flips all spins on a plaquette if they satisfy a certain orientation
(clockwise/anticlockwise in Fig. 1). Here, $\hat{S}^{+/-}_{\textbf{r},\mu}$
denotes the spin raising/lowering operators for the spin joining sites r and
$\textbf{r}+\hat{\mu}$, where $\hat{\mu}=\hat{i}/\hat{j}$ is a lattice vector.
The first term in Eq. (1) is a kinetic energy term, whereas the second is a
potential energy term Rokhsar and Kivelson (1988). Spin configurations are
characterized by the set of static charges on all sites of the lattice and the
charge on each site $q_{\textbf{r}}$ is determined by the eigenvalue equation
of the generator of gauge transformations (which evaluates the difference in
incoming and outgoing arrows from a site, see Fig. 1),
$G_{r}\ket{\psi}=q_{r}\ket{\psi}$ Brenes  _et al._ (2018). For a spin-$1/2$
Hamiltonian $q_{\textbf{r}}\in\\{-2,-1,0,1,2\\}$. This model has been shown to
host disorder-free localized physics and a continuous localization-
delocalization transition Karpov  _et al._ (2021); Chakraborty  _et al._
(2022). Crucially, the charges are static in our model.

We consider the problem with Eq. (1) as the effective Hamiltonian. Without any
static background charges we get square ice, i.e. the local constraint of two
in and two out spins on every vertex. On increasing charge density beyond a
critical density $C_{\rm{dcr}}$, there is a disorder-free localized phase
Karpov  _et al._ (2021); Chakraborty  _et al._ (2022). If one thinks of these
charges as being thermally activated, the disorder-free localized phase would
correspond to a high temperature phase (recently proposed in a different
context for 1D Halimeh  _et al._ (2022)). Temperature-dependent dynamics for
gauge theories has also been studied in context of the Kitaev model
Metavitsiadis  _et al._ (2017).
Our disorder-free localization can be studied using a percolation model, where
one starts from a classical random spin configuration, which fixes the initial
charge configuration, and then tracks all possible spin configurations induced
by single plaquette flips (1st term in Eq. (1)). Plaquette dynamics in such a
model is constrained by the static charge environment and certain charge
configurations will render a subset of plaquettes permanently frozen. If the
initial state hosts a sufficient density of background charges to fragment the
lattice into disconnected finite clusters of kinetically active plaquettes,
localization ensues (see Karpov  _et al._ (2021); Chakraborty  _et al._ (2022)
for percolation procedure).

Effective exact diagonalization using cluster expansion in the localized
phase: We employ an effective exact diagonalization procedure to calculate the
quantum dynamics of our 2D interacting system deep in the localized phase. The
percolation model gives us the real-space cluster structure, which in the
localized phase comprises only finite disconnected clusters. In contrast, in
the delocalized phase there is an infinite percolating cluster coexisting with
the finite clusters.
The decomposition of the whole lattice into only finite clusters in the
localized phase makes the quantum problem numerically tractable Chakraborty
_et al._ (2022). To calculate spatially averaged expectation values of time-
dependent correlation functions for the lattice, we thus use a cluster
expansion approach. We calculate the quantum dynamics individually for each
cluster using exact diagonalization. We then sum the results for all clusters
upto a certain size after which we find that including larger clusters adds
negligible weight to our observables. Plaquettes in a cluster are connected
via subsequent single plaquette flips, whereas plaquettes belonging to
different clusters (disjoint in real space) are uncorrelated. Since we
calculate dynamical expectation values of two point correlators, the only non-
zero contributions come from pairs of spins belonging to the same cluster.
Note that our cutoff procedure works only deep in the localized phase where
the mean cluster size of the corresponding non-percolating phase is not too
large. Closer to the transition, the mean cluster size increases and the
contribution from finite but large clusters becomes important, hence the
problem goes beyond the realm of our exact diagonalization numerics.
Frozen spins also contribute a delta function response exactly at zero
frequency. However, since our results pertain to non-zero frequencies, frozen
spins are only included in the denominator in Eqs. (2)(1st equation) and (3)
for calculating global spectral functions.
Besides the global response, we also calculate local spectral functions,
spatially averaged over a finite size cluster. To do so we average over all
spins belonging to all clusters of a particular size, over multiple disorder
realizations.

Information freezing due to fragmentation: We use our method to first
calculate unequal time commutators which are related to the susceptibility via
the Kubo formula. Such commutators (also their relatives, the out-of-time-
ordered correlators) quantify information spreading, a quantity related to
entanglement growth, and have been a measure of chaos and Hilbert-space
fragmentation in quantum many-body systems Maldacena  _et al._ (2016); Sala
_et al._ (2020); Khemani  _et al._ (2020); von Keyserlingk  _et al._ (2018).
We calculate

$\begin{split}&C_{G}(d,t)=\frac{1}{N}\sum_{\begin{subarray}{c}\textbf{r}_{i},\textbf{r}_{j}\in\mathcal{L}\\\
|\textbf{r}_{i}-\textbf{r}_{j}|=d\end{subarray}}\bra{\psi}\big{[}\hat{S^{z}_{\textbf{r}_{i}}}(t),\hat{S^{z}}_{\textbf{r}_{j}}(0)\big{]}\ket{\psi}\\\
&=\frac{1}{N_{d}}\sum_{|\rm{cl}|<=P_{c}}\sum_{\begin{subarray}{c}\textbf{r'}_{i},\textbf{r'}_{j}\in\rm{cl}\\\
|\textbf{r'}_{i}-\textbf{r'}_{j}|=d\end{subarray}}\bra{\psi_{\rm{cl}}}[\hat{S^{z}_{\textbf{r'}_{i}}}(t),\hat{S^{z}}_{\textbf{r'}_{j}}(0)]\ket{\psi_{\rm{cl}}}\\\
&C_{L}(d,t)=\frac{1}{N^{\prime}_{d}}\sum_{|\rm{cl}|=P}\sum_{\begin{subarray}{c}\textbf{r'}_{i},\textbf{r'}_{j}\in\rm{cl}\\\
|\textbf{r'}_{i}-\textbf{r'}_{j}|=d\end{subarray}}\bra{\psi_{\rm{cl}}}[\hat{S^{z}_{\textbf{r'}_{i}}}(t),\hat{S^{z}}_{\textbf{r'}_{j}}(0)]\ket{\psi_{\rm{cl}}}\end{split}$
(2)

where $\hat{S^{z}}_{\textbf{r}_{i}}$ and $\hat{S^{z}}_{\textbf{r}_{j}}$ are
spins at $\textbf{r}_{i}$ and $\textbf{r}_{j}$,
$|\textbf{r}_{i}-\textbf{r}_{j}|=d$, $\ket{\psi_{cl}}$ is a Haar random state
drawn from the Hilbert space spanned by all possible configuration states of
the cluster $\rm{cl}$ and $N_{d}$ is the number of pairs of spins in the
lattice $\mathcal{L}$ separated by distance $d$, excluding the pairs of spins
that belong to all clusters of size $|cl|>P_{c}$ as per our method ($P_{c}=18$
for this paper). $N^{\prime}_{d}$ is the number of pairs of spins separated by
distance $d$ in all clusters with $P$ plaquettes.
We analyze global information spreading represented by $C_{G}(d,t)$ in Eq.
(2), related to entanglement growth in the full system as well as local
information spreading $C_{L}(d,t)$, related to the same in a finite size
cluster.
Fig. 2a and b display a stark difference between the two cases in the
$|\lambda/J|\ll 1$ regime. There is a signal up to large distances with a
light cone like structure for information spread within a cluster, affirming
the ergodic nature of large finite size clusters in this regime Chakraborty
_et al._ (2022). However, for the whole lattice, there is an exponential decay
with a correlation length of around 3 to 4 lattice spacings, i.e. next-
nearest-neighbour plaquettes, implying information freezing and saturation of
entanglement growth. Such a qualitative difference is even more apparent in
the long-time average of the commutators with distance in Fig. 2c. Note that
such information freezing is present for all values of $|\lambda/J|$ and is a
generic feature of all disorder-free localized lattice gauge theories, where
kinetic constraints induced by static charges, cause real-space fragmentation.

Figure 2: Local and global information spreading quantified by
$J^{2}/\lambda^{2}|C_{L/G}(t)|$ for $|\lambda/J|=0.2$. a) Information
spreading within a 18 plaquette cluster b) Information freezing in the full
system c) Long time average vs distance for a) and b).

High temperature spectral functions We also calculate high-temperature
($T\gg\rm{max}(|\lambda|,|J|)$) spectral functions, which are relevant for
neutron scattering experiments on frustrated magnets that might realize a U(1)
QLM type effective theory. To do so, one employs arguments from dynamical
quantum typicality (DQT) Steinigeweg  _et al._ (2014); Elsayed and Fine (2013)
which state that a pure state drawn from a Haar-random distribution gives the
same expectation values as the entire statistical ensemble. One can then
define correlation functions as

$\begin{split}S^{zz}&(d,t,T)=\frac{1}{N}\sum_{\begin{subarray}{c}\textbf{r}_{i},\textbf{r}_{j}\in\mathcal{L}\\\
|\textbf{r}_{i}-\textbf{r}_{j}|=d\end{subarray}}\dfrac{\bra{\Phi_{\beta(t)}}\hat{S^{z}_{\textbf{r}_{i}}}\ket{\varphi_{\beta(t)}}}{\braket{\Phi_{\beta}(0)}{\Phi_{\beta}(0)}}+\epsilon\\\
&=\frac{1}{N^{\prime}_{d}}\sum_{\begin{subarray}{c}|cl|\\\
<=P_{c}\end{subarray}}\sum_{\begin{subarray}{c}\textbf{r}_{i},\textbf{r}_{j}\in
cl\\\
|\textbf{r}_{i}-\textbf{r}_{j}|=d\end{subarray}}\dfrac{\bra{\Phi_{\beta_{cl}}(t)}\hat{S^{z}_{\textbf{r}_{i}}}\ket{\varphi_{\beta_{cl}}(t)}}{\braket{\Phi_{\beta_{cl}}(0)}{\Phi_{\beta_{cl}}(0)}}+\epsilon\end{split}$
(3)

where $\ket{\Phi_{\beta_{cl}}(t)}=e^{-iH_{cl}t-\beta
H_{cl}/2}\ket{\psi_{cl}}$,
$\ket{\varphi_{\beta_{p}}(t)}=e^{-iH_{cl}t}\hat{S^{z}_{\textbf{r}_{j}}}e^{-\beta
H_{cl}/2}\ket{\psi_{cl}}$,$\beta=1/T$, $H_{cl}$ is the Hamiltonian in Eq. (1)
but with the sum over all plaquettes belonging to cluster $cl$, and $\epsilon$
denotes the error which is exponentially small in the number of thermally
occupied eigenstates. The typicality approach also works very well in a
localized phase for large Hilbert spaces Steinigeweg  _et al._ (2016).
Taking the space-time Fourier transform of Eq. (3), gives the transverse
component of the dynamical structure factor. As shown in Figs. 1 and 3a, there
are dominant delta-function like peaks at low frequencies amidst a background
of featureless response and a vanishing response as $\omega\rightarrow 0$.
This is the hallmark of disorder-free localization and distinguishes it from
paramagnetic response Conlon and Chalker (2009). Such response remains robust
on increasing system size and for all values of $\lambda$ and $J$ since the
dominant peaks arise from the formation of small clusters (see next section)
deep in the localized phase due to fragmentation induced by the static
charges. The density of such clusters saturates at a constant value as a
function of background charge density deep in the localized phase, smoothly
decays across the localization transition and then vanishes as a polynomial
deep in the delocalized phase, as in standard percolation. Crucially, on
crossing the localization transition, an infinite ergodic cluster emerges
which has a continuous spectral response Chakraborty  _et al._ (2022). Hence,
on crossing the transition point from the localized side, there is a shift in
weight from the discrete peaks to the continuous response.
However, the infinite cluster exhibits long-time hydrodynamic behaviour with
U(1) charge conservation. Assuming diffusion in 2D, i.e. a temporal decay of
$t^{-d/2}=1/t$, one would expect a logarithmic divergence of the spectral
response as $\omega\rightarrow 0$. The localized phase, however, has no such
contribution. Moreover, since the localization physics is dictated by the
small clusters, there is a vanishing spectral response as $\omega\rightarrow
0$ due to their discrete and non-degenerate spectra.

Besides the above global signatures, local spectral functions, such as the
temporal Fourier transform of onsite correlators, $S^{zz}_{0}(\omega)$,
spatially averaged over spins belonging to all clusters of a particular size,
evolves from a smooth curve with short bumps in the $|\lambda/J|\ll 1$ regime
to a discrete set of peaks in the $|\lambda/J|\gg 1$ regime, as in Fig. 3b.
Such a difference indicates a crossover from ergodic to non-ergodic behaviour
within large but finite clusters as predicted in Ref. Chakraborty  _et al._
(2022).

Figure 3:  Evidence of disorder-free localization from high temperature
($T=5\text{max}(|\lambda|,|J|)$) spatially averaged spectral functions. a) The
$q=0$ component of the transverse component of the dynamic structure factor, A
= max($|\lambda|$,$|J|$) and B = min($|\lambda|$,$|J|$), cutoff size
$P_{c}=18$ and $|\lambda/J|=5.0$ , (inset) same quantity as main plot, but for
$|\lambda/J|=0.2$ b) On-site spectral function averaged over all spins
belonging to all 18 plaquette clusters for $|\lambda/J|=5.0$ and (inset)
$|\lambda/J|=0.2$

Analytic estimate of dominant peaks in spectral functions: The position of the
discrete peaks observed in the spectral functions above can be accurately
estimated from the local physics of clusters of few plaquettes which form two
and three level systems. Consider a two level system, for example a single
flippable plaquette, the eigenvalues of such a system will be $-\lambda+J$ and
$-\lambda-J$. Such a system should have a discrete peak at $\omega=J/\pi$.
Similarly, for a two-plaquette system with 3 states ($11/10/01$; $1$ for
flippable and $0$ for not flippable), the transition between the closest
spaced pair is given by a discrete peak at $\omega=J^{2}/(\lambda\pi)$ in the
$|\lambda/J|\gg 1$ regime. The dominance of these peaks arising from such
local physics implies: i) On the global scale the dominant response deep in
the localized phase is from the small clusters and ii) In the $|\lambda/J|\gg
1$ regime, the physics of larger clusters is effectively dominated by local
motifs which form two- or three-level systems. Since within a cluster there is
no real-space fragmentation, the local motifs arise purely due to interference
effects. However, this apparent non-ergodicity is a finite-size effect with
dominant discrete peaks over a large finite-size regime, the spectral weight
of which eventually decreases on increasing size giving way to a continuous
spectrum in the infinite-size limit Chakraborty  _et al._ (2022). The former
contribution, however, remains discrete in the infinite-size limit, as the
presence of a finite density of charges will always induce real-space
fragmentation and form small clusters. These clusters have the same spectra
irrespective of their position, hence even after spatial averaging there
remains a sharp response in all regimes of the Hamiltonian.
Discussion and outlook: In this work we have presented numerical results for
spatially averaged spectral functions in 2D as characteristic evidence for
disorder-free localization. Local spectroscopic signatures have also been
proposed for conventional disordered MBL, as a way to distinguish between
strong and weak MBL Nandkishore  _et al._ (2014); Johri  _et al._ (2015).
However, most sharp spectral features are blurred out on spatial averaging,
with a soft spectral gap at zero frequency being the only distinct consequence
of both kinds of MBL in the presence of random disorder. For the disorder-free
localized phase, there are dominant sharp peaks in the dynamical structure
factor as well as vanishing response as $\omega\rightarrow 0$, both of these
can, in principle, be observed in spectroscopic experiments. Our work can be
translated to higher dimensions in which such sharp spectral features could
serve as evidence for quantum spin ice materials in 3D, where the analog of
the charges would be static defects. Our results on the unequal-time
commutators further motivate lattice gauge theories, possibly realized in
quantum simulators, as a fertile platform to probe information spreading and
chaos. The problem of lattice gauge theory with mobile matter adds a dynamic
charge component to the problem, hence an important future direction would be
to explore how, in that case, the signatures presented in this work would be
modified.

Acknowledgements.— We thank Hongzheng Zhao and Adam McRoberts for valuable
discussions. This project has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (Grant Agreement No. 853443), and M.H. further acknowledges support
by the Deutsche Forschungsgemeinschaft via the Gottfried Wilhelm Leibniz Prize
program. This work was in part supported by the Deutsche
Forschungsgemeinschaft under grants SFB 1143 (project-id 247310070) and the
cluster of excellence ct.qmat (EXC 2147, project-id 390858490).

## References

  * Fradkin (2013) E. Fradkin,  _Field theories of condensed matter physics_ (Cambridge University Press, 2013).
  * Moessner and Moore (2021) R. Moessner and J. E. Moore,  _Topological phases of matter_ (Cambridge University Press, 2021).
  * Moessner and Sondhi (2003) R. Moessner and S. L. Sondhi, Phys. Rev. B 68, 184512 (2003).
  * Hermele  _et al._ (2004) M. Hermele, M. P. A. Fisher, and L. Balents, Phys. Rev. B 69, 064404 (2004).
  * Smith  _et al._ (2017) A. Smith, J. Knolle, D. L. Kovrizhin, and R. Moessner, Phys. Rev. Lett. 118, 266601 (2017).
  * Brenes  _et al._ (2018) M. Brenes, M. Dalmonte, M. Heyl, and A. Scardicchio, Phys. Rev. Lett. 120, 030601 (2018).
  * Karpov  _et al._ (2021) P. Karpov, R. Verdel, Y.-P. Huang, M. Schmitt, and M. Heyl, Phys. Rev. Lett. 126, 130401 (2021).
  * Chakraborty  _et al._ (2022) N. Chakraborty, M. Heyl, P. Karpov, and R. Moessner, Phys. Rev. B 106, L060308 (2022).
  * Conlon and Chalker (2009) P. H. Conlon and J. T. Chalker, Phys. Rev. Lett. 102, 237206 (2009).
  * Papaefstathiou  _et al._ (2020) I. Papaefstathiou, A. Smith, and J. Knolle, Phys. Rev. B 102, 165132 (2020).
  * Li  _et al._ (2021) W.-H. Li, X. Deng, and L. Santos, Phys. Rev. Lett. 127, 260601 (2021).
  * Zhu and Heyl (2021) G.-Y. Zhu and M. Heyl, Phys. Rev. Research 3, L032069 (2021).
  * Hodson  _et al._ (2021) T. Hodson, J. Willsher, and J. Knolle, Phys. Rev. B 104, 045116 (2021).
  * Schweizer  _et al._ (2019) C. Schweizer, F. Grusdt, M. Berngruber, L. Barbiero, E. Demler, N. Goldman, I. Bloch, and M. Aidelsburger, Nature Physics 15, 1168 (2019).
  * Yang  _et al._ (2020) B. Yang, H. Sun, R. Ott, H.-Y. Wang, T. V. Zache, J. C. Halimeh, Z.-S. Yuan, P. Hauke, and J.-W. Pan, Nature 587, 392 (2020).
  * Zhou  _et al._ (2022) Z.-Y. Zhou, G.-X. Su, J. C. Halimeh, R. Ott, H. Sun, P. Hauke, B. Yang, Z.-S. Yuan, J. Berges, and J.-W. Pan, Science 377, 311 (2022).
  * Semeghini  _et al._ (2021) G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, _et al._ , Science 374, 1242 (2021).
  * Halimeh  _et al._ (2021) J. C. Halimeh, H. Zhao, P. Hauke, and J. Knolle, arXiv:2111.02427 (2021).
  * Chandrasekharan and Wiese (1997) S. Chandrasekharan and U.-J. Wiese, Nuclear Physics B 492, 455 (1997).
  * Horn (1981) D. Horn, Physics Letters B 100, 149 (1981).
  * Orland and Rohrlich (1990) P. Orland and D. Rohrlich, Nuclear Physics B 338, 647 (1990).
  * Wiese (2013) U.-J. Wiese, Annalen der Physik 525, 777 (2013).
  * Kitaev (2003) A. Y. Kitaev, Annals of Physics 303, 2 (2003).
  * Rokhsar and Kivelson (1988) D. S. Rokhsar and S. A. Kivelson, Phys. Rev. Lett. 61, 2376 (1988).
  * Moessner and Sondhi (2001) R. Moessner and S. L. Sondhi, Phys. Rev. B 63, 224401 (2001).
  * Shannon  _et al._ (2004) N. Shannon, G. Misguich, and K. Penc, Phys. Rev. B 69, 220403 (2004).
  * Halimeh  _et al._ (2022) J. C. Halimeh, P. Hauke, J. Knolle, and F. Grusdt, arXiv:2206.11273 (2022).
  * Metavitsiadis  _et al._ (2017) A. Metavitsiadis, A. Pidatella, and W. Brenig, Phys. Rev. B 96, 205121 (2017).
  * Maldacena  _et al._ (2016) J. Maldacena, S. H. Shenker, and D. Stanford, Journal of High Energy Physics 2016, 1 (2016).
  * Sala  _et al._ (2020) P. Sala, T. Rakovszky, R. Verresen, M. Knap, and F. Pollmann, Phys. Rev. X 10, 011047 (2020).
  * Khemani  _et al._ (2020) V. Khemani, M. Hermele, and R. Nandkishore, Phys. Rev. B 101, 174204 (2020).
  * von Keyserlingk  _et al._ (2018) C. W. von Keyserlingk, T. Rakovszky, F. Pollmann, and S. L. Sondhi, Phys. Rev. X 8, 021013 (2018).
  * Steinigeweg  _et al._ (2014) R. Steinigeweg, J. Gemmer, and W. Brenig, Phys. Rev. Lett. 112, 120601 (2014).
  * Elsayed and Fine (2013) T. A. Elsayed and B. V. Fine, Phys. Rev. Lett. 110, 070404 (2013).
  * Steinigeweg  _et al._ (2016) R. Steinigeweg, J. Herbrych, F. Pollmann, and W. Brenig, Phys. Rev. B 94, 180401 (2016).
  * Nandkishore  _et al._ (2014) R. Nandkishore, S. Gopalakrishnan, and D. A. Huse, Phys. Rev. B 90, 064203 (2014).
  * Johri  _et al._ (2015) S. Johri, R. Nandkishore, and R. N. Bhatt, Phys. Rev. Lett. 114, 117401 (2015).


#  Edge states of a diffusion equation in one dimension:
Rapid heat conduction to the heat bath

Shusei Makino1, Takahiro Fukui1, Tsuneya Yoshida2, and Yasuhiro Hatsugai2
1Department of Physics, Ibaraki University, Mito 310-8512, Japan  2Institute
of Physics, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8571,
Japan

###### Abstract

We propose a one-dimensional (1D) diffusion equation (heat equation) for
systems in which the diffusion constant (thermal diffusivity) varies
alternately with a spatial period $a$. We solve the time evolution of the
field (temperature) profile from a given initial distribution, by
diagonalising the Hamiltonian, i.e., the Laplacian with alternating diffusion
constants, and expanding the temperature profile by its eigenstates. We show
that there are basically phases with or without edge states. The edge states
affect the heat conduction around heat baths. In particular, rapid heat
transfer to heat baths would be observed in a short time regime, which is
estimated to be $t<10^{-2}$s for $a\sim 10^{-3}$m system and $t<1$s for $a\sim
10^{-2}$m system composed of two kinds of familiar metals such as titanium,
zirconium and aluminium, gold, etc. We also discuss the effective lattice
model which simplifies the calculation of edge states up to high energy. It is
suggested that these high energy edge states also contribute to very rapid
heat conduction in a very short time regime.

###### pacs:

##  I Introduction

The bulk-edge correspondence Hatsugai (1993a, b) is recognized as one of
fundamental concepts in physics Hasan and Kane (2010); Qi and Zhang (2011).
The original idea was proposed to answer the question of whether the quantum
Hall effect (QHE) is a bulk or edge-state property Hatsugai (1993a, b), but it
has now become a way to more broadly characterize topological insulators.
Indeed, it has been playing a central role in the development of topological
insulators and superconductors Kane and Mele (2005); Qi et al. (2008).
Remarkably, it has been extended to classical systems such as photonic
crystals Raghu and Haldane (2008); Haldane and Raghu (2008); Wang et al.
(2009); Ozawa et al. (2019); Kivshar (2019),  phononic systems Prodan and
Prodan (2009); Savin et al. (2010); Kane and Lubensky (2013); Kariyado and
Hatsugai (2015); Süsstrunk and Huber (2015); Chien et al. (2018); Yoshida and
Hatsugai (2019); Kivshar (2019),  electrical circuits Albert et al. (2015);
Lee et al. (2018); Helbig et al. (2020); Yoshida et al. (2020), hydrodynamics
Delplace et al. (2017); Sone and Ashida (2019), and so on.

Recently, the importance of edge states was also pointed out in diffusion
phenomena Yoshida and Hatsugai (2021). This implies that edge states would
have a large effect not only on propagating wave systems but also diffusive
systems such as heat conduction. In particular, in Ref. Yoshida and Hatsugai
(2021), a lattice model was introduced by coarsely discretizing the diffusion
equation, and their prediction was indeed observed experimentally Hu et al.
(2021); Qi et al. (2021). Here, our question is whether the heat conduction in
a continuous medium that is not necessarily protected by symmetry is affected
by edge states or not. This attempt is also aimed at understanding the role of
edge states in more mundane phenomena that are not necessarily of direct
relationship with topology.

In this paper, we examine the one-dimensional (1D) diffusion equation with
position-dependent diffusion constant, especially paying attention to the role
of edge states. In general, one dimensional systems would show anomalous heat
conduction due to anomalous diffusion Li and Wang (2003). However, we assume
systems governed by the normal diffusion as well as the normal Fourier or Fick
law, and hence, consider the conventional diffusion equation. Nevertheless,
the diffusion equation shows a characteristic behavior when systems allow edge
states.

To examine the effect of edge states, we propose 1D diffusion equation with
periodic array of two distinct diffusion constants.  Most coarse
discretization of such an equation with respect to space leads to the Su-
Schriffer-Heeger (SSH) Hamiltonian Su et al. (1979), which was already studied
in Yoshida and Hatsugai (2021). In this case, the relationship between edge
states and bulk topology is manifest. In this paper, we investigate the
opposite limit, a continuous equation for a continuous medium directly, in
which the bulk-edge correspondence is rather vague. Namely, a continuous
medium allows various boundary conditions, and the bulk topological invariant
protected by inversion symmetry does not necessarily guarantee the edge states
on boundaries which do not respect inversion symmetry. Nevertheless, even
without inversion symmetry, we find that there appear many edge states, and
the diffusion of initial heat distributions given near a boundary is
accelerated by edge states in a short time regime. These edge states are
surely related to those that appear for a boundary condition respecting
inversion symmetry. Thus, we claim that even if boundary conditions breaks
symmetry, some of them disappear but others remain persistently.

This paper is organized as follows. In Sec. II, we derive a diffusion equation
with an $x$-dependent diffusion constant. We regard it as an imaginary-time
Schrödinger equation and solve its initial value problem by using the complete
set of eigenstates. To this end, we solve the eigenvalue equation of the
Hamiltonian first by the Bloch techniques for the bulk system, and next by the
Fourier series expansion under Dirichlet boundary condition for a finite
system in contact with the heat baths. In the latter system, we find the
appearance of edge states localized at the boundary, which are not present in
the bulk system. Solving the initial value problem, we find rapid heat
transfer to the heat bath at the boundary, which can be understood as the
consequence of edge states. In Sec. III, we next discretize the diffusion
equation and derive an effective lattice model. This enables us to obtain the
edge states very simply up to very high energies. In Sec. IV, we give summary
and discussion including the experimental feasibility.

##  II Diffusion equation in one-dimension

We consider generic systems described by a one-dimensional diffusion equation
whose diffusion constant takes two values periodically in space. We assume
purely one-dimensional arrays of different materials, or some layered systems
stacked in one direction but uniform in the directions perpendicular to it.

For a while, we assume a diffusion equation with a position-dependent
diffusion constant $D(x)$ generically. Let $\phi(t,x)$ be the local field at
time $t$. We assume the Fick law

$\displaystyle j(t,x)=-D(x)\partial_{x}\phi(t,x),$ (1)

where $j(t,x)$ stands for the current density. Then, the continuity equation
reads

$\displaystyle\partial_{t}\phi(t,x)+\partial_{x}j(t,x)=0.$ (2)

These equations lead to the following diffusion equation,

$\displaystyle\partial_{t}\phi(t,x)-\partial_{x}D(x)\partial_{x}\phi(t,x)=0.$
(3)

In the following, we regard Eq. (3) as an imaginary-time Schrödinger equation
$\partial_{t}\psi+H\psi=0$ and consider the eigenstates of the Hamiltonian
operator

$\displaystyle
H\phi(x)\equiv-\partial_{x}D(x)\partial_{x}\phi(x)=\varepsilon\phi(x).$ (4)

Note that the Hamiltonian $H$ in Eq. (4) is Hermitian.

In what follows, we are mainly interested in the heat conduction. Then, to be
precise, it is necessary to take into account the $x$-dependence of two
constants characterizing the materials, i.e., the heat capacity and thermal
conductivity. As a result, the Fourier law with a $x$-dependent thermal
conductivity and the continuity equation with a $x$-dependent heat capacity
lead to more complicated equation than Eq. (4). In particular, the Hamiltonian
becomes non-Hermitian in general. However, for the sake of simplicity in this
paper, we consider only thermal diffusivity $D(x)$ as a $x$-dependent
parameter, ignoring the derivative of the heat capacity. Even with this
simpler equation, it is possible to analyze edge states, since they have
topological origin, albeit indirectly, as we will see. This guarantees the
robustness of edge states even if we have resort to any approximations, as
long as they are small.

More basically, if one derives the diffusion equation in inhomogeneous systems
microscopically, one starts with the Langevin equation, and derives the
Fokker-Planck equation. Here, it should be noted that the Langevin equation
allows several interpretations related to the examined process and the nature
of the noise, leading to different diffusion equations Leibovich and Barkai
(2019); dos Santos et al. (2020). Among them, Itô and Stratonovich types yield
non-hermitian Hamiltonians generically, whereas in this paper, we adopt the
Hänggi-Klimontovich type Hanggi (1982); Klimontovich (1990) which allows
hermitian Hamiltonian (4) for simplicity. We also mention that the Hamiltonian
(4) can be interpreted as a quantum mechanical kinetic Hamiltonian with a
position-dependent mass da Costa et al. (2020).

###  II.1 Bulk spectrum

We now assume $D(x+a)=D(x)$. Then, the Bloch theorem states that the
eigenfunction can be written as

$\displaystyle\phi(x)=e^{ikx}u_{k}(x),$ (5)

where $-\pi/a<k<\pi/a$ and the Bloch state $u_{k}(x)$ is periodic,
$u_{k}(x+a)=u_{k}(x)$. For $u_{k}(x)$, the eigenvalue equation (4) becomes

$\displaystyle-\left[(\partial_{x}D)(\partial_{x}+ik)+D(\partial_{x}+ik)^{2}\right]u_{k}=\varepsilon_{k}u_{k}.$
(6)

Now, let us expand the periodic functions $D(x)$ and $u_{k}(x)$ in the Fourier
series,

$\displaystyle D(x)=\sum_{n}e^{ik_{n}x}d_{n},\quad
u_{k}(x)=\sum_{n}e^{ik_{n}x}u_{k,n},$ (7)

with $k_{n}=2\pi n/a$ $(n=0,\pm 1,\pm 2,\cdots)$. Substituting these into Eq.
(6), we have

$\displaystyle\sum_{m}{\cal H}_{nm}u_{m}=\varepsilon u_{n},$ (8)

where the $k$-dependence of $u_{k,n}$ has been suppressed, and the Hamiltonian
is given by

$\displaystyle{\cal
H}_{nm}\equiv\left[k_{n-m}(k+k_{m})+(k+k_{m})^{2}\right]d_{n-m}.$ (9)

Thus, the Fourier coefficient $d_{n}$ determines the diffusion phenomena in
the present system.

---
Figure 1:  Schematic illustration of the system composed of alternating
diffusion constants $D_{\rm a}$ and $D_{\rm b}$. The dotted line indicates the
basic unit cell. The dashed line indicates the alternative unit cell with
inversion symmetry, which will be used for topological discussions in Sec.
III.1.2.

To be concrete, let us assume

$\displaystyle D(x)=\left\\{\begin{array}[]{ll}D_{\rm
a},&na<x<\left(n+\frac{1}{2}\right)a,\\\ D_{\rm
b},&\left(n+\frac{1}{2}\right)a<x<(n+1)a,\\\ \end{array}\right.$ (12)

as illustrated in Fig. 1. This can be expressed in the Fourier series such
that

$\displaystyle D(x)=\bar{D}+\frac{\epsilon}{i\pi}\sum_{n={\rm
odd}}\frac{1}{n}e^{ik_{n}x}\equiv\sum_{n}d_{n}e^{ik_{n}x},$ (13)

where $\bar{D}\equiv(D_{\rm a}+D_{\rm b})/2$ and $\epsilon=D_{\rm a}-D_{\rm
b}$.

In this paper, we assume one typical diffusion constant $d$, and consider
systems composed of materials with $d$ and $10d$. The effect of the ratio
$D_{\rm a}/D_{\rm b}$ will be discussed in Sec. IV. In Fig. 2, we show some
spectra of the Hamiltonian (9). The Hamiltonian of the uniform system with
$D_{\rm a}=D_{\rm b}\equiv d$ ($\epsilon=0$) is a simple Laplacian with the
dispersion,

$\displaystyle\varepsilon=dk^{2}\equiv\frac{1}{\tau}\tilde{k}^{2},\quad\tau\equiv\frac{a^{2}}{\pi^{2}d},\quad\tilde{k}\equiv\frac{ak}{\pi},$
(14)

as seen in Fig. 2 (a). Once $\epsilon\neq 0$ is introduced, spectrum has a gap
at multiples of $k=\pm\pi/a$. As a result, the lowest band of $D_{\rm
a}>D_{\rm b}$ system approaches the band of the uniform system with smaller
diffusion constant $D_{\rm a}=D_{\rm b}=d$, not that of uniform systems with
$D_{\rm a}=D_{\rm b}=10d$.

---
Figure 2:  Spectra of the model consisting of two kinds of diffusion
constants, Eq. (12). (a) $(D_{\rm a},D_{\rm b})=(d,d)$, and (b) $(D_{\rm
a},D_{\rm b})=(10d,d)\mbox{ or }(d,10d)$. The inset shows the low energy
spectra, where the black and red lines are the same as those in panels (a) and
(b). The black dotted line is for the uniform case $(D_{\rm a},D_{\rm
b})=(10d,10d)$.

So far we have discussed the gapped spectrum of the system with $(D_{\rm
a},D_{\rm b})=(10d,d)$. It may be needless to say that the spectrum for the
opposite series of the diffusion constant, $(D_{\rm a},D_{\rm b})=(d,10d)$, is
completely the same. Switch between $D_{\rm a}$ and $D_{\rm b}$ is achieved by
the transformations

  * (r)

reflection: $x\rightarrow-x$

  * (t)

translation: $x\rightarrow x\pm a/2$.

These symmetries are for the bulk systems, but approximately apply to the
finite systems, as seen in Sec. II.2.

###  II.2 Finite system with boundaries

Let us consider the system with length $L$ and impose the Dirichlet boundary
condition which implies, in the case of heat conduction, a system sandwiched
by heat baths at both ends. See Fig. 3. We assume that the field $\phi(x)$
describes the difference of the temperatures $\phi(x)=T(x)-T_{0}$, where
$T(x)$ and $T_{0}$ stand for the temperatures of the material and of the heat
baths, respectively. Therefore, we set $\phi(x)=0$ for $x<0,L<x$. We also
assume that within the temperature range of interest, the heat conduction in
each layer follows the usual diffusion equation with thermal diffusivity
$D_{\rm a}$ or $D_{\rm b}$.

---
Figure 3:  Schematic illustration of the finite size system of length $L$. At
the boundaries, the Dirichlet boundary condition $\phi(0)=\phi(L)=0$ is
imposed.

Then, we can expand $\phi(x)$ in the following Fourier series

$\displaystyle\phi(x)=\sum_{n>0}\phi_{n}\sin K_{n}x,$ (15)

where $\phi_{n}$ is real, and $K_{n}=\pi n/L$. Substituting Eq. (15) as well
as Eq. (7) into Eq. (4), we have the Hamiltonian

$\displaystyle{\cal
H}_{nm}^{L}=\bar{D}K_{n}^{2}\delta_{nm}{\color[rgb]{0,0,0}-}\pi_{nm}\frac{{\color[rgb]{0,0,0}2}\epsilon}{\pi}\left(\frac{2\pi}{a}S^{(1)}_{nm}K_{m}-S^{(2)}_{nm}K_{m}^{2}\right).$
(16)

Derivations and the definitions of matrix elements are given in Appendix A. In
numerical calculations below, we introduce a cutoff $n,m\leq 30$ for the
Hamiltonian (16).

---
|
Figure 4:  (a) Spectra of the finite size system with length $L=11a$. The
vertical grid lines are guide for the eye, separating the bulk bands. On some
of these lines, the edge states appear. Red dots are for $(D_{\rm a},D_{\rm
b})=(10d,d)$, whereas the black ones are $(D_{\rm a},D_{\rm b})=(d,d)$ for the
lower series and $(D_{\rm a},D_{\rm b})=(10d,10d)$ for the upper series. Lower
two panels show the eigenstate of (b) lowest edge state with $\varepsilon_{\rm
edge}\sim 3.22/\tau$ and (c) the 15th extended state with $\varepsilon_{\rm
edge}\sim 4.32/\tau$. Red lines and blue lines are for $(D_{\rm a},D_{\rm
b})=(10d,d)$ and $(D_{\rm a},D_{\rm b})=(d,10d)$, respectively.

In Fig. 4 (a), we show the spectra of the systems with length $L=11a$, solving
the eigenvalue equation ${\cal
H}^{L}\phi_{\alpha}=\varepsilon_{\alpha}\phi_{\alpha}$. Here, we emphasize
that the system with $(D_{\rm a},D_{\rm b})=(10d,d)$ gives exactly the same
spectrum, including the edge states, as the system with $(D_{\rm a},D_{\rm
b})=(d,10d)$: The differences lies in the wave functions. To see this, we show
in Fig. 4 (b), the eigenfunctions $\phi_{\alpha}(x)$ ($\alpha=11$),
corresponding to the lowest edge states. The red and blue curves are for the
systems with $(D_{\rm a},D_{\rm b})=(10d,d)$ and $(D_{\rm a},D_{\rm
b})=(d,10d)$, respectively, corresponding to Fig. 3(a) and 3(b). In the cases
$D_{\rm a}>D_{\rm b}$ and $D_{\rm a}<D_{\rm b}$, the edge states appear on the
left and right ends, respectively. For reference sake, we also show one
example of eigenfunction of the extended state in Fig. 4 (c). The
eigenfunctions of the edge states and extended states have a clear
distinction: While the extended states are related with each other by
reflection (r) as well as translation (t), the edge states are related only by
reflection (t).

###  II.3 Time evolution

Let $\phi_{n\alpha}$ be the normalized $\alpha$-th eigenstate of the
Hamiltonian ${\cal H}^{L}$ in Eq. (16). Since the Hamiltonian is Hermitian,
the eigenfunctions form a complete orthonormal basis,

$\displaystyle\sum_{n}\phi_{\alpha
n}^{T}\phi_{n\beta}=\delta_{\alpha\beta},\quad\sum_{\alpha}\phi_{n\alpha}\phi_{\alpha
m}^{T}=\delta_{nm},$ (17)

from which it follows

$\displaystyle\sum_{\alpha}\phi_{\alpha}(x)\phi_{\alpha}(y)=\frac{L}{2}\delta(x-y).$
(18)

Let $\phi_{\rm i}(x)$ be a given initial distribution of $\phi(t=0,x)$. Then,
the time evolution of $\phi(t,x)$ is induced by the Hamiltonian such that

$\displaystyle\phi(t,x)$ $\displaystyle=e^{-Ht}\phi_{\rm
i}(x)=\int_{0}^{L}dye^{-Ht}\delta(x-y)\phi_{\rm i}(y)$
$\displaystyle=\sum_{\alpha}e^{-\varepsilon_{\alpha}t}\phi_{\alpha}(x)(\phi_{\alpha},\phi_{\rm
i}),$ (19)

where $\displaystyle(\psi,\chi)=(2/L)\int_{0}^{L}dx\psi(x)\chi(x)$ stands for
the inner product of two real functions.

|
---|---
Figure 5:  The profiles $\phi(t=0.5\tau,x)$ starting from two kinds of initial
distributions (a) $(\ell_{1},\ell_{2})=(5,6)$ and (b) $(0,1)$ in Eq. (22)
denoted by green lines. The system size is $L=11a$. The red and blue lines are
for $(D_{\rm a},D_{\rm b})=(10d,d)$ and $(D_{\rm a},D_{\rm b})=(d,10d)$,
whereas black and black-dashed lines are for $(D_{\rm a},D_{\rm b})=(d,d)$ and
$(10d,10d)$, respectively.

In what follows, we give some numerical results for the initial state

$\displaystyle\phi_{\rm
i}(x)=\left\\{\begin{array}[]{ll}1&\ell_{1}a<x<\ell_{2}a\\\
0&\mbox{(others)}\end{array}\right.,$ (22)

where $\ell_{j}=0,1,2,\cdots$, implying that a sequence of several unit cells
has $\phi=1$ in the background of $\phi=0$.

####  II.3.1 Role of edge states

We first see that conventional bulk diffusion properties are indeed governed
by a smaller diffusion constant $d$ rather the larger one $10d$. To this end,
let us give an initial nonzero field at the center of the system,
$(\ell_{1},\ell_{2})=(5,6)$. Then, edge states should have nothing to do with
the diffusion as far as $t\ll L^{2}/d\sim 10^{3}\tau$. As shown in Fig. 5 (a),
starting from such an initial state, the profile of $\phi$ at finite $t$ is,
regardless of $D_{\rm a}<D_{\rm b}$ or $D_{\rm a}>D_{\rm b}$, similar to the
uniform system with smaller diffusion constant $(D_{\rm a},D_{\rm b})=(d,d)$
denoted by the black solid curve. The red and blue curves are slightly shifted
to the opposite directions. Physically, it is quite natural, while
mathematically, it is induced by (t) symmetry of the extended eigenfunctions
shown in Fig. 4 (c).

In contrast, if one gives a finite field $\phi_{\rm i}$ at the boundary unit
cell in contact with the heat bath, the diffusion of $\phi$ depends strongly
on $D_{\rm a}<D_{\rm b}$ or $D_{\rm a}>D_{\rm b}$. The red curve in Fig. 5 (b)
shows more rapid diffusion than the blue curve. On one hand, such a difference
seems quite natural, since the initial nonzero field flows directly to the
heat bath through the diffusion constant next to the bath. Indeed, in
$0<x<a/2$, the red (blue) curve is just on the dashed (solid) curve which is
the profile with the larger (smaller) diffusion constant. On the other hand,
from the point of view of the formula (19), the difference between the red and
blue curves is attributed to the edge states, since in the two systems, (i)
the spectra are exactly the same and (ii) bulk states are related with each
other by the translation (t), whereas only the edge states break the
translation (t).

|
---|---
Figure 6:  (a) Total internal energy and (b) edge internal energy as functions
of $t$. Colors of curves are the same as those in Fig. 5. Green dashed line
stands for $-\varepsilon_{\rm edge}t+{\rm const}$, where $\varepsilon_{\rm
edge}\sim 3.3/\tau$ is the energy of the lowest edge state, as indicated in
Fig. 4.

To clarify this more quantitatively, let us define “internal energy” by

$\displaystyle\Phi(t,x_{1},x_{2})=\int_{x_{1}}^{x_{2}}dx^{\prime}\phi(t,x^{\prime}).$
(23)

In Fig. 6, we show $\Phi_{\rm tot}(t)\equiv\Phi(t,0,L)$ and $\Phi_{\rm
cell}(t)\equiv\Phi(t,0,a)$ for the system in Fig. 5 (b). Figure 6 (a)
manifestly shows that the heat transfer to the heat bath is governed by the
leftmost diffusion constants. To see the effects of the edge states, we
calculate the diffusion from the unit cell at the left edge, shown in Fig. 6
(b). In the case of $(D_{\rm a},D_{\rm b})=(d,10d)$, there are no low-energy
edge states, implying that the diffusion occurs through the bulk extended
states only. Contrary to this, in the case of $(D_{\rm a},D_{\rm b})=(10d,d)$,
quite rapid diffusion is induced at very short time $t/\tau\ll 0.1$. This is
mainly due to high-energy edge states, as will be discussed in Sec. III.
Around $t/\tau\sim 0.1$, the lowest energy edge state with the energy
$\varepsilon\sim 3.2/\tau$ dominates the main diffusion which is suggested by
the coincidence with the green dashed line.

####  II.3.2 Rapid heat conduction

Since the localization length of the lowest edge state shown in Fig. 4 (b) is
about $2$ unit cells (See also Sec. III), the effect of the edge state would
be more pronounced if we start from a wider distribution of the initial
temperature near the boundary.

|
---|---
Figure 7:  (a) Same temperature profile as those in Fig. 5, but starting from
the initial distribution $(\ell_{1},\ell_{2})=(0,5)$. Panel (b) shows the
temperature profile of the unit cell at the left end in (a).

In Fig. 7, the temperature profile starting from the initial distribution
$(\ell_{1},\ell_{2})=(0,5)$ is shown. Remarkably, rapid heat transfer from the
left end to the heat bath is observed for the $(D_{\rm a},D_{\rm b})=(10d,d)$
system. Indeed, comparing the red curve with the dotted curve, we see that the
temperature at each point of the left end cell with $(D_{\rm a},D_{\rm
b})=(10d,d)$ is about half of a uniform system with $(D_{\rm a},D_{\rm
b})=(10d,10d)$. Since the temperature profile at $0<x<a/2$ is determined by
the balance between the heat escaping into the left heat bath and the heat
coming in from the right $a/2<x$, the temperature profile of the uniform
system with the dashed line which is higher than the red curves is due to
larger heat transfer from the right than that of the alternating system. This
fact, that heat far from the boundary influences the boundary temperature
profile, suggests that the temperature profile is determined by the extended
states in the uniform system. On the other hand, the temperature profile in
$0<x<a/2$ of the alternating system denoted by the red curve is not affected
by the heat far from the boundary and is quite similar to that in Fig. 5 (b),
as a result. This suggests that such temperature profile is dominated by the
edge states.

##  III Effective lattice model

So far we have discussed that the edge states modify the short-time heat
conduction near the boundaries in contact with heat baths. Eigenvalues and
eigenstates under the Dirichlet boundary condition have been obtained using
the Fourier series expansion. In this approach, one cannot judge an eigenstate
to be a bulk state or an edge state unless one has a look at the profile of
the eigenfunction. On the other hand, for tight-binding models, it is very
simple to obtain only the edge states separated from the bulk states Dwivedi
and Chua (2016); Duncan et al. (2018); Kunst et al. (2017, 2018, 2019a,
2019b); Pletyukhov et al. (2020); Fukui (2020). In this section, in order to
discuss the edge states of the diffusion equation more simply, we derive an
effective tight-binding Hamiltonian by discretizing the continuum Hamiltonian
in Eq. (4).

---
Figure 8:  Schematic illustration of the discretized system. The dashed-line
stands for the unit cell.

To this end, let us define the lattice labeled by $x=ja_{0}$ ($j=0,\pm
1,\cdots$), where $a_{0}$ is a lattice constant, and introduce

$\displaystyle\phi_{j}\equiv\phi(ja_{0}),\quad D_{j}\equiv D(ja_{0}),$ (24)

where $a=2a_{0}n$. On the lattice, we replace the differential operator in Eq.
(4) into difference operators defined by $\partial
f_{j}=(f_{j+1}-f_{j})/a_{0}$ and $\partial^{*}f_{j}=(f_{j}-f_{j-1})/a_{0}$.
The Hamiltonian in Eq. (4) can be discretized such that

$\displaystyle H\phi_{j}$
$\displaystyle\equiv-\partial^{*}D_{j}\partial\phi_{j}$
$\displaystyle=\left[-D_{j-1}\phi_{j-1}+(D_{j-1}+D_{j})\phi_{j}-D_{j}\phi_{j+1}\right]/a_{0}^{2}.$
(25)

It may be convenient to regard $2n$ sites as a unit cell, as illustrated in
Fig. 8. Namely, we label the sites such that

$\displaystyle\phi_{j}=\phi_{\alpha,J},\quad j=2nJ+\alpha.$ (26)

Then, the Hamiltonian becomes $2n\times 2n$ matrix operator labelled by
$\alpha$.

Since the unit cell can be chosen by an arbitrary $2n$ sequential sites, let
us define more generic Hamiltonian such that

$\displaystyle
H=\frac{1}{a_{0}^{2}}\left(\begin{array}[]{ccccccc}D_{2n}+D_{1}&-D_{1}&&&&&-D_{2n}\delta^{*}\\\
-D_{1}&D_{1}+D_{2}&-D_{2}&&&&\\\ &-D_{2}&&&&&\\\ &&&\ddots&&&\\\
&&&&&-D_{2n-2}&\\\ &&&&-D_{2n-2}&D_{2n-2}+D_{2n-1}&-D_{2n-1}\\\
-D_{2n}\delta&&&&&-D_{2n-1}&D_{2n-1}+D_{2n}\end{array}\right),$ (34)

where $\delta$ and $\delta^{*}$ stand for the forward and backward shift
operators defined by $\delta f_{J}=f_{J+1}$ and $\delta^{*}f_{J}=f_{J-1}$. The
bulk spectrum can be determined by setting $\delta\rightarrow e^{ik}$ and
$\delta^{*}\rightarrow e^{-ik}$. The wavenumber $k$ of this lattice model,
expressed in terms of the wavenumber $k$ of the Bloch state in Eq. (5),
implies $ak$. Here, we choose $D_{1},D_{2},\cdots,D_{n}=D_{\rm a}$, and
$D_{n+1},\cdots,D_{2n}=D_{\rm b}$ for the effective lattice Hamiltonian in
Sec. II.

###  III.1 Edge states

For tight-binding models, the transfer matrix method is the standard technique
to derive the edge states embedded in the bulk Hatsugai (1993a, b).
Alternatively, the edge states located at the left end and at the right end
can be separately computed based on the Hermiticity of the Hamiltonian, as
proposed in Ref. Fukui, 2020. In this paper, we use the latter method to
compute the edge states, which may be useful for Hamiltonians with large
dimensions such as the discretized Hamiltonian of the present system. This
method is briefly summarized in Appendix B.

####  III.1.1 Edge state Hamiltonian

As discussed in Appendix B, the edge states localized at the left end are
described by the $(2n-1)\times(2n-1)$ Hamiltonian obtained by neglecting the
$2n$th row and column of the Hamiltonian (34),

$\displaystyle H_{\rm
es}=\frac{1}{a_{0}^{2}}\left(\begin{array}[]{cccccc}D_{2n}+D_{1}&-D_{1}&&&&\\\
-D_{1}&D_{1}+D_{2}&-D_{2}&&&\\\ &-D_{2}&&&&\\\ &&&\ddots&&\\\
&&&&&-D_{2n-2}\\\ &&&&-D_{2n-2}&D_{2n-2}+D_{2n-1}\\\ \end{array}\right).$ (41)

The eigenstates of this Hamiltonian are not necessarily true edge states of
the Hamiltonian (34): We can choose the edge states by requiring that the
localization length $1/\kappa$ should be positive,

$\displaystyle
e^{-\kappa}=\left|\frac{D_{2n-1}\chi_{2n-1}}{D_{2n}\chi_{1}}\right|<1.$ (42)

The edge states at the right end can be also derived in a similar way Fukui
(2020).

---
Figure 9:  (a) Energies of edge states localized at the left end [eigenvalues
of the edge state Hamiltonian (41)] with their momenta defined by $k$ in Eq.
(70) are denoted by colored dots on the background of the bulk spectra for the
model with $(D_{\rm a},D_{\rm b})=(10d,d)$ denoted by lines. Red (gray) dots
show the energies of the edge states (not) satisfying the localization
condition (42). (b) Same as (a) but for the states localized at the right end.
There are no states satisfying the localization condition, implying no edge
states at the right end. Inset shows the states in (a) in the low energy
region.

In Fig. 9, we plot the eigenvalues of the left and right edge state
Hamiltonians by colored dots on the background of the bulk bands computed by
the Hamiltonian (34). As in the case of the continuum model in Sec. II, the
lattice model with $D_{\rm a}>D_{\rm b}$ shows the edge states only at the
left end. Table 1 shows the localization length of the edge states. The lowest
edge state has $\sim 2$ unit cell localization length, implying that the
boundary heat conduction is not affected by the initial heat in $x>2$, as
indeed observed in Sec. II.3.2.

---
Figure 10:  Same as Fig. 9 up to very high energies. We use $n=100$, so that
there can potentially be 199 edge states.

In addition to the edge states in Fig. 9, there can be many other edge states
up to very high energies. We extend Fig. 9 up to very high energies in Fig.
10. The system with $(D_{\rm a},D_{\rm b})=(10d,d)$ is indeed the case,
whereas the system with $(D_{\rm a},D_{\rm b})=(d,10d)$ never show edge states
in all energy regime. Such difference affects the heat conduction at very
short time discussed in Sec. II.3.1.

Table 1: Energies and localization lengths $1/\kappa$ of the edge states denoted by the red dots in Fig. 9. State # is from the lowest to higher energy. State # | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
$\varepsilon$ | 3.02 | 8.59 | 17.83 | 59.10 | 77.59 | 103.89
$1/\kappa$ | 2.20 | 0.88 | 3.23 | 1.89 | 0.98 | 7.36

####  III.1.2 Topological properties

In this section, we argue that the edge states obtained so far have intimate
relationship with the topological property of the bulk system. For the
continuum model, we have chosen the area surrounded by the dotted line in
Fig.1 as the unit cell. Correspondingly, we have chosen the unit cell of the
lattice model to match the continuum model. In this case, the model has broken
inversion symmetry. However, the unit cell in Fig. 1 and corresponding lattice
model has inversion symmetry. Then, the Berry phase of each band is quantized
as $0$ or $\pi$, which serves as the topological invariant Ryu and Hatsugai
(2002). It should be noted that bulk spectrum does not depend on the choice of
a specific unit cell. The choice of a unit cell means the choice of the
boundaries.

---
Figure 11:  Same as Fig. 9, only the unit cell is changed to have inversion
symmetry, as illustrated in Fig. 1. The solid lines and dashed lines stand for
the bulk bands with Berry phase $\pi$ and $0$, respectively. The red and blue
dots stand for the edge states at the left and right ends, respectively,
satisfying the localization condition. The top two panels (a) are for $(D_{\rm
a},D_{\rm b})=(10d,d)$, while the bottom two panels (b) are for $(D_{\rm
a},D_{\rm b})=(d,10d)$.

In Fig. 11, we show the spectra of the lattice model with inversion symmetry.
Each bulk band is distinguished as a solid curve (Berry phase $\pi$) or a
dashed curve (Berry phase $0$). With inversion symmetry, the left end and the
right end are completely the same, so that the same edge states appear on both
ends. In this system, we can check the bulk-edge correspondence, namely, an
edge state appears if the sum of Berry phases of the lower bulk bands are
$\pi$ modulo $2\pi$. Comparing Figs. 9, 10, and 11, we conclude that although
the edge states in Figs. 9 and 10 are not directly protected by inversion
symmetry, their origin lies in the symmetry-protected edge states: Symmetry
breaking boundaries make some of them disappear, whereas others remain.

##  IV Summary and discussion

We have examined the 1D diffusion equation with the SSH-like alternating
diffusion constants. Since the diffusion equation can be regarded as an
imaginary-time Schödinger equation, we have calculated the eigenvalues and
eigenstates of the Hamiltonian which is the spatial part of the diffusion
equation. We have shown that the SSH-like structure yields spectral gaps, and
if boundaries are introduced, there appear edge states within the bulk gaps.
These edge states have a significant effect on the diffusion process in a
short-time regime near boundaries: Rapid heat conduction to the heat bath is
expected if a boundary allows edge states.

---
Figure 12:  Energies of edge states at the left end denoted by red dots. The
background liens are bulk bands. Panels (a), (b), and (c) are the data
obtained for $D_{\rm a}/D_{\rm b}=3$, $30$, and $100$, respectively. The
energies and localization lengths of the lowest edge states are (a) 1.78 and
2.17, (b) 3.64 and 3.92, and (c) 3.85 and 6.60.

Finally, let us discuss the experimental feasibility. The model includes
basically two parameters, $D_{\rm a}$ and $D_{\rm b}$. We assume $D_{\rm
a}>D_{\rm b}\equiv d$, for simplicity. In this case, the model has edge states
at the left end, which are characterized well by rearranged two parameters
$\tau=a^{2}/(\pi^{2}d)$ and $D_{\rm a}/D_{\rm b}$. The key parameter is
$\tau=a^{2}/(\pi^{2}d^{2})$ rather than their ratio $D_{\rm a}/D_{\rm b}$. To
see this, we show the bulk and edge states energies in Fig. 12, changing the
value of $D_{\rm a}$ such that $D_{\rm a}/D_{\rm b}=3,30,100$ (and $10$ in the
inset in Fig. 9). The energies and the localization lengths vary, but
qualitative difference is rather small.

Next, let us consider the material-dependence of the parameter $\tau$. The
titanium and zirconium are candidates of materials with smaller thermal
diffusivities, $d\sim 7\times 10^{-6}$ m2/s, and hence, fixing the length of
the unit cell as $a=10^{-3}$ m, we have $\tau=1.4\times 10^{-2}$ s. Combined
with materials having larger thermal diffusivities such as carbon ($D_{\rm
a}=2.3d$), iron ($2.5d$), lead ($3.4d$), aluminum ($12d$), and gold ($14d$),
rapid heat transfer through edge states discussed in Sec. II.3.2 would be
observed up to $t\sim 0.5\tau\sim 10^{-2}$ s for samples with $a=10^{-3}$ m,
and $t\sim 1$ s for samples with $a\sim 10^{-2}$ m.

###### Acknowledgements.

This work was supported in part by Grants-in-Aid for Scientific Research
Numbers 17K05563, 17H06138, JP21K13850, and JP20H04627 from the Japan Society
for the Promotion of Science.

##  Appendix A Fourier series expansion under the Dirichlet boundary condition

According to the Dirichlet boundary condition (15), It may be convenient to
rewrite Eq. (13) with respect to $\sin k_{l}$,

$\displaystyle D(x)=\bar{D}+\frac{2\epsilon}{\pi}\sum_{l={\rm
odd}>0}\frac{1}{l}\sin k_{l}{\color[rgb]{0,0,0}x}.$ (43)

Substituting Eq. (43) as well as (15) into Eq. (4), and using

$\displaystyle\int_{0}^{L}\sin K_{n}x\sin K_{m}xdx=(L/2)\delta_{nm},$
$\displaystyle\int_{0}^{L}\sin K_{n}x\cos
K_{m}xdx=\frac{2K_{n}}{K_{n}^{2}-K_{m}^{2}}\pi_{nm},$ (44)

where $\pi_{nm}=[1-(-)^{n+m}]/2$, we obtain the Hamiltonian (16). The matrix
elements associated with the couplings with $d_{l}$ are defined by

$\displaystyle S^{(1)}_{nm}=\frac{2}{\pi}\sum_{l={\rm
odd}>0}\left[\frac{n+m}{(n+m)^{2}-(2Nl)^{2}}+\frac{n-m}{(n-m)^{2}-(2Nl)^{2}}\right],$
$\displaystyle S^{(2)}_{nm}=\frac{4N}{\pi}\sum_{l={\rm
odd}>0}\left[\frac{1}{(n+m)^{2}-(2Nl)^{2}}-\frac{1}{(n-m)^{2}-(2Nl)^{2}}\right],$
(45)

with $N=L/a$.

##  Appendix B Edge state Hamiltonian

The eigenvalue equation of the Hamiltonian (25) can be written as

$\displaystyle\left({\cal K}\delta^{*}+{\cal V}+{\cal
K}^{\dagger}\delta\right)\phi_{J}=\varepsilon\phi_{J},$ (46)

where $\cal K$ and $\cal K^{\dagger}$ are matrices proportional to
$\delta^{*}$ and $\delta$, and $\cal V$ is the remaining part of the
Hamiltonian in Eq. (34). As noted in Sec. II, the continuum Hamiltonian is
Hermitian. As a result, the discretized Hamiltonian (34) is also Hermitian as
far as the system is defined on the infinite line. However, if the system has
a boundary, the Hamiltonian is not necessarily Hermitian because of the
boundary term due to the summation by parts. As proposed in Ref. Fukui, 2020,
imposing the Hermiticity condition naturally leads to a theoretical framework
that allows us to discuss only the edge states in isolation from the bulk
states.

Assume that the system has a boundary at $J=1$ and is defined on the semi-
infinite line $J\geq 1$. Then, Hermiticity is guaranteed by introducing the
reference state $\phi_{0}$ and requiring

$\displaystyle{\cal K}\phi_{0}=0.$ (47)

It is readily find that the state $\phi_{0}^{T}=(\chi,0)$ satisfies Eq. (47),
where $\chi$ is a vector with $2n-1$ components. Taking $\phi_{0}$ as a
reference state, we can obtain the edge states at the left end by assuming the
Bloch-like sates $\phi_{J}=\phi_{0}e^{iKJ}$ with $K=k+i\kappa$. Note that $k$
and $1/\kappa$ are the momentum and the localization length of the edge state,
respectively. The eigenvalue equation then becomes

$\displaystyle\frac{1}{a_{0}^{2}}\left(\begin{array}[]{cccccc|c}D_{2n}+D_{1}&-D_{1}&&&&&-D_{2n}e^{-iK}\\\
-D_{1}&D_{1}+D_{2}&-D_{2}&&&&\\\ &-D_{2}&&&&&\\\ &&&\ddots&&&\\\
&&&&&-D_{2n-2}&\\\ &&&&-D_{2n-2}&D_{2n-2}+D_{2n-1}&-D_{2n-1}\\\ \hline\cr-
D_{2n}e^{iK}&&&&&-D_{2n-1}&D_{2n-1}+D_{2n}\end{array}\right)\left(\begin{array}[]{c}\\\
\\\ \\\ \chi\\\ \\\ \\\ \hline\cr 0\\\
\end{array}\right)=\varepsilon\left(\begin{array}[]{c}\\\ \\\ \\\ \chi\\\ \\\
\\\ \hline\cr 0\\\ \end{array}\right).$ (69)

Therefore, the edge states are eigenstates of the upper $(2n-1)$ dimensional
matrix in Eq. (41) with the constraint,

$\displaystyle
e^{iK}=e^{ik}e^{-\kappa}=-\frac{D_{2n-1}\chi_{2n-1}}{D_{2n}\chi_{1}},$ (70)

which follows from the condition to be satisfied by the lowest component in
Eq. (69). This equation gives the localization length $1/\kappa$ and the
momentum $k$ of the edge states by the wave functions.

## References

  * Hatsugai (1993a) Y. Hatsugai, Physical Review Letters 71, 3697 (1993a), URL http://link.aps.org/doi/10.1103/PhysRevLett.71.3697.
  * Hatsugai (1993b) Y. Hatsugai, Physical Review B 48, 11851 (1993b), URL https://link.aps.org/doi/10.1103/PhysRevB.48.11851.
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Reviews of Modern Physics 82, 3045 (2010).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Reviews of Modern Physics 83, 1057 (2011).
  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Physical Review Letters 95, 146802 (2005).
  * Qi et al. (2008) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Physical Review B 78, 195424 (2008).
  * Raghu and Haldane (2008) S. Raghu and F. D. M. Haldane, Physical Review A 78, 033834 (2008), URL https://link.aps.org/doi/10.1103/PhysRevA.78.033834.
  * Haldane and Raghu (2008) F. D. M. Haldane and S. Raghu, Physical Review Letters 100, 013904 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.100.013904.
  * Wang et al. (2009) Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, Nature 461, 772 EP (2009), URL https://doi.org/10.1038/nature08293.
  * Ozawa et al. (2019) T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, et al., Reviews of Modern Physics 91, 015006 (2019), URL https://link.aps.org/doi/10.1103/RevModPhys.91.015006.
  * Kivshar (2019) Y. Kivshar, Low Temperature Physics 45, 1026 (2019), URL https://doi.org/10.1063/1.5121273.
  * Prodan and Prodan (2009) E. Prodan and C. Prodan, Physical Review Letters 103, 248101 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.103.248101.
  * Savin et al. (2010) A. V. Savin, Y. S. Kivshar, and B. Hu, Physical Review B 82, 195422 (2010), URL https://link.aps.org/doi/10.1103/PhysRevB.82.195422.
  * Kane and Lubensky (2013) C. L. Kane and T. C. Lubensky, Nature Physics 10, 39 EP (2013), URL https://doi.org/10.1038/nphys2835.
  * Kariyado and Hatsugai (2015) T. Kariyado and Y. Hatsugai, Scientific Reports 5, 18107 EP (2015), URL https://doi.org/10.1038/srep18107.
  * Süsstrunk and Huber (2015) R. Süsstrunk and S. D. Huber, Science 349, 47 (2015), URL https://science.sciencemag.org/content/sci/349/6243/47.full.pdf.
  * Chien et al. (2018) C.-C. Chien, K. A. Velizhanin, Y. Dubi, B. R. Ilic, and M. Zwolak, Physical Review B 97, 125425 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.125425.
  * Yoshida and Hatsugai (2019) T. Yoshida and Y. Hatsugai, Physical Review B 100, 054109 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.054109.
  * Albert et al. (2015) V. V. Albert, L. I. Glazman, and L. Jiang, Physical Review Letters 114, 173902 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.173902.
  * Lee et al. (2018) C. H. Lee, S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, and R. Thomale, Communications Physics 1, 39 (2018), URL https://doi.org/10.1038/s42005-018-0035-2.
  * Helbig et al. (2020) T. Helbig, T. Hofmann, S. Imhof, M. Abdelghany, T. Kiessling, L. W. Molenkamp, C. H. Lee, A. Szameit, M. Greiter, and R. Thomale, Nature Physics 16, 747 (2020), URL https://doi.org/10.1038/s41567-020-0922-9.
  * Yoshida et al. (2020) T. Yoshida, T. Mizoguchi, and Y. Hatsugai, Physical Review Research 2, 022062 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.022062.
  * Delplace et al. (2017) P. Delplace, J. B. Marston, and A. Venaille, Science 358, 1075 (2017), URL https://science.sciencemag.org/content/sci/358/6366/1075.full.pdf.
  * Sone and Ashida (2019) K. Sone and Y. Ashida, Physical Review Letters 123, 205502 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.205502.
  * Yoshida and Hatsugai (2021) T. Yoshida and Y. Hatsugai, Scientific Reports 11, 888 (2021).
  * Hu et al. (2021) H. Hu, S. Han, Y. Yang, D. Liu, H. Xue, G.-G. Liu, Z. Cheng, Q. J. Wang, S. Zhang, B. Zhang, et al., _Observation of topological edge states in thermal diffusion_ (2021), eprint 2107.05811.
  * Qi et al. (2021) M. Qi, D. Wang, P.-C. Cao, X.-F. Zhu, C.-W. Qiu, H. Chen, and Y. Li, _Localized heat diffusion in topological thermal materials_ (2021), eprint 2107.05231.
  * Li and Wang (2003) B. Li and J. Wang, Physical Review Letters 91, 044301 (2003).
  * Su et al. (1979) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Physical Review Letters 42, 1698 (1979), URL https://link.aps.org/doi/10.1103/PhysRevLett.42.1698.
  * Leibovich and Barkai (2019) N. Leibovich and E. Barkai, Physical Review E 99, 042138 (2019), URL https://link.aps.org/doi/10.1103/PhysRevE.99.042138.
  * dos Santos et al. (2020) M. A. F. dos Santos, V. Dornelas, E. H. Colombo, and C. Anteneodo, Physical Review E 102, 042139 (2020), URL https://link.aps.org/doi/10.1103/PhysRevE.102.042139.
  * Hanggi (1982) P. Hanggi, Physical Review A 25, 1130 (1982), URL https://link.aps.org/doi/10.1103/PhysRevA.25.1130.
  * Klimontovich (1990) Y. L. Klimontovich, Physica A: Statistical Mechanics and its Applications 163, 515 (1990), URL https://www.sciencedirect.com/science/article/pii/037843719090142F.
  * da Costa et al. (2020) B. G. da Costa, I. S. Gomez, and M. A. F. dos Santos, EPL (Europhysics Letters) 129, 10003 (2020), URL http://dx.doi.org/10.1209/0295-5075/129/10003.
  * Dwivedi and Chua (2016) V. Dwivedi and V. Chua, Physical Review B 93, 134304 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.134304.
  * Duncan et al. (2018) C. W. Duncan, P. Öhberg, and M. Valiente, Physical Review B 97, 195439 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.195439.
  * Kunst et al. (2017) F. K. Kunst, M. Trescher, and E. J. Bergholtz, Physical Review B 96, 085443 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.085443.
  * Kunst et al. (2018) F. K. Kunst, G. van Miert, and E. J. Bergholtz, Physical Review B 97, 241405 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.241405.
  * Kunst et al. (2019a) F. K. Kunst, G. van Miert, and E. J. Bergholtz, Physical Review B 99, 085426 (2019a), URL https://link.aps.org/doi/10.1103/PhysRevB.99.085426.
  * Kunst et al. (2019b) F. K. Kunst, G. van Miert, and E. J. Bergholtz, Physical Review B 99, 085427 (2019b), URL https://link.aps.org/doi/10.1103/PhysRevB.99.085427.
  * Pletyukhov et al. (2020) M. Pletyukhov, D. M. Kennes, J. Klinovaja, D. Loss, and H. Schoeller, Physical Review B 101, 165304 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.101.165304.
  * Fukui (2020) T. Fukui, Physical Review Research 2, 043136 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.043136.
  * Ryu and Hatsugai (2002) S. Ryu and Y. Hatsugai, Physical Review Letters 89, 077002 (2002), URL http://link.aps.org/doi/10.1103/PhysRevLett.89.077002.