# Nonlinear non-Hermitian higher-order topological laser

Motohiko Ezawa  Department of Applied Physics, University of Tokyo, Hongo
7-3-1, 113-8656, Japan

###### Abstract

We investigate topological lasers in combination of nonlinear, non-Hermitian
and topological lattice systems based on a quench dynamics starting from one
site. We consider explicitly the topological laser in the Su-Schrieffer-Heeger
(SSH) model with two topological edge states and the second-order topological
laser in the breathing Kagome lattice with three topological corner states.
Once we stimulate any one site, after a delay, all sites belonging to the
topological edge or corner states begin to emit stable laser light depending
on the density of states, although no wave propagation is observed from the
stimulated site. It is intriguing that the profile of topological edge or
corner states is observable by measuring the intensity of lasing. The
phenomenon occurs due to a combinational effect of linear non-Hermitian loss
terms and nonlinear non-Hermitian gain terms in the presence of the
topological edge or corner states.

##  I Introduction

Topological physics is one of the most essential concepts found in recent
fundamental physicsHasan ; Qi . Recently, it is ubiquitously found in various
systems in photonicKhaniPhoto ; Hafe2 ; Hafezi ; WuHu ; TopoPhoto ; Ozawa16 ;
Ley ; KhaniSh ; Zhou ; Jean ; Ota18 ; Ozawa ; Ota19 ; OzawaR ; Hassan ; Ota ;
Li ; Yoshimi ; Kim ; Iwamoto21 , acousticProdan ; TopoAco ; Berto ; Xiao ; He
; Abba ; Xue ; Ni ; Wei ; Xue2 , mechanicalLubensky ; Chen ; Nash ; Paul ; Sus
and electric circuitTECNature ; ComPhys ; Hel ; Lu ; YLi ; EzawaTEC ; EzawaLCR
; EzawaSkin  systems. Among them, topological photonics is most extensively
studied theoretically and experimentally. One of the reasons is that it is
possible to observe real-time and real-space dynamics. Another merit is that
topological photonics has opened a new field of topological physics, i.e. non-
Hermitian topology and nonlinear topology and their combination. Non-Hermitian
effects are introduced by a loss and a gain of photons. On the other hand,
nonlinear effects are introduced by the Kerr effect or a stimulated emission
effect.

A topological laser is a prominent application of topological physicsHarrari ;
Bandres ; Schome ; Weimann ; Jean ; Ota18 ; Parto ; Zhao ; Malzard ;
MalzardOpt ; Zhong . It utilizes topological edge states for the coherent
laser emission. Thanks to the topological protection, the topological laser is
robust against the randomness and the defects of the sample, which is
favorable for future laser applications. A topological laser is an ideal play
ground to investigate nonlinear non-Hermitian topological physics. The loss of
photons and gain from stimulated emissions constitute the non-Hermitian terms.
The nonlinear effect is included in the gain term, which represents the
saturation of the gain.

Higher-order topological phases are extension of topological phasesFan ;
Science ; APS ; Peng ; Lang ; Song ; Bena ; Schin ; EzawaKagome . There emerge
topological corner states in the second-order topological phase. A simplest
example is given by the breathing Kagome lattice, where three topological
corner states appear in triangle geometry. This model is a natural
generalization of the Su-Schrieffer-Heeger (SSH) model to two dimensions.
Since the model requires only positive hoppings, it is realized in various
systems including photonicHassan ; Li , acousticXue ; Ni  and electric circuit
systemsWuKagome .

In this paper, we analyze a quench dynamics of a nonlinear-non-Hermitian
topological laser and a higher-order topological laser by stimulating any one
site. We study explicitly the SSH model for a topological laser and the
breathing Kagome model for a higher-order topological laser. In the SSH model
with two topological edge states, once we stimulate any one site, after a
delay, all sites belonging to the two topological edge states begin to emit
stable laser light depending on the density of states (DOS), although no wave
propagation is observed from the stimulated site. Thus, the profile of a
topological edge state together with the DOS is observable by laser intensity.
The strength of the laser light is identical for the two edges because of
reflection symmetry. This is also the case for the second-order topological
laser in the breathing Kagome lattice with three topological corner states.
Here, the system has trigonal symmetry. The phenomenon occurs due to a
combinational effect of linear non-Hermitian loss terms and nonlinear non-
Hermitian gain terms in the presence of the topological edge or corner states.

##  II Topological laser

###  II.1 Model

We consider a coupled-ring system made of active resonatorsHarrari  . The
dynamics of a laser system is governed byHarrari

$i\frac{d\psi_{n}}{dt}=\sum_{nm}M_{nm}\psi_{m}-i\gamma\left(1-\xi\frac{P_{n}}{1+\left|\psi_{n}\right|^{2}/\eta}\right)\psi_{n},$
(1)

where $\psi_{n}$ is the amplitudes of the site $n$, where $n=1,2,3,\cdots,N$
in the system composed of $N$ sites; $M_{nm}$ describes a hopping matrix;
$\gamma$ represents the loss in each resonator; $\gamma\xi$ represents the
amplitude of the optical gain via stimulated emission; $\eta$ represents the
nonlinearity; $P_{n}$ stands for the spatial profile of the pump. The system
turns into the linear model in the limit $\eta\rightarrow\infty.$ On the other
hand, $\gamma$ controls the non-Hermicity. The system turns into a Hermitian
model for $\gamma=0$. We call the term proportional to $\gamma$ the loss term
and the term proportional to $\gamma\xi$ the nonlinear gain term.

We take

$P_{n}=\sum_{\bar{n}}\delta_{n,\bar{n}},$ (2)

where $\bar{n}$ runs over the edge or corner sites. Namely, optical gains are
introduced only at the edge or corner sites.

We are interested in the case where $M_{nm}$ represents a tight-binding model
possessing a topological phase. We explicitly consider the SSH model
illustrated in Fig.1(a1)$\sim$(a3), where $\bar{n}$ takes values at the left
and right edges, and the breathing Kagome model illustrated in
Fig.1(b1)$\sim$(b3), where $\bar{n}$ takes values at the top, bottom-left and
bottom-right corners. See Appendix for the topological charges in these
models.

It is possible to solve Eq.(1) numerically for explicit system parameters, as
we do later. However, to reach a deeper understanding of the phenomena, an
analytical study is indispensable. Since this is impossible for general system
parameters, we make an analytical study for special cases.

Figure 1: (a1)$\sim$(a3) Illustration of a dimerized lattice with (a1)
$\kappa_{A}=0$, (a2) $\kappa_{A}\kappa_{B}\neq 0$ and (a3) $\kappa_{B}=0$.
(b1)$\sim$(b3) Illustration of a breathing Kagome lattice with (b1)
$\kappa_{A}=0$, (b2) $\kappa_{A}\kappa_{B}\neq 0$ and (b3) $\kappa_{B}=0$. A
line (triangle) contains many small segments (triangles). At the edges
(corners) of the chain (triangle), there are two (three) isolated atoms for
$\kappa_{A}=0$, while there are dimer (trimer) states for $\kappa_{B}=0$. They
are marked by dotted circles. Lattice sites are numbered from $n=1$ to $N$ as
indicated in (a2) and (b2).

###  II.2 Edge or corner dynamics

We first consider the dynamics of an edge or corner site when it is perfectly
isolated as in Fig.1(a1) or (b1). This is the case where $M_{nm}=0$ for the
edge or corner sites. For instance, this is realized by setting $\kappa_{A}=0$
in Eq.(12) or $\lambda=-1$ in Eq.(15) for the SSH model.

The dynamics is governed by isolated equations,

$\frac{d\psi_{\bar{n}}}{dt}=-\gamma\left(1-\xi\frac{1}{1+\left|\psi_{\bar{n}}\right|^{2}/\eta}\right)\psi_{\bar{n}}.$
(3)

Solving this equation numerically with the initial conditions
$\psi_{\bar{n}}=1$ and $\psi_{\bar{n}}=0.1$, we show the results in Fig.2(a)
and (b), respectively. It is intriguing that the dynamics of the topological
edge or corner state does not depend on the initial condition of
$\psi_{\bar{n}}$. The saturated value of $\psi_{\bar{n}}$ as
$t\rightarrow\infty$ is identical for all isolated edge or corner states. This
can be understood analytically as follows.

We solve Eq.(3) for nontrivial stationary solutions. The stationary solution
for an edge or corner site $\bar{n}$ is given by

$\xi\frac{1}{1+\left|\psi_{\bar{n}}\right|^{2}/\eta}=1.$ (4)

Hence, the nontrivial solution reads

$\lim_{t\rightarrow\infty}\left|\psi_{{}_{\bar{n}}}\right|^{2}=\eta(\xi-1),$
(5)

where it is necessary that $\xi>1$. We have the trivial solution,
$\lim_{t\rightarrow\infty}\left|\psi_{{}_{\bar{n}}}\right|^{2}=0$, for
$\xi\leq 1$.

As a result, there are only two stable solutions in Eq.(3). One is the ground
state mode $\psi_{{}_{\bar{n}}}=0$, and the other is the stimulated mode
$\left|\psi_{{}_{\bar{n}}}\right|^{2}=\eta\left(\xi-1\right)$ for $\xi>1$. The
initial condition determines which state is realized.

The nonlinear term is essential to have a stationary solution. Indeed, we
obtain a linear theory in the limit $\eta\rightarrow\infty$ , where the
amplitude $\left|\psi_{n}\right|^{2}$ diverges.

The state remains to be real for a real initial condition, and Eq.(3) is
simplified

$\frac{d\psi}{-\psi_{\bar{n}}+\xi\frac{1}{1+\psi_{\bar{n}}^{2}/\eta}\psi_{\bar{n}}}=\gamma
dt,$ (6)

which is solved as

$\frac{2\log\psi_{\bar{n}}-\xi\log\left[\gamma\left(\psi_{\bar{n}}^{2}/\eta+1-\xi\right)\right]}{2\left(\xi-1\right)}=\gamma\left(t+t_{0}\right).$
(7)

The state evolution $\psi_{\bar{n}}\left(t\right)$ is given by the inverse of
this equation.

Figure 2: Dynamics of the isolated equation (3) under the initial conditions
(a) $\psi_{1}(0)=1$ and (b) $\psi_{1}(0)=0.1$. The solutions reach the same
stationary solution irrespective of the initial conditions. Red color
indicates the solution with $\xi=2$, while cyan color indicates that with
$\xi=0$. We have set $\eta=1$ and $\gamma=1/2$.

###  II.3 Bulk dynamics

We next analyze the dynamics of the bulk site, where there is no nonlinear
non-Hermitian term because $P_{n}=0$. Then, the dynamics is governed by the
linear equation

$i\frac{d\psi_{n}}{dt}=\sum_{nm}M_{nm}\psi_{m}-i\gamma\psi_{n}.$ (8)

With the use of a solution of the linear equation

$i\frac{d\psi_{n}^{0}}{dt}=\sum_{nm}M_{nm}\psi_{m}^{0},$ (9)

the solution of Eq.(8) is written in the form

$\psi_{n}=e^{-\gamma t}\psi_{n}^{0}.$ (10)

It means that the amplitude exponentially decays as a function of time in the
bulk.

Figure 3: Saturated amplitude $|\psi_{1}|$ at (a1) the edge site in the SSH
model, and (b1) the corner site in the breathing Kagome model, where
$\gamma=1/2$ for magenta, $\gamma=1/4$ for cyan and $\gamma=1/8$ for orange
colors. (a2) The energy spectrum in the SSH model made of a finite chain with
topological edge states in red. (b2) The energy spectrum of the breathing
Kagome model in triangle geometry with topological corner states in red. The
horizontal axis is $\lambda$. We have set $\kappa=1$, $\eta=1$ and $\xi=2$.

###  II.4 Quench dynamics

Quench dynamics starting from one site is a good signal to detect whether the
system is topological or trivialQWalk . It is also applicable to various
nonlinear systemsTopoToda ; MechaRot ; NLPhoto ; Sine . Let us study a quench
dynamics starting from one site indexed by $m$,

$\psi_{n}\left(t\right)=\delta_{n,m}\quad\text{at}\quad t=0,$ (11)

where $m=1$ represents the left-edge site or the top-corner site as in
Fig.1(a2) and (b2).

Figure 4: Density plot of the time evolution of the amplitude $|\psi_{n}|$ as
a function of $n$ in (a1) the topological phase with $\lambda=-0.5$ and (b1)
the trivial phase with $\lambda=0.5$ in the SSH model. The amplitude
$|\psi_{n}|$ as a function of $t$ for various $n$ in (a2) the topological
phase and (b2) the trivial phase, where red curves indicate $n=1$, while cyan
curves indicate $n=N$ in (a2) and (b2). The saturated amplitude $|\psi_{n}|$
as a function of $n$ in (a3) the topological phase and (b3) the trivial phase
after enough time. We have set $\kappa=1$, $\eta=1$, $\gamma=1/2$ and $\xi=2$.
We take a sample with $N=20$.

##  III Nonlinear non-Hermitian SSH model

###  III.1 Model

A topological laser based on the SSH model has been discussedSchome ; Weimann
; Jean ; Ota18 ; Parto ; Zhao ; Malzard ; MalzardOpt . We consider the case
where the hopping matrix is governed by the SSH matrix in Eq.(1). The matrix
is explicitly given by

$\displaystyle M_{nm}$ $\displaystyle=$
$\displaystyle-\left(\kappa_{A}+\kappa_{B}\right)\delta_{n,m}+\kappa_{A}\left(\delta_{2n,2m-1}+\delta_{2m,2n-1}\right)$
(12) $\displaystyle+\kappa_{B}\left(\delta_{2n,2m+1}+\delta_{2m,2n+1}\right).$

The explicit equations for a finite chain with length $N$ are given by

$\displaystyle i\frac{d\psi_{2n-1}}{dt}$ $\displaystyle=$
$\displaystyle\kappa_{A}\left(\psi_{2n}-\psi_{2n-1}\right)+\kappa_{B}\left(\psi_{2n-2}-\psi_{2n-1}\right)$
(13)
$\displaystyle-i\gamma\left(1-\xi\frac{\delta_{n,1}}{1+\left|\psi_{2n-1}\right|^{2}/\eta}\right)\psi_{2n-1},$
$\displaystyle i\frac{d\psi_{2n}}{dt}$ $\displaystyle=$
$\displaystyle\kappa_{B}\left(\psi_{2n+1}-\psi_{2n}\right)+\kappa_{A}\left(\psi_{2n-1}-\psi_{2n}\right)$
(14)
$\displaystyle-i\gamma\left(1-\xi\frac{\delta_{n,N}}{1+\left|\psi_{2n}\right|^{2}/\eta}\right)\psi_{2n}.$

It is convenient to introduce the coupling strength $\kappa$ and the
dimerization parameter $\lambda$ by

$\kappa_{A}=\kappa\left(1+\lambda\right),\quad\kappa_{B}=\kappa\left(1-\lambda\right),$
(15)

with $|\lambda|\leq 1$. The isolated edge limit is realized at $\lambda=-1$ in
Fig.1(a1).

The present SSH model (1) has the same topological structure as in the
original SSH model by the following reasoning. First, the contribution of the
nonlinear gain term ($\propto\gamma\xi$) is negligible in the bulk since it
exists only at the edge site. Next, the loss term ($\propto\gamma$) only shift
the matrix as

$\overline{M}_{nm}\equiv M_{nm}-i\gamma\delta_{nm}.$ (16)

Eq.(1) is rewritten in the form of the linear model for the bulk
($n\neq\bar{n}$),

$i\frac{d\psi_{n}}{dt}=\sum_{nm}\overline{M}_{nm}\psi_{m}.$ (17)

Then, the topological properties are determined by the matrix
$\overline{M}_{nm}$, and they are identical to those in the original SSH
model, as explained in Appendix: See Eq.(32). We show the band structure for a
finite chain in Fig.3(a2). The system is topological for $\lambda<0$ with the
emergence of the topological edge states marked in red, while it is trivial
for $\lambda>0$.

Figure 5: Density plot of the time evolution of the amplitude $|\psi_{n}|$ as
a function of $n$ in (a1) the topological phase with $\lambda=-0.5$ and (b1)
the trivial phase with $\lambda=0.5$ in the SSH model. We start from the site
$n=7$. We have set $\kappa=1$, $\eta=1$, $\gamma=1/2$ and $\xi=2$. We take a
sample with $N=20$.

###  III.2 Quench dynamics

We numerically solve Eqs.(13) and (14) under the initial condition (11) by
taking $m=1$. We show the time evolution of $\left|\psi_{n}\right|$ in Fig.4.
The quench dynamics is significantly different between the topological and
trivial phases. In the topological phase, the amplitude $|\psi_{1}|$ at the
left edge is rapidly saturated as shown in the red curve in Fig.4(a2), while
there is a delay in the saturation of the amplitude $|\psi_{N}|$ at the right
edge as shown in the cyan curve in Fig.4(a2). The delay implies the
propagation of a wave along a chain, although the propagation of the wave is
invisible in Fig.4(a1), as is consistent with Eq.(10). It means that a wave
with a tiny amplitude transfers the information of the excitation to the right
edge and induces a stimulated emission. On the other hand, in the trivial
phase, the amplitude $|\psi_{1}|$ rapidly decreases as shown in Fig.4(a2).
Furthermore, there is no excitation $|\psi_{N}|$ at the right edge.

Fig.4(a3) shows a spatial profile of the amplitude $\left|\psi_{n}\right|$
after enough time, which is the DOS for a pair of topological edge states in
nonlinear non-Hermitian system.

We find that there is no reflection of the propagating wave by the right edge.
This is due to the loss term ($\propto\gamma$) in the bulk. It is highly
contrasted with the case of the Hermitian model. In addition, the amplitudes
$\left|\psi_{1}\right|$ and $\left|\psi_{N}\right|$ are always identical. It
is due to reflection symmetry $x\longleftrightarrow N-x$ in the right-hand
side of Eqs.(13) and (14) because $d\psi_{n}/dt=0$ for the stationary
solution. This is confirmed in Eq.(5) explicitly for the limit ($\lambda=-1$)
of the isolated edge states, and numerically for any value of $\lambda$.

We show the saturated amplitude $|\psi_{1}|$ as a function of the dimerization
$\lambda$ in Fig.3(a1). It is finite for the topological phase although it
deviates from 1 other than $\lambda=-1$ due to the hopping term. On the other
hand, it is almost zero for the trivial phase. These features correspond to
the emergence or the absence of the topological edge states as shown in
Fig.3(a2). Namely, the quench dynamics well signatures the topological phase
transition although there are nonlinear non-Hermitian terms.

We also study the dynamics under the initial condition (11) by taking the site
$m$ in the bulk. The result is shown in Fig.5 by choosing $m=7$. All sites
belonging to the two topological edge states are stimulated after a delay with
the intensity depending on the DOS. The timing of the stimulation is
determined by the distance from the initial site.

Figure 6: (a) Dynamics of the amplitude $|\psi_{n}|$ of the dimer in the SSH
model (18) and (19). (b) The dynamics of the trimer in the breathing Kagome
model (25) and (26). We have set $\kappa_{A}=1$, $\eta=1$, $\gamma=1/2$ and
$\xi=2$.

Figure 7: Time evolution of the amplitude $\left|\psi_{n}\right|$ of the
breathing Kagome model in (a) the topological phase with $\lambda=-0.5$, (b)
the metal phase with $\lambda=0.2$ and (c) the trivial phase with
$\lambda=0.5$. The color density indicates the amplitude $|\psi_{n}|$. (a5)
and (c5) give the profiles of the topological corner states and the trimer
states together with their DOS, respectively. We have set $\kappa=1$,
$\eta=1$, $\gamma=1/2$ and $\xi=2$. We take a triangle with $N=108$.

###  III.3 Dimer states

We study the trivial phase ($\lambda>1$). In particular, we may solve the
equations of motion analytically in the dimer limit ($\lambda=1$). In this
case, we obtain a closed set of equation for the sites $1$ and $2$,

$\displaystyle i\frac{d\psi_{1}}{dt}$
$\displaystyle=\kappa_{A}\left(\psi_{2}-\psi_{1}\right)-i\gamma\left(1-\xi\frac{1}{1+\left|\psi_{1}\right|^{2}/\eta}\right)\psi_{1},$
(18) $\displaystyle i\frac{d\psi_{2}}{dt}$
$\displaystyle=\kappa_{A}\left(\psi_{1}-\psi_{2}\right).$ (19)

The stationary solutions of (18) and (19) are either the trivial one

$\psi_{1}=\psi_{2}=0,$ (20)

or a nontrivial one

$\psi_{1}=\sqrt{\eta\left(\xi-1\right)},\qquad\psi_{2}=\text{const}.$ (21)

We note that the right-hand side of Eq.(19) is not necessary to be zero
because only the phase rotates if it is not zero. It is a dynamical problem
depending on the initial condition which stationary solution is actually
chosen.

We show numerical solution with the initial condition $\psi_{1}=1$ and
$\psi_{2}=0$ in Fig.6(a). The amplitudes $\left|\psi_{1}\right|$ and
$\left|\psi_{2}\right|$ exponentially decreases to zero with oscillations.
Hence, the stationary solutions are $\psi_{1}=\psi_{2}=0$. Furthermore, we
have found numerically that $\left|\psi_{1}\right|$ is almost zero in the
trivial phase as in Fig.3(a). As far as we have checked numerically, there is
no nontrivial dimer solution in the trivial phase of the SSH model.

##  IV Nonlinear non-Hermitian breathing Kagome model

###  IV.1 Quench dynamics

We proceed to investigate the system where the hopping matrix $M_{nm}$
describes the breathing Kagome lattice, whose lattice structure is illustrated
in Fig.1(b). The hopping matrix and the topological number are given by (33)
and (38) in Appendix. There are topological, trivial and metal phases in the
breathing Kagome modelEzawaKagome  as shown in Fig.3(b2).

By solving Eq.(1) under the initial condition (11) with the choice of $m=1$,
we show the time evolution of the amplitude $\left|\psi_{n}\right|$ in Fig.7
and in Fig.8, where $n=1$ denotes the top-corner site. Fig.7 displays a global
picture how the stimulated signal at site $n=1$ propagates all over the sites
as time evolves. Fig.8 displays a detailed evolution along the left side of
the lattice. We also show the saturated amplitude as a function of $\lambda$
in Fig.3(b1). The saturated amplitudes are identical at the three corner sites
due to trigonal symmetry of the breathing Kagome lattice, as found in Fig.7.

We comment that Fig.7(a5) shows a spatial profile of the saturated amplitude
$\left|\psi_{n}\right|$, which is the DOS for three topological corner states
in nonlinear non-Hermitian system.

These figures show typically different behaviors in the topological, trivial
and metal phases. First, in the topological phase, the amplitudes at the other
two corner sites increase after a delay to the saturated amplitude as shown in
Fig.8(a2). Second, in the metal phase, the amplitude at the corner site
rapidly decreases as in Fig.8(b2). These features are very much similar to
those in the the SSH model.

On the other hand, in the trivial phase, there is a significant difference
between the SSH model and breathing Kagome model. In the SSH model, the
amplitude rapidly decreases as in Fig.4(b2). However, this is not the case for
the breathing Kagome model. The amplitude at the top corner site suddenly
decreases to the saturated value, while those at the bottom corner sites
increase to the saturated value after a delay as in Fig.8(c2). Furthermore,
the saturated amplitude $\left|\psi_{1}\right|$ depends on the magnitude of
$\xi$ as shown in Fig.3(b1). We pursue the reason why the amplitudes do not
vanish in the trivial phase. We will see that it is due to the formation of
the trimer state.

Figure 8: Density plot of the time evolution of the amplitude $|\psi_{n}|$
along one side of the breathing Kagome lattice in (a1) topological phase with
$\lambda=-0.5$, (b1) metal phase with $\lambda=0.2$, and (c1) trivial phase
with $\lambda=0.5$. Time evolution of the amplitude $\left|\psi_{n}\right|$ of
the breathing Kagome model for various $n$ in (a2) topological phase with
$\lambda=-0.5$, (b2) metal phase with $\lambda=0.2$, and (c2) trivial phase
with $\lambda=0.5$. We have set $\kappa=1$, $\eta=1$, $\gamma=1/2$ and
$\xi=2$. We take a triangle with $N=108$.

###  IV.2 Trimer states

We consider the limit $\lambda=1$, where the system is decomposed into a set
of trimers as in Fig.1(b3). The equations of motion is explicitly given by

$\displaystyle i\frac{d\psi_{1}}{dt}$
$\displaystyle=\kappa_{A}(\psi_{2}+\psi_{3})-i\gamma\left(1-\xi\frac{1}{1+\left|\psi_{1}\right|^{2}/\eta}\right)\psi_{1},$
(22) $\displaystyle i\frac{d\psi_{2}}{dt}$
$\displaystyle=\kappa_{A}\left(\psi_{1}+\psi_{3}\right),$ (23) $\displaystyle
i\frac{d\psi_{3}}{dt}$
$\displaystyle=\kappa_{A}\left(\psi_{1}+\psi_{2}\right).$ (24)

Without loss of generality we may set $\psi_{2}=\psi_{3}$ and obtain

$\displaystyle i\frac{d\psi_{1}}{dt}$
$\displaystyle=2\kappa_{A}\psi_{2}-i\gamma\left(1-\xi\frac{1}{1+\left|\psi_{1}\right|^{2}/\eta}\right)\psi_{1},$
(25) $\displaystyle i\frac{d\psi_{2}}{dt}$
$\displaystyle=\kappa_{A}\left(\psi_{1}+\psi_{2}\right).$ (26)

The stationary solutions are either the trivial one

$\psi_{1}=\psi_{2}=0,$ (27)

or a nontrivial one

$\psi_{1}=\sqrt{\eta\left(\xi-1\right)},\qquad\psi_{2}=\text{const},$ (28)

which is identical to the stationary solution (21) in the SSH model. There is
no $\gamma$ dependence in the stationary solution (28), which agrees with in
Fig.3(b1).

We have found analytically the trimer state in the limit of $\lambda=1$ in the
trivial phase. A numerical solution is given as a function of time $t$ in
Fig.6(b), where $\left|\psi_{1}\right|$ and $\left|\psi_{2}\right|$ approach
two saturated values (28), as is consistent with Fig.3(b1). Fig.3(b1) suggests
that the trimer state is formed also away from the limit $\lambda=1$ depending
on the value of $\gamma$. Indeed, Fig.7(c5) shows a spatial profile of three
trimer states at $\lambda=0.5$. The formation of trimer states implies the
presence of a nontopological laser in the breathing Kagome model. These
behaviors are contrasted with the dimer state in the SSH model, where
$\left|\psi_{1}\right|=0$ in the trivial phase as in Fig.3(a1).

##  V Discussion

A topological laser provides a unique arena of topology, non-Hermicity and
nonlinearity. We have studied topological lasing in the SSH lattice and the
breathing Kagome lattice. They have topological edge or corner states in the
topological phase. When any one site is stimulated, all sites belonging to the
topological edge or corner states begin to emit stable laser light depending
on the DOS. The results would be universal for the physics of topological
lasing with the use of topological edge or corner states.

Non-Hermicity and nonlinearity play essential roles for the stabilization of
lasing. Without the nonlinearity, there is no stable lasing because the
amplitude exponentially grows or decays. The amplifier plays an essential role
of the lasing. If there is no loss term, there should be reflection wave at
the right edge site or the bottom corner sites, which are absent in a
topological laser.

On the other hand, nontopological lasing is not universal. Indeed, the quench
dynamics is different between the SSH model and the breathing Kagome model in
the trivial phase. There is no stable laser emission in the SSH model.
However, once one site is stimulated, all three corner states emit stable
laser lights in the breathing Kagome model due to the formation of trimer
states at the three corners.

In conclusion, we have found stable laser emission to occur in the topological
phase. Our results show that topological lasers provide an ideal play ground
of nonlinear non-Hermitian topological physics, where the profile of a
topological edge or corner state is observable by measuring the intensity of
lasing.

The author is very much grateful to N. Nagaosa for helpful discussions on the
subject. This work is supported by the Grants-in-Aid for Scientific Research
from MEXT KAKENHI (Grants No. JP17K05490 and No. JP18H03676). This work is
also supported by CREST, JST (JPMJCR16F1 and JPMJCR20T2).

##  Appendix A Topological number in the SSH model

The hopping matrix of the SSH model (12) is given by

$\overline{M}\left(k\right)=-\left(\kappa_{A}+\kappa_{B}+i\gamma\right)I_{2}+M_{0}\left(k\right)$
(29)

in the momentum space, with

$q\left(k\right)=\kappa_{A}+\kappa_{B}e^{-ik},$ (30)

and

$M_{0}\left(k\right)=\left(\begin{array}[]{cc}0&q\left(k\right)\\\
q^{\ast}\left(k\right)&0\end{array}\right).$ (31)

The topological number in the original SSH model is given by the Berry phase

$\Gamma=\frac{1}{2\pi}\int_{0}^{2\pi}A\left(k\right)dk,$ (32)

where
$A\left(k\right)=-i\left\langle\phi(k)\right|\partial_{k}\left|\phi(k)\right\rangle$
is the Berry connection with $\phi(k)$ the eigenfunction of
$\overline{M}\left(k\right)$. Note that the diagonal term in Eq.(17) with
Eq.(12) does not contribute to the topological charge because the wave
function $\phi(k)$ does not depend on the diagonal term. Hence, the present
model (1) has the same phases as the original SSH model. The system is
topological for $\lambda<0$ and trivial for $\lambda>0$.

##  Appendix B Breathing Kagome model

The hopping matrix of the breathing Kagome model is given byEzawaKagome

$M\left(\mathbf{k}\right)=-\left(\begin{array}[]{ccc}0&h_{12}&h_{13}\\\
h_{12}^{\ast}&0&h_{23}\\\ h_{13}^{\ast}&h_{23}^{\ast}&0\end{array}\right)$
(33)

with

$\displaystyle h_{12}$
$\displaystyle=\kappa_{A}e^{i\left(k_{x}/2+\sqrt{3}k_{y}/2\right)}+\kappa_{B}e^{-i\left(k_{x}/2+\sqrt{3}k_{y}/2\right)},$
(34) $\displaystyle h_{23}$
$\displaystyle=\kappa_{A}e^{i\left(k_{x}/2-\sqrt{3}k_{y}/2\right)}+\kappa_{B}e^{i\left(-k_{x}/2+\sqrt{3}k_{y}/2\right)},$
(35) $\displaystyle h_{13}$
$\displaystyle=\kappa_{A}e^{ik_{x}}+\kappa_{B}e^{-ik_{x}}$ (36)

in the momentum space, where we have introduced two hopping parameters
$\kappa_{A}$ and $\kappa_{B}$ corresponding to upward and downward triangles
in Fig.1(b).

##  Appendix C Topological number in the breathing Kagome model

The topological number in the breathing Kagome model is defined byEzawaKagome

$\Gamma=3\left(p_{x}^{2}+p_{y}^{2}\right),$ (37)

where

$p_{i}=\frac{1}{S}\int_{\text{BZ}}A_{i}d^{2}\mathbf{k},$ (38)

with
$A_{i}=-i\left\langle\phi(\mathbf{k})\right|\partial_{k_{i}}\left|\phi(\mathbf{k})\right\rangle$
being the Berry connection with $x_{i}=x,y$, and $S=8\pi^{2}/\sqrt{3}$ being
the area of the Brillouin zone; $\phi(\mathbf{k})$ the eigenfunction of
$M\left(\mathbf{k}\right)$. We obtain $\Gamma=0$ for $1\geq\lambda>1/3$, which
is the trivial phase with no topological corner states. On the other hand, we
obtain $\Gamma=1$ for $-1\leq\lambda<-1/3$, which is the topological phase
with the emergence of three topological corner states. Finally, $\Gamma$ is
not quantized for $-1/3<\lambda<1/3$, which is the metal phase.

## References

  * (1) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (2) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (3) A. B. Khanikaev, S. H. Mousavi, W.-K. Tse, M. Kargarian, A. H. MacDonald, G. Shvets, Nature Materials 12, 233 (2013).
  * (4) M. Hafezi, E. Demler, M. Lukin, J. Taylor, Nature Physics 7, 907 (2011).
  * (5) M. Hafezi, S. Mittal, J. Fan, A. Migdall, J. Taylor, Nature Photonics 7, 1001 (2013).
  * (6) L.H. Wu and X. Hu, Phys. Rev. Lett. 114, 223901 (2015).
  * (7) L. Lu. J. D. Joannopoulos and M. Soljacic, Nature Photonics 8, 821 (2014).
  * (8) T. Ozawa, H. M. Price, N. Goldman, O. Zilberberg and I. Carusotto Phys. Rev. A 93, 043827 (2016).
  * (9) D. Leykam and Y. D. Chong, Phys. Rev. Lett. 117, 143901 (2016).
  * (10) A. B. Khanikaev and G. Shvets, Nature Photonics 11, 763 (2017).
  * (11) X. Zhou, Y. Wang, D. Leykam and Y. D. Chong, New J. Phys. 19, 095002 (2017).
  * (12) P. St-Jean, V. Goblot, E. Galopin, A. Lemaitre, T. Ozawa, L. Le Gratiet, I. Sagnes, J. Bloch and A. Amo, Nature Photonics 11, 651 (2017).
  * (13) Y. Ota, R. Katsumi, K. Watanabe, S. Iwamoto and Y. Arakawa, Communications Physics 1, 86 (2018)
  * (14) T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg and L. Carusotto, Rev. Mod. Phys. 91, 015006 (2019).
  * (15) Y. Ota, F. Liu, R. Katsumi, K. Watanabe, K. Wakabayashi, Y. Arakawa and S. Iwamoto, Optica 6, 786 (2019).
  * (16) T. Ozawa and H. M. Price, Nature Reviews Physics 1, 349 (2019).
  * (17) A. E. Hassan, F. K. Kunst, A. Moritz, G. Andler, E. J. Bergholtz, M. Bourennane, Nature Photonics 13, 697 (2019).
  * (18) Y. Ota, K. Takata, T. Ozawa, A. Amo, Z. Jia, B. Kante, M. Notomi, Y. Arakawa, S.i Iwamoto, Nanophotonics 9, 547 (2020).
  * (19) M. Li, D. Zhirihin, D. Filonov, X. Ni, A. Slobozhanyuk, A. Alu and A. B. Khanikaev, Nature Photonics 14, 89 (2020).
  * (20) H. Yoshimi, T. Yamaguchi, Y. Ota, Y. Arakawa and S. Iwamoto, Optics Letters 45, 2648 (2020).
  * (21) M. Kim, Z. Jacob and J. Rho, Light: Science and Applications 9, 130 (2020).
  * (22) S. Iwamoto, Y. Ota and Y. Arakawa, Optical Materials Express 11, 319 (2021).
  * (23) E. Prodan and C. Prodan, Phys. Rev. Lett. 103, 248101 (2009).
  * (24) Z. Yang, F. Gao, X. Shi, X. Lin, Z. Gao, Y. Chong and B. Zhang, Phys. Rev. Lett. 114, 114301 (2015).
  * (25) P. Wang, L. Lu and K. Bertoldi, Phys. Rev. Lett. 115, 104302 (2015).
  * (26) M. Xiao, G. Ma, Z. Yang, P. Sheng, Z. Q. Zhang and C. T. Chan, Nat. Phys. 11, 240 (2015).
  * (27) C. He, X. Ni, H. Ge, X.-C. Sun,Y.-B. Chen1 M.-H. Lu, X.-P. Liu, L. Feng and Y.-F. Chen, Nature Physics 12, 1124 (2016).
  * (28) H. Abbaszadeh, A. Souslov, J. Paulose, H. Schomerus and V. Vitelli, Phys. Rev. Lett. 119, 195502 (2017).
  * (29) H. Xue, Y. Yang, F. Gao, Y. Chong and B.Zhang, Nature Materials 18, 108 (2019).
  * (30) X. Ni, M. Weiner, A. Alu and A. B. Khanikaev, Nature Materials 18, 113 (2019).
  * (31) M. Weiner, X. Ni, M. Li, A. Alu, A. B. Khanikaev, Science Advances 6, eaay4166 (2020).
  * (32) H. Xue, Y. Yang, G. Liu, F. Gao, Y. Chong and B. Zhang, Phys. Rev. Lett. 122, 244301 (2019).
  * (33) C. L. Kane and T. C. Lubensky, Nature Phys. 10, 39 (2014).
  * (34) B. Gin-ge Chen, N. Upadhyaya and V. Vitelli, PNAS 111, 13004 (2014).
  * (35) L. M. Nash, D. Kleckner, A. Read, V. Vitelli, A. M. Turner and W. T. M. Irvine, PNAS 112, 14495 (2015).
  * (36) J. Paulose, A. S. Meeussen and V. Vitelli, PNAS 112, 7639 (2015).
  * (37) R. Susstrunk, S. D. Huber, Science 349, 47 (2015).
  * (38) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, R. Thomale, Nat. Phys. 14, 925 (2018).
  * (39) C. H. Lee , S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling and R. Thomale, Communications Physics, 1, 39 (2018).
  * (40) T. Helbig, T. Hofmann, C. H. Lee, R. Thomale, S. Imhof, L. W. Molenkamp and T. Kiessling, Phys. Rev. B 99, 161114 (2019).
  * (41) Y. Lu, N. Jia, L. Su, C. Owens, G. Juzeliunas, D. I. Schuster and J. Simon, Phys. Rev. B 99, 020302 (2019).
  * (42) Y. Li, Y. Sun, W. Zhu, Z. Guo, J. Jiang, T. Kariyado, H. Chen and X. Hu, Nat. Com. 9, 4598 (2018).
  * (43) M. Ezawa, Phys. Rev. B 98, 201402(R) (2018).
  * (44) M. Ezawa, Phys. Rev. B 99, 201411(R) (2019).
  * (45) M. Ezawa, Phys. Rev. B 99, 121411(R) (2019).
  * (46) G. Harari, M. A. Bandres, Y. Lumer, M. C. Rechtsman, Y. D. Chong, M. Khajavikhan, D. N. Christodoulides, M. Segev, Science 359, eaar4003 (2018)
  * (47) M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren, Mordechai Segev, Demetrios N. Christodoulides, Mercedeh Khajavikhan Science 359, 1231 (2018)
  * (48) H. Schomerus, Opt. Lett. 38, 11 (2013)
  * (49) S. Weimann, M. Kremer, Y. Plotnik, Y. Lumer, S. Nolte, K. G. Makris, M. Segev, M. C. Rechtsman and A. Szameit, Nat. Mat. 16, 433 (2017)
  * (50) M. Parto, S. Wittek, H. Hodaei, G. Harari, M. A. Bandres, J. Ren, M. C. Rechtsman, M. Segev, D. N. Christodoulides and M. Khajavikhan, Phys. Rev. Lett. 120, 113901 (2018)
  * (51) H. Zhao, P. Miao, M. H. Teimourpour, S. Malzard, R. El-Ganainy, H. Schomerus and L. Feng, Nat. Com. 9, 981 (2018)
  * (52) S. Malzard and H. Schomerus, New. J. Phys. 20, 063044 (2018)
  * (53) S. Malzard, E. Cancellieri, and H. Schomerus Optics Express 26, 22506 (2018)
  * (54) H. Zhong, Y. V. Kartashov, A. Szameit, Y. Li, C. Liu and Y. Zhang, APL Photonics 6, 040802 (2021)
  * (55) F. Zhang, C.L. Kane and E.J. Mele, Phys. Rev. Lett. 110, 046404 (2013).
  * (56) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, 10.1126/science.aah6442.
  * (57) F. Schindler, A. Cook, M. G. Vergniory, and T. Neupert, in APS March Meeting (2017).
  * (58) Y. Peng, Y. Bao, and F. von Oppen, Phys. Rev. B 95, 235143 (2017).
  * (59) J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Phys. Rev. Lett. 119, 246401 (2017).
  * (60) Z. Song, Z. Fang, and C. Fang, Phys. Rev. Lett. 119, 246402 (2017).
  * (61) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Phys. Rev. B 96, 245115 (2017).
  * (62) F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig, and T. Neupert, Science Advances 4, eaat0346 (2018).
  * (63) M. Ezawa, Phys. Rev. Lett. 120, 026801 (2018).
  * (64) J. Wu, X. Huang, J. Lu, Y. Wu, W. Deng, F. Li, and Z. Liu, Phys. Rev. B 102, 104109 (2020).
  * (65) M. Ezawa, Phys. Rev. B 100, 165419 (2019).
  * (66) M. Ezawa, cond-mat/arXiv:2105.10851.
  * (67) M. Ezawa, J. Phys. Soc. Jpn. 90, 114605 (2021).
  * (68) M. Ezawa, arXiv:2110.06578.
  * (69) M. Ezawa, arXiv:2110.15602


# Transport Regimes of Underdamped Brownian Particles in a Tilted Washboard
Potential

Trey Jiron  Dept.of Physical and Environmental Sciences, Colorado Mesa
University
Grand Junction, CO.  Marygrace Prinster  Dept.of Mathematics and Statistics,
Colorado Mesa University
Grand Junction, CO.  Jarrod Schiffbauer<EMAIL_ADDRESS>Dept.of
Physical and Environmental Sciences, Colorado Mesa University
Grand Junction, CO.

###### Abstract

Abstract

In this paper, a comprehensive examination of the temperature- and bias-
dependent diffusion regimes of underdamped Brownian particles is presented. A
temperature threshold for a transition between anomalous and normal diffusive
behaviors is located, yielding a new phase diagram for the system. In the low-
temperature regime, the system exhibits an apparent negative differential
mobility due to persistent, long-time subdiffusion at low-bias; at high
temperature (or critical bias,) the system rapidly approaches normal diffusion
below an intermediate barrier height, $U_{0}\sim 6k_{B}T$. By consideration of
numerical results, comparison to the overdamped case, and the related Kramers
multistable escape problems, it is demonstrated that the low-bias non-
monotonic temperature dependence of the diffusivity, persistent subdiffusion,
and negative differential mobility can be traced to inertial effects, which
are evident in the oscillatory modes of the velocity power spectra at low
bias. In the giant diffusion regime, the velocity power spectra exhibit
coupling between the “locked” and “running” states, with a characteristic
frequency corresponding to the principal frequency of the limit cycles of a
damped, driven plane pendulum near critical bias. Non-linear second harmonic
generation, corresponding to oscillatory transient anomalous diffusivity, is
observed with increasing bias and decreasing temperature, further emphasizing
that the low-noise diffusion problem converges to noise-free dynamics,
complementing analytic results for the average velocity [L. Cheng and N.K.
Yip, Physica D, 2015].

††preprint: APS/123-QED

##  I I. Introduction and Background

The problem of Brownian motion in a superposition of a periodic landscape with
a bias force, the so-called “tilted washboard” potential, has been used to
study a wide range of mesoscale non-equilibrium phenomena. It describes driven
diffusive transport in diverse physical systems including colloids [1, 2, 3],
super-ionic conductors [4], biomolecular motors [5], diffusion of atoms on
solid surfaces [6, 7], and solid-state Josephson junctions [8]. A range of
results has been obtained demonstrating a rich phenomenology in a relatively
simple system, such as diffusion enhancement [9, 10], non-monotonic
temperature dependence of the diffusivity [11, 12, 13], and in some cases,
regions of negative differential mobility, i.e. a local minimum in the
diffusivity with respect to the bias force [14]. Such effects are of more than
purely fundamental interest, having applications in fields such as
biomolecular sensing and separations [15], colloidal assembly and active
materials [16, 17, 18, 19, 20].
Starting with Stratonovich [21], the existence of tilt-dependent regimes of
“running” vs. “locked” solutions suggested that such systems may be described
by a two-state Markov process. Extensive theoretical studies have since been
carried out of the Fokker-Planck equation governing diffusion of probability
for the whole ensemble [22, 23, 24]. From the point of view of Langevin
equations for the over-damped system, Reimann et al [9, 10] demonstrated giant
diffusion enhancement is a universal effect, and furthermore that the enhanced
diffusion peak becomes infinitely high and narrow in the limit $T\rightarrow
0$ for the overdamped case. Using a modified Lifson-Jackson formula for the
overdamped case, Berezhkovskii and Dagdug [14], found that the shape of the
potential controls the effective diffusion, resulting in a single maximum (for
sinusoidal potentials), a local minimum prior to an asymptotic approach to
saturated enhancement (for “cusped” potentials,) and a monotonic approach to
saturation (for square-well potentials.) Lindner and Sokolov [11] studied the
temperature dependence in the underdamped case and demonstrated a non-
monotonic temperature-dependence of the enhanced diffusion for a finite range
of bias forces. Similar effects are seen for time-periodic forcing in spatial
potentials [12]. This is explored in detail in a recent paper by Spiechowicz
and Luczka [13] in the underdamped case, where it is argued that the non-
monotonic temperature dependence arises due to a temperature-driven change
evident in the long-time stationary velocity probability distribution.
While Brownian motion in a tilted periodic potential is well-studied, there is
little unifying or clarifying discussion of the dynamics in these different
regimes. Here, we aim to remedy this by recourse to detailed numerical
studies, enabling their analysis using the power-spectra of the velocity to
gain more insight into the dynamics. By comparison to the zero-noise power
spectrum, the results demonstrate numerically the convergence of the full,
time-dependent dynamics towards those of the noise-free case. This complements
an important analytical result; the long-time average velocities and
transition times in the underdamped tilted periodic potential have been shown
to converge to those of the noise-free system in the limit of low noise by
Cheng and Yip [25] at moderate to high bias, i.e. in the region near the
giant-diffusion effect.
The remainder of the paper is structured as follows: In Section II the chosen
normalization and numerical implementation of the Langevin equation are
explained. In Section III, initial simulation results and qualitative
observations are introduced for particle diffusivity over a range of bias
force along with a new phase diagram in well-depth vs. temperature. To connect
the range of diffusive behaviors to the underlying dynamics, the zero-bias
dynamics, jump dynamics, and zero-noise system are discussed separately before
considering the role of bias across temperature on the full dynamics of the
system. Section IV then reiterates and discusses the main conclusions, along
with possible applications and future work.

##  II II. Model and simulation

We employ a standard Langevin equation for a single Brownian particle of mass
$m$ immersed in a bath at uniform temperature $T$, with the particle subject
to a biased cosine potential, e.g. the “tilted washboard” potential
$U(x)=U_{o}\cos{(\kappa x)}+F_{o}x$,

$m\ddot{x}+\gamma\dot{x}-\kappa U_{o}\sin{\kappa x}+F_{o}-\sqrt{2\gamma
k_{B}T}\xi(t)=0$ (1)

where $\gamma$ is the viscous friction of the bath and the noise term obeys
$\langle\xi(t)\rangle=0$ and $\langle\xi(t)\xi(0)\rangle=\hat{\delta}(t)$. The
length is scaled by the inverse wavenumber $\kappa^{-1}$ and the time by the
viscous relaxation time, $m/\gamma$, so that the governing equation becomes

$\ddot{y}+\dot{y}-\alpha\sin{y}+\beta-\delta\eta(t)=0$ (2)

with a dimensionless Gaussian noise $\eta(t)$.
The three dimensionless parameters are the coupling to the bath,
$\delta=\kappa\sqrt{2mk_{B}T}/\gamma$; the external bias, $\beta=m\kappa
F_{o}/\gamma^{2}$; and the potential landscape,
$\alpha=m\kappa^{2}U_{o}/\gamma^{2}$. It will be helpful to the reader to note
that $\delta\sim{T^{1/2}}$. These three dimensionless numbers represent ratios
of the characteristic energies of respective processes in the system to
viscous dissipation. Here, for a temperature-independent friction, the scaling
places us in the intermediate-high-damping regime of the corresponding Kramers
problem.
Equation 1 is reduced to a pair of coupled, first-order equations by defining
$\dot{y}=v$ and explicit trajectories $y(t)$ are obtained numerically using
RK4 integration implemented in Python. Due to the possible long run times of
each trajectory and the large number of trajectories required for any
meaningful statistical information to be extracted, efficient calculation
requires parallel processing on an Intel Xeon “Haswell” cluster. The effective
diffusivity is then computed directly from the variance of the trajectories,

$D=\lim_{t\to\infty}\frac{\langle y^{2}(t)\rangle-\langle
y(t)\rangle^{2}}{2t}$ (3)

with ensembles of $N=1000-2000$ trajectories run for $10^{6}$ to $10^{7}$
time-steps and total system times of $10^{4}$ to $10^{5}$ (dimensionless time,
normalized by the viscous relaxation time) averaging over every $10^{th}$ step
of the last $10^{4}$ points in each trajectory. A well-defined constant
diffusivity does not exist for all parameter ranges; there are regimes of
persistent, anomalous diffusion [13]. Thus, it is convenient to introduce the
time-dependent diffusion coefficient,

$D(t)=\frac{\langle y(t)\rangle^{2}-\langle y^{2}(t)\rangle}{2t}$ (4)

which will facilitate our discussion in super- and sub-diffusive regimes.

##  III III. Results and Discussion

Figure 1: Effective diffusivity vs. bias parameter, $\beta$, normalized by
Einstein-Sutherland value. Top (a) figure shows curves for $\alpha=1$ and a
range of $\delta$. Inset shows un-normalized diffusivity in neighborhood of
giant diffusion maximum. Note that negative differential mobility is only
evident in the log scaled y-axis. Bottom (b) figure shows curves at
$\delta=0.25$ for several values of $\alpha$.

The effective diffusivity vs. bias, $\beta$ is shown in Figs. 1. The top
figure shows the diffusivity for $\alpha=1.0$ and several values of $\delta$,
as indicated. Let us make a few qualitative observations regarding the
response of the diffusion coefficient. The well-known giant diffusion effect
[10] and non-monotonic temperature response of the giant diffusion [11]
(emphasized in the inset) are evident. Another curious response is also
evident, as emphasized by the vertical log scale; below the onset of giant
diffusion, the effective diffusivity decreases in response to increasing bias,
i.e. corresponding to a negative differential mobility, and then increases
again prior to the onset of giant diffusion. While the existence of “locked”
trajectories is well-known [21], this effective decrease in diffusivity with
applied bias does not result in a complete or immediate vanishing of the
diffusivity. Rather, at least at sufficiently low $\delta$ ($\sim T^{1/2}$),
there is a slow approach to a minimum diffusivity, before rising again with
the onset of diffusion enhancement. A non-monotonic minimum for cusped
potentials has been reported in the overdamped case [14], the phenomenon has
not been previously reported for the cosine potential in either underdamped or
overdamped regimes. This effect has a pronounced temperature dependence, in
particular, disappearing above some approximate temperature threshold. Similar
phenomena occur in the transition from static-to-kinetic friction, grain
boundary hysteresis in shape memory alloys and other stick-slip phenomena
[26]. At high bias, the diffusivity is dominated by the frictional response,
with a weak temperature dependence (for temperature-independent friction.)
What is noteworthy is that the width of the transition from a region of
apparent minimum effective diffusivity through the giant diffusion maximum
displays similar behavior for decreasing $\delta$ as it does for increasing
$\alpha$. As $\delta$ decreases or $\alpha$ increases, there is a widening and
lowering of the apparent “locked” plateau of low but non-zero diffusivity, as
well as a “sharpening” and (finite) increase of the diffusion maximum peak.
To better characterize the nature of the diffusion for various parameter
ranges, the exponent, $q$, in the mean-square deviation at long times, e.g.
$MSD\sim t^{q}$ is calculated. As usual, $q\sim 1$ denotes normal diffusion,
whereas $q<1$ denotes sub-diffusive and $q>1$ super-diffusive relaxation (with
the ballistic limit, $q=2$.) To do this, we sample the MSD at 4 time-intervals
of increasing length, i.e. the last $n\Delta t=1000n$ times up to $t=10,000$,
and display the average of fitted exponents, along with error bars indicating
standard deviation in the exponents across time intervals, for a range of
$\delta,\beta$ values in Fig. 2.

Figure 2: Fitted exponents from mean square deviation of trajectories at long
times. The black states lie below the transition $\delta_{c}\sim 0.57$
(denoted by the red markers) for $\alpha=1$ and a range of bias; blue states
lie above.

Several interesting observations can be made from Fig. 2. First, overall, the
states above the rough threshold range of $\delta\sim 0.5-0.57$ for $\alpha=1$
generally reach normal diffusion by the last interval of $\Delta t=1000$ for
all bias values, while those below only do so for biases corresponding to
diffusive amplification (and higher). Second, for non-zero bias across all
parameter ranges below the onset of diffusive amplification, i.e.
$0<\beta<\beta_{c}$ for some $\beta_{c}$, there is some degree of persistent,
long-time sub-diffusion due to the bias itself, corresponding to the apparent
negative differential mobility. As the bias approaches a critical value, here
$\beta\sim 1$, the system rapidly approaches normal diffusion for $\delta$ in
the range $0.125$ to $1.0$, with a marked reduction in standard deviation
about $q=1$ as compared to sub-critical forcing.
The above results demonstrate that the overall behavior of diffusion
coefficient as a function of $\beta$ changes as $\delta$ increases.
Corresponding to the transition from a largely subdiffusive regime to one of
essentially normal diffusion, we note that the slope of the diffusion
coefficient with respect to $\beta$ changes from negative to positive around
$\delta\sim 0.5$. By examining the decrease in diffusivity with temperature
($\delta$) for fixed well-depth ($\alpha$), as exhibited in Fig. 1a and b, it
is clear that the apparent negative differential mobility occurs when the well
depth becomes sufficiently large compared to the thermal energy. Normalizing
the diffusion coefficient by the zero-bias value (see inset of Fig. 3), the
transition occurs about the curve shown in Fig. 3, which can be anticipated to
go as $\delta^{2}$.

Figure 3: Phase diagram in $\alpha-\delta$ plane. The dashed curves correspond
to $m=2.25$ (black) and $m=1.87$ (blue); the solid black curve is the fit for
$m=2$. The region above the curve(s) corresponds to the intermediate-high
damping regime in terms of the Kramers turnover. Inset shows diffusivity
normalized by zero-bias value used to obtain bracketing values of $\delta$ for
$\alpha=0.25,0.5$. Brackets for the points at $\alpha=1.0,1.5,2.0$ were
similarly obtained.

We verify our intuition by fitting the high and low points bracketing the
transition to a power-law curve,

$\alpha=\lambda\delta^{m}$ (5)

where the exponent $m$ is found to lie between $1.87$ and $2.25$
(corresponding to the dashed curves in Fig. 3.) Proceeding to assume
$\alpha=\lambda\delta^{2}+\alpha_{o}$, we find $\lambda=3.209$ for the middle
curve, which indicates that, roughly speaking, $U_{o}\sim 6k_{B}T$. Note that
this temperature-dependent transition in diffusivity is evident in the zero-
bias diffusivity of Fig. 1a. As $\delta$ increases, the lowest diffusivity in
the plot occurs near $\delta=0.5$, with the value rising back up as the
temperature increases.

###  III.1 Zero-bias dynamics

In Fig. 1a, it is apparent that even the zero-bias diffusivity responds non-
monotonically to changes in temperature ($\sim\delta^{2}$). The fact that low,
but non-zero bias exhibits a decrease in the diffusivity, and that the range
of the non-monotonic temperature response extends from $0$ to above the
critical bias, suggests that a closer examination of the zero bias response is
warranted. Note both underdamped and overdamped systems exhibit a distinct
change in response at zero bias with respect to $\delta$. The overdamped
response can be seen in Fig. 4, where the Lifson-Jackson formula [27]

$\frac{D_{LJ}}{D_{o}}=\frac{(2\pi)^{2}}{\int_{-\pi}^{\pi}dy\exp{\left(-\frac{\alpha}{\delta^{2}}\sin{y}\right)}\int_{-\pi}^{\pi}dy\exp{\left(\frac{\alpha}{\delta^{2}}\sin{y}\right)}}$
(6)

is used to plot the zero-bias normalized diffusivity vs. $\delta$ for
$\alpha=0.25,0.5,1.0,2.0$. The diffusivity varies very slowly until
approaching $\alpha\sim 3\delta^{2}$, where the diffusivity undergoes a rapid,
however monotonic, increase with temperature. Considering the un-biased
periodic potential, Ferrando et al [28] used a numerical solution of the
corresponding Klein-Kramers equation, along with analytical approximations of
the Melnikov type [29], to demonstrate that the jump rates and jump length
distributions had strong dependence on temperature. While this temperature-
dependence (relative to barrier height) is present across a wide range of
friction, the rates and probabilities for multiple jumps exhibit a distinct
maximum in the underdamped regime. It is unsurprising that both overdamped and
underdamped dynamics exhibit a transition around the same temperature.
However, note that the overdamped case shows monotonic temperature response,
in contrast to the zero-bias values shown in Fig. 1 in which $\delta\sim
0.5-0.57$ shows a minimum in zero-bias diffusivity at long time. Thus, this
non-monotonic response at low bias is tied intimately to inertial effects.

Figure 4: Lifson-Jackson formula diffusivity (Eqn. 6) vs $\delta$ for several
values of $\alpha$. Rapid increase of $D_{LJ}/D_{o}$ shows critical $\delta$
behavior; but the overall response is monotonic

The time-dependent (zero-bias) diffusion coefficient for the underdamped case
with cosine potential is plotted in Fig. 5 for several values of $\delta$. The
non-monotonicity with respect to temperature is persistent at long times,
consistent with the subdiffusion, and apparent negative differential mobility.
For lower values of $\delta$, the diffusivity decreases with longer simulation
times as the “locked” state waiting time increases, leading to persistent
subdiffusion as the noise decreases. For example, at $\tau=100,000$, the
diffusivity is about an order of magnitude smaller than at $\tau=10,000$,
suggesting that the diffusivity decreases roughly as $1/\tau$. However, as
$\delta$ increases, the diffusivity settles much more rapidly (see Fig. 5).
Note the contrast between the non-monotonic temperature dependence exhibited
by the underdamped system, as opposed to the monotonic temperature response of
the overdamped diffusivity exhibited in Fig. 4.

Figure 5: Zero-bias normalized diffusion coefficient vs. time for several
temperatures (indicated) in the $\alpha=1$ case. Inset shows same figure at
very long times. Lowest $\delta$ values exhibit persistent sub-diffusion at
long times.

The negative differential mobility occurs in a distinct region of the
parameter space for low $\beta$ and $\delta$, corresponding to a persistent
long-time sub-diffusive regime. The same temperature threshold occurs in the
overdamped system, but there the zero bias diffusivity increases monotonically
with temperature. Generally speaking (though not visible on the resolution of
the figures), the normalized overdamped diffusivity approaches zero far faster
with decreasing $\delta$ than does the zero-bias diffusivity of the
underdamped system _below the temperature threshold,_ while at higher
temperatures, the overdamped system tends to have a higher normalized
diffusivity at zero bias. A deeper understanding of the situation may be
obtained by considering the jump dynamics explicitly.

###  III.2 Jump dynamics

Noise-assisted escape from a periodic well involves inherently multistable
dynamics[28], and the symmetry-breaking effect of the bias introduces another
complication. Our assertion is that the decrease in diffusivity with
increasing bias and corresponding subdiffusion can be attributed to asymmetry
in the forward and reverse-biased jump rates and a concomitant decrease in the
overall jump rate at low bias. To test this assumption, the forward and
reverse-biased jump rates are calculated, averaging over the well indices for
forward and reverse jumps separately.
Calculating directional well jumps from particle trajectories, naively one
expects the biased system to exhibit a lower reverse-biased jump rate compared
to the forward-biased rate. One also expects that the forward biased jump rate
will increase as the reverse jump rate decreases. The reader will note that in
Fig. 6b, at low $\beta$ and low $\delta$, this expectation is borne out.
However, it can be seen there is a maximum in the forward-biased jump rate at
low $\beta$, and the corresponding value of bias for this maximum rate
increases with temperature. Beyond this maximum, the forward-biased jump rate
decreases to a local minimum, and corresponding with the onset of the giant
diffusion effect, this rate increases dramatically. Thus, at $\delta$ below
the identified transition, the forward-biased jump rate possesses both a local
maximum (at very low bias,) and a local minimum immediately prior to the onset
of giant diffusion. This local maximum and minimum move closer together with
increasing temperature and annihilate for $\delta$ in the vicinity of $0.57$
such that the forward biased jump rate increases monotonically at higher
$\delta$.
The reverse-biased jump rate is no less intriguing. It is worth pointing out
that, above critical bias, the forward jump rate is a universal function of
bias and does not depend on temperature, while the temperature-dependence of
the reverse-biased jump rate is unexpected. At low temperatures, the reverse-
bias rate decreases monotonically with increasing bias, as one would expect.
However, nearing the transition temperature, the reverse jump rate develops a
local maximum _in the vicinity of the giant diffusion maximum, approaching
critical $\beta$_, exceeding the zero bias reverse jump rate. This maximum
increases with respect to both the global reverse jump rate and the zero bias
jump rate; but critically it exceeds the zero bias value by almost an order of
magnitude around $\delta=0.75$, coincident with the maximum un-normalized
giant diffusion (see inset of Fig. 1a.) Above this, the maximum becomes less
pronounced as the temperature continues to increase. As will be explained in
the following section, this corresponds to bias-driven coupling between
“locked” and “running” states across the I-III bifurcation, which broadens
with respect to bias as the temperature increases.
Keeping in mind that $\beta$ represents the characteristic work done by the
external force as the particle is forced through one period, increasing
$\beta$ increases the kinetic energy, bringing the lowest energy states closer
to the top of the barrier and increasing the kinetic energy of those already
at or above the barrier top. Eventually, at sufficiently large $\beta$ for a
given $\delta$, this results in an increase in particle current more than
sufficient to offset any reduction of jump rate due to asymmetry and the
diffusivity begins to rise again with increasing $\beta$

Figure 6: Mean directional jumps averaged over the full ensemble for
$\alpha=1$ (a) Reverse-biased (-) and forward-biased (+) jumps for a range of
$\delta$ and (b) forward jumps at low $\beta$ showing the low-bias maximum at
the lowest $\delta$. Since all simulations are performed for the same maximum
time, $t_{f}=10,000$, these may be interpreted as jump rates without loss of
generality (multiplying by the appropriate $10^{-4}$ factor.)

###  III.3 Zero-noise dynamics

The peculiar features of the jump dynamics, particularly the decrease in
forward-biased jump rate and existence of strong reverse-biased motion, can be
better understood with recourse to the behavior of the system in the absence
of thermal fluctuations,

$\ddot{y}+\dot{y}-\alpha\sin{y}+\beta=0.$ (7)

The model, equivalent to the noise-free Josephson junction or the damped,
driven pendulum, has a rich oscillatory dynamics typical of inertial systems
with bifurcations between three regions in Fig. 7 (see for example [30] and
references therein.) Region I contains both stable sinks and unstable saddles,
corresponding to the well bottoms and tops. Region III contains limit cycle
oscillations, and the Jacobian corresponding to the first-order system
represented by Eqn. 7 admits stable spirals as well as saddle nodes in Region
II. The role of thermal fluctuations is particularly important around the
bifurcations. Here, for $\beta$ and $\alpha$ sufficiently close, fluctuations
can couple trajectories on either side of the line. This results in the giant
diffusion effect itself by maximizing the spread in trajectories which
initially started near the same phase point. For lower $\delta$ values,
thermal fluctuations make the region III limit cycles available to phase
points starting in Region I and the diffusion maximum occurs just below the
$\alpha=\beta$ line. For the remainder of the paper, we confine our attention
to the dynamics approaching and exceeding the I-III bifurcation. As the long-
time diffusivity for $\alpha=1.5$ (see Fig. 1b) has similar overall behavior
(at least for lower $\delta$,) we leave investigation of the dynamics at the
I-II and II-III bifurcations for later study.

Figure 7: Bifurcation diagram for the noise-free dynamical system, Eqn. 7.
This is essentially the same as Fig.(8.5.10) Strogatz [30], but in the scaling
used here. The solid blue line indicating the homoclinic bifurcation diverges
from the infinite period bifurcation at $\alpha=5/4$.

###  III.4 Bias- and temperature- dependent dynamics

To facilitate discussion of the role of bias by providing some insight into
the time-dependent dynamics of the stochastic system, the velocity power
spectrum is introduced:

$a(\omega)=\left|\int_{0}^{\infty}dt\exp{(-i\omega t)}v(t)\right|^{2}$ (8)

The power spectra for a range of $\delta$ and several bias values below the
giant diffusion effect are shown in Fig. 8, with solid lines indicating zero
bias and biased states indicated by interrupted lines, as explained in the
caption. The zero and lowest bias states exhibit a distinct maximum, with the
maximum approaching the frequency of the oscillations about the well minimum
as $\delta$ (temperature) increases. Note that the overall frequency
distribution shifts energy to lower frequencies as the bias increases. This
can be seen in the corresponding linearized damped oscillator problem, where
the frequency of oscillations about the (stretched) minimum is given by the
(imaginary) root $\sqrt{1-4\alpha\sqrt{1-\beta^{2}/\alpha^{2}}}/2$. Thus,
below critical bias, energy shifts to lower frequencies with increasing bias.
Since $2/\delta$ gives the typical time-scale for thermal transition between
wells, velocity oscillations corresponding to frequencies less than the
inverse value indicate well-to-well, i.e. sink-over-saddle, transitions. Thus,
as $\delta$ increases, there is more energy being made available to the low-
frequency jumps, while simultaneously the thermal velocity is increasing.
Still, for $\delta$ below about $0.5-0.57$ (for $\alpha=1$,) this is
negligible in comparison to the energy being stored in the vibrational modes
inside the wells. An important observation is that the simple-harmonic
oscillatory modes in the “locked” states arise due to the underlying
oscillatory dynamics enabled by the inertial term in the zero-noise system. As
$\delta$ increases above the threshold and $\beta$ increases towards $\alpha$,
the power spectrum increasingly resembles the $1/f$ “pink noise” spectrum, a
characteristic of noise-driven transitions [31].

Figure 8: Power spectra for several cases for $\alpha=1$ at zero bias (solid
lines) and biased cases indicated by broken lines as follows: $\beta=0.25$
(dashed), $\beta=0.5$ (dot-dashed), and $\beta=0.75$ (dotted).

Figure 9: Time-dependent diffusion coefficients for several cases for
$\alpha=1$ at early times. Long-time behavior is reflected in Fig.1. Top (a)
shows decrease of initial diffusion coefficient with increasing bias for
$\delta=0.25$, typical of temperatures below the critical. Figs (b) and (c)
show the transition-region temperature ($\delta=0.57$)and the high-temperature
($\delta=1.0$) behavior, respectively.

A comparison of Fig. 5 and Fig. 9 shows that the increase in tilt and increase
in temperature can both result in reduced peak anomalous diffusivity and
faster subdiffusive decay to normal diffusion. The peak diffusion coefficient
at the end of the ballistic-dynamical transient ($t\sim 10$) decreases with
increasing bias for all cases. But the effect is more drastic at lower
temperatures.
The diffusion coefficients and power spectra as $\beta$ approaches $\alpha$ in
an infinite-period (I-III) bifurcation are shown for several values of
$\delta$ in Fig. 10. As $\delta$ increases, the power spectrum approaches
$1/f$, as mentioned above. However, below $\delta\sim 0.57$, a distinct peak
emerges from a flat band of frequencies with decreasing $\delta$. At the
lowest $\delta$ values shown, these peaks indicate the shift of energy to
intra-well oscillations as temperature decreases for bias near critical. This
manifests as a “ringing” in the early time-dependent diffusion reminiscent of
the transient, decaying oscillations of the damped linear harmonic oscillator
subject to an impulsive load. Here, however, the frequency corresponds to
region III non-linear limit cycle oscillations thermally coupled to the sinks
and saddles of region I.

Figure 10: (a) Power spectra and (b) diffusion coefficient for several
temperatures near critical bias, i.e. the bifurcation line $\alpha=\beta=1$.
The top curves in the power spectra (for $\delta=1$) are for
$\beta=1,1.25,1.5$ (solid, dotted, dashed, respectively) and show little
variation from a $1/f$ “pink noise” spectrum.
Figure 11: Power spectra for $\delta=0.25$ ($\alpha=1$) case at several beta
showing state just above max ($\beta=1$) through $\beta=1.5$, where second
harmonic generation is evident

In Fig. 11 frequency peaks increase in height with increasing bias,
demonstrating a shift of energy to intra-well oscillations with bias exceeding
critical at low temperature. As Figure 12 demonstrates, corresponding peaks
for the zero-noise system can be expected to come from frequencies
approximately equal to the nondimensional bias value. However, it is notable
in Fig. 11 that peak frequencies at low noise seem to shift lower than their
corresponding zero-bias values with decreasing bias. At $\beta>\alpha$, the
peak frequencies correspond very well to zero-noise peaks, and the non-linear
second harmonic emerges from the flat band of frequencies above this. The
power spectra for $\beta=1.5$ is shown in Fig. 12 along with the time-
dependent diffusion coefficient for a range of $\delta$. The zero-noise power
spectrum is also shown, where the dynamics clearly appear to converge as
$\delta\rightarrow 0$. The limit cycle dynamics of the damped, driven pendulum
provides a useful, mechanistic understanding of the nature of the “running”
states. Importantly, the velocity has periodic fluctuations about a well-
defined average, with multiple time-scales associated with typical non-linear
oscillations. These can be seen in Figs. 11 and 12. However, the distinctions
between the meanings of the coordinates aside, the Brownian system cannot
follow true limit cycles, as the trajectories are nowhere differentiable, and
hence not unique. This allows for coupling bands of nearby trajectories via
thermally-induced transitions, hence broader peaks and shoulders in comparison
to the sharper spectrum of zero-noise dynamics (see Fig. 12.)
This suggests, rather than a phase space containing distinct multi-stable
“running” velocity states, the topology more closely resembles that of the
damped, driven pendulum (taking obvious liberty with the meaning of $y(t)$.)
The diffusion coefficient displays strong oscillatory ringing at the lowest
$\delta$, which is, in principle, experimentally observable in the mean-square
deviation at early times for a system with a very high barrier, $U_{o}\sim
32k_{B}T$ and larger for $\alpha=1$, driven above giant diffusion.

Figure 12: (a) Power spectra with zero-noise case shown (dashed) and (b) time-
dependent diffusion coefficient for $\alpha=1$, $\beta=1.5$, and several
$\delta$ showing onset of bifurcation/node coalescence and subsequent second
harmonic generation.

##  IV IV. Discussion and Conclusions

In this paper, we have revisited the classic problem of underdamped Brownian
motion in a tilted washboard potential via extensive numerical simulations,
complementing our discussion with a review of the corresponding zero-noise
dynamics and the related mulitstable Kramers problem. A picture emerges of
stochastic dynamics in a complex three-dimensional parameter space, where now,
in addition to the previously known friction vs. bias phase diagram [13], a
well-depth vs. temperature phase diagram is located. Anomalous diffusive
behavior at low temperature and low bias results in non-monotonic temperature
response and negative differential diffusivity; below we clarify how these
arise from the inertia of the Brownian particle. At very low bias, a curious
temperature-dependent maximum in the forward-biased jump rate is demonstrated,
which annihilates a temperature-dependent subcritical minimum in forward-
biased jump rate as the temperature increases across this transition. The low-
temperature dynamics of the system above the critical bias reveal ringing and
second-harmonic generation in the velocity power spectra, hallmarks of
stochastic resonant dynamics which demonstrate a convergence of the full
dynamics with the zero-temperature limit to complement Cheng and Yip’s
analytic work on the average velocity [25].
The phase diagram of Fig. 3 depicts low- and high-temperature– equivalently,
high- and low-barrier– driven jump-diffusion regimes. For the unbiased
periodic potential, this transition corresponds to the barrier heights below
which (or equivalently, temperatures $\delta$ above which) the frequency of
the linearized damped oscillations about the shifted well minimum approach the
magnitude of the typical inverse time-of-flight for a particle escaping a well
at thermal velocity. At lower barriers/higher temperatures, the time-of-flight
decreases, the jump rates increase (see Fig. 6), and the probability for
multiple jumps increases [28]. In our scaling, this is consistent with the
system approaching the regime of very high damping for all values of forcing
[28, 22]111This holds for temperature-independent friction and well-depth.,
i.e., above about $U_{o}\sim 6k_{B}T$ where, in terms of the Kramers problem,
inertial effects become negligible [28, 22]. Above this, the observed low-bias
negative differential mobility disappears; the mean forward-biased jumps at
zero bias exceed unity and a maximum in the uphill jump appears in the
vicinity of the giant diffusion maximum. Our analysis further allows the phase
diagram of Speciowicz and Luczka [13] to be identified with regimes of
persistent subdiffusion at low bias versus a relatively rapid approach to
normal diffusion above critical bias, corresponding to the well-known giant
diffusion effect [9, 10].
It is notable that the dynamics of the system are distinct in each region of
the $\alpha$-$\beta$ and $\alpha$-$\delta$ parameter spaces of Fig. 7 and Fig.
3. For sub-critical $\beta$, the system possesses long-lived “locked”
trajectories, wherein the particles with too low a kinetic energy to cross the
barrier remain trapped and oscillate in a single well for a very long time
between thermally-activated jumps. At the lowest values of $\beta$ well below
the thermal transition, the energy distribution is distorted Maxwellian. The
most energetic mode corresponds to the frequency of simple harmonic
oscillations about the stretched well minimum, and as seen in Fig. 8, the
frequency decreases concurrent with the increased energy distributed to these
low frequency modes. Thus the negative differential mobility and decrease to a
minimum in jump rate (see also Fig. 6) may be also attributed to inertia; the
additional energy input due to increasing bias is effectively stored in the
low-frequency simple harmonic oscillator modes corresponding to the
linearized, zero-temperature inertial dynamics.
Based on the previous discussion of the crossover of the system into the very
high damping regime, it appears that the inertia of the Brownian particle is
responsible for the persistent subdiffusion, and hence non-monotonic
temperature response for $\beta<\alpha$. This occurs without access to stable
“running” states until $\beta\rightarrow\alpha$, with rare well-to-well and
transitions increasing with $\delta$. On the other hand, the increase in
$\delta$ beyond about $\sqrt{\alpha/3}$ results an increase in the jump rate
(see Fig. 6a), causing more rapid decay of the diffusion coefficient after the
initial superdiffusive transient (see Fig. 5). So, while the inertial effects
themselves are the main cause of the persistent subdiffusion, the rapid decay
is facilitated by the decrease in barrier height relative to temperature
leading to increasing multiple jump probability and rate, which itself is a
rapidly changing function of friction in the underdamped regime [28]. It is
worth remarking that the friction-bias (or equivalently well depth-bias) and
well depth-temperature phase diagrams suggest at least some phenomenological
connection between the biased diffusion in a periodic potential and kinetic
phase transitions such as glassy behavior and gelation in colloids [32, 33].
Below a temperature threshold, there are infinite-waiting time “locked” states
that can only couple thermally to the “running” states around a threshold
bias; at higher temperatures, thermal coupling occurs for any non-zero bias.
By considering the velocity power spectra in these various regimes of driven
diffusion, the role of the stability of the underlying zero-noise fixed point
dynamics becomes clear. The giant diffusion effect corresponds to coupling of
trajectories across the I-III bifurcation via thermal fluctuations. Above the
bifurcation, characteristic frequencies emerge which show the convergence of
full, time-dependent dynamics towards those of the noise-free case. This
manifests as transient “ringing” in the time-dependent diffusivity and non-
linear harmonic generation with decreasing temperature evident in the power-
spectra for bias above the bifurcation. The evident ringing behavior is of
course also reminiscent of inertial systems; as the non-inertial, zero-
temperature dynamics do not possess these same oscillatory modes. Overall,
this response is consistent with a wide array of stochastic resonant phenomena
in periodically-driven systems [34, 35], and critically, it extends previous
analytic results concerning the average velocity [25].
Taken together, and contrasted against the overdamped case [14], our analysis
suggests inertial effects play a crucial role in the persistent subdiffusion,
negative differential mobility, low-bias non-monotonic temperature dependent
diffusivity, and even in the specifics of the stochastic resonant dynamics of
driven diffusive systems. It is usually taken for granted that inertial
effects for diffusion in a potential can be neglected at reasonably long
times, or equivalently, when the potential doesn’t vary appreciably over the
characteristic Kramers diffusion length [22]. However, in cases of periodic
forcing, not only is it feasible that potentials may indeed vary greatly over
such a typical length, but moreover any such steep variation occurs twice per
period, contributing non-negligibly to the long-time averages.
As stated previously the problem of a Brownian particle in a titled washboard
potential provides a dynamic model for many non-equilibrium statistical
systems. One particularly well-suited system for a driven Brownian particle is
that of an active colloid in a periodic optical array. In principle, with the
tunability of optical traps, it should be possible to achieve $U_{o}\geq
32k_{B}T$ to observe the transient ringing in the particle mean-square
displacement. While this may or may not be observable in a practical colloid-
in-liquid experiment, colloids in gasses or dusty colloidal plasmas [19, 17,
18] would be an ideal platform. Such an experiment could broadly test the
correspondence between the zero-noise dynamics and the high-barrier dynamics.
Beyond that, some of the phenomena elucidated herein may find applications in
fields such as colloidal assembly or in understanding collective effects in
active colloids, where inertial effects have begun to attract attention [19,
17, 18, 20]. From a theoretical perspective, while we have considered the
velocity dynamics near the I-III bifurcation, the nature of the velocity
dynamics near the I-II and II-III may also yield interesting behavior, and the
origin of the curious maximum in the forward-biased jump rate at low bias
warrants further study, possibly using analytic approaches similar to those of
Melnikov [29, 22].

###### Acknowledgements.

This material is based upon work supported by the National Science Foundation
under Grant Nos. 2001078 and 2040600 and the work utilized resources from the
University of Colorado Boulder Research Computing Group, which is supported by
the National Science Foundation (awards ACI-1532235 and ACI-1532236), the
University of Colorado Boulder, and Colorado State University.

## References

  * Evstigneev  _et al._ [2008] M. Evstigneev, O. Zvyagolskaya, S. Bleil, R. Eichhorn, C. Bechinger, and P. Reimann, Phys. Rev. E 77, 041107 (2008).
  * Volpe and G. [2013] G. Volpe and V. G., Am. J. Phys. 81, 224 (2013).
  * Volpe  _et al._ [2008] G. Volpe, S. Perrone, J. M. Rubi, and D. Petrov, Phys. Rev. E 77, 051107 (2008).
  * Fulde  _et al._ [1975] P. Fulde, L. Pietrono, W. R. Schneider, and S. Str$\ddot{a}$ssler, Phys. Rev. Lett. 35, 1776 (1975).
  * Peskin and Elston [2006] C. S. Peskin and T. C. Elston, SIAM J. Appl. Math 60, 842 (2006).
  * Ferrando  _et al._ [1992] R. Ferrando, R. Spadacini, and G. E. Tommei, Phys. Rev. B 45, 444 (1992).
  * Lindenberg  _et al._ [2005] K. Lindenberg, A. M. Lacasta, J. M. Sancho, and A. H. Romero, New J. Physics 7, 1 (2005).
  * Levi  _et al._ [1978] M. Levi, F. C. Hoppensteadt, and W. L. Miranker, Quart. Appl. Math. 36, 167 (1978).
  * Reimann  _et al._ [2001] P. Reimann, C. Van den Broeck, H. Linke, P. Hanggi, J. M. Rubi, and A. Perez-Madrid, Phys. Rev. Lett. 87, 010602 (2001).
  * Reimann  _et al._ [2002] P. Reimann, C. Van den Broeck, H. Linke, P. Hanggi, J. M. Rubi, and A. Perez-Madrid, Phys. Rev. E 65, 031104 (2002).
  * Lindner and Sokolov [2016] B. Lindner and I. M. Sokolov, Phys. Rev. E 93, 042106 (2016).
  * Marchenko  _et al._ [2018] I. G. Marchenko, I. I. Marchenko, and A. V. Zhiglo, Phys. Rev. E 97, 012121 (2018).
  * Spiechowicz and Luczka [2020] J. Spiechowicz and J. Luczka, Phys. Rev. E 101, 032123 (2020).
  * Berezhkovskii and Dagdug [2019] A. M. Berezhkovskii and L. Dagdug, J. Chem. Phys. 151, 131102 (2019).
  * Slapik  _et al._ [2019] A. Slapik, J. Luczka, P. Hänggi, and J. Spiechowicz, Phys. Rev. Lett 122, 070602 (2019).
  * Wang [2019] J. Wang, IOP Conf. Ser. Earth Environ Sci. 237, 052005 (2019).
  * Lowen [2020] H. Lowen, J. Chem. Phys. 152, 040901 (2020).
  * Lisin  _et al._ [2021] E. Lisin, O. Vaulina, I. Lisina, and O. Petrov, Phys. Chem. Chem. Phys 23, 16248 (2021).
  * Mandal  _et al._ [2019] S. Mandal, B. Liebchen, and H. Lowen, Phys. Rev. Lett 123, 228001 (2019).
  * De Karmakar and Ganesh [2020] S. De Karmakar and R. Ganesh, Phys. Rev. E 101, 032121 (2020).
  * Stratonovich [1957] R. L. Stratonovich, Soviet Physics Doklady Vol. 2, 416 (1957).
  * Coffey  _et al._ [1996] W. T. Coffey, Y. P. Kalmykov, and J. T. Waldron,  _“The Langevin Equation With Applications in Physics, Chemistry, and Electrical Engineering”_ (World Scientific, 1996).
  * Risken [1989] H. Risken,  _“The Fokker-Planck Equation Methods of Solution and Applications”_ (Springer-Verlag, 1989).
  * Jung and Risken [1984] P. Jung and H. Risken, Z. Phys. B Cond. Mat. 54, 357 (1984).
  * Cheng and Yip [2015] L. Cheng and N. K. Yip, Physica D 297, 1 (2015).
  * Bhattacharya [1999] K. Bhattacharya, Proc. Math. Phys. Eng. Sci. 455, 757 (1999).
  * Lifson and Jackson [1962] S. Lifson and J. L. Jackson, J. Chem. Phys. 36, 2410 (1962).
  * Ferrando  _et al._ [1993] R. Ferrando, R. Spadacini, and G. E. Tommei, Phys. Rev. E 48, 02437 (1993).
  * Melnikov [1991] V. I. Melnikov, Phys. Rep. 209, 1 (1991).
  * Strogatz [2015] S. H. Strogatz,  _“Nonlinear Dynamics and Chaos: wth Applications to Physics, Biology, Chemistry, and Engineering”_ (Westview Press, 2015).
  * Yan  _et al._ [2016] Y. Yan, Y. Wang, S. Senapati, J. Schiffbauer, G. Yossifon, and H.-C. Chang, Phys. Rev. E 94, 022613 (2016).
  * [32] P. Segre, V. Prasad, A. Schofield, and D. A. Weitz, Phys. Rev. Lett 86, 6042.
  * [33] M. Sztucki, T. Narayanan, G. Belina, A. Moussaid, F. Pignon, and H. Hoekstra, Phys. Rev. E 74, 051504.
  * Jung [1993] P. Jung, Phys. Rep. 234, 175 (1993).
  * Schmid and Hänngi [2005] G. Schmid and P. Hänngi, Physica A 351, 95 (2005).