Table Of ContentDecay to the nonequilibrium steady state of
6 the thermal diffusion in a tilted periodic
0
0
potential
2
n
a
J
0
2 T.Monnaia, A.Sugitab, J.Hirashimab, and K.Nakamurab
]
h
c aDepartment of Applied Physics ,Waseda University, 3-4-1 Okubo,
e
Shinjuku-ku,Tokyo 169-8555,Japan
m
- bDepartment of Applied Physics, Osaka City University, 3-3-138 Sugimoto,
t
a Sumiyoshi-ku, Osaka 558-8585, Japan
t
s
.
t
a
m
-
d Abstract
n
o
We investigate asymptotic decay phenomenon towards the nonequilibrium steady
c
[ state of the thermal diffusion in the presence of a tilted periodic potential. The pa-
2 rameter dependenceof thedecay rateis revealed byinvestigating theFokker-Planck
v (FP) equation in the low temperature case under the spatially periodic boundary
1
condition (PBC). We apply the WKB method to the associated Schro¨dinger equa-
3
4 tion. While eigenvalues of the non-Hermitian FP operator are complex in general,
1 in a small tilting case the imaginary parts of the eigenvalues are almost vanishing.
0
Then the Schro¨dinger equation is solved with PBC. The decay rate is analyzed in
6
0 the context of quantum tunneling through a triple-well effective periodic potential.
/ In a large tilting case, the imaginary parts of the eigenvalues of FP operator are
t
a
crucial. We apply the complex-valued WKB method to the Schro¨dinger equation
m
with the absorbing boundary condition, finding that the decay rate saturates and
-
d depends only on the temperature, the potential periodicity and the viscous con-
n
stant. The intermediate tilting case is also explored. The analytic results agree well
o
c with the numerical data for a wide range of tilting.
:
v
i
X
Key words: decay rate, thermal diffusion, tilted periodic potential, WKB analysis,
r Fokker-Planck equation
a
Email addresses: [email protected] (T.Monnai),
[email protected](A.Sugita),
[email protected](K.Nakamura).
Preprint submitted to Elsevier Science 6 February 2008
1 Introduction
Thermal diffusion process in a periodic potential in the presence of external
force appears in many situations. In solid-state physics, it appears in diffu-
sion of ions or molecules on crystal surfaces[1], Josephson junctions[2], motion
of fluxons in superconductors[3], rotation of molecules in solids[4], superionic
conductors[5], charge density waves[6], to mention a few. As a ratchet system,
it is also applied for particle selection[7], biophysical processes[8] and intra-
cellular transport[9][10]. Most of the researches in this context are concerned
with the steady state. For example, it is shown that a symmetric dichotomous
noise causes non-zero steady state current in such a system[11]. The existence
of an optimal diffusion coefficient vs. external force is also reported[12].
Such a diffusion process can be described by the over-damped Langevin equa-
tion with a tilted periodic potential,
∂ 2πW
ηx˙(t) = U(0)(x)+ x+ξ(t),
−∂x L
ξ(t) = 0, ξ(t)ξ(t′) = 2ηθδ(t t′). (1)
h i h i −
Here, η, θ = k T and W are the viscous coefficient, temperature and the
B
external force, respectively. ξ(t) is the Langevin force which satisfies the fluc-
tuation dissipation theorem. The periodic potential U(0)(x + L) = U(0)(x)
with period L is tilted by the external tilting force F = 2πW and we call the
L
potential U(x) U(0)(x) 2πWx as a tilted potential.
≡ − L
In order to investigate the probability density P(x,t) that the particle is found
in a position x at time t, it is useful to transform the Langevin equation into
the corresponding Fokker-Planck(FP) equation,
∂ 1 ∂ ∂U(x) ∂
P(x,t) = +θ P(x,t). (2)
∂t η∂x ∂x ∂x!
Hereafter we solve FP equation (2) under the (spatially) periodic boundary
condition(PBC),P(x+L,t) = P(x,t),whichissatisfiedbythesteady stateat
t . The analytic expression for the steady state Pst(x) lim P(x,t)
t→∞
→ ∞ ≡
under PBC is well-known[9,13,14]:
x+L
η
Pst(x) = N e−U(θx) dyeU(θy),
θ
Zx
2
L x+L −1
θ U(y)−U(x)
N = dx dye θ , (3)
η
Z0 Zx
where N is the normalization constant chosen so that LdxPst(x) = 1. The
0
steady state has a non-zero probability current Jst(x) given by
R
1 ∂
Jst U(x) θ Pst(x) = N 1 e−2πθW . (4)
≡ η − ∂x! −
(cid:16) (cid:17)
The above steady state is a nonequilibrium state. In this paper, we study the
low-temperature relaxation dynamics towards the steady state in the presence
of a tilted periodic potential. In particular, asymptotic behaviors of the decay
phenomenon is investigated by theeigenfunction expansion method.Our main
interest lies in parameter dependence of the decay rate.
This paper is organized as follows. In section 2 we present the framework
of our analysis, relevant numerical data for eigenvalues of the FP operator,
and the main result. Section 3 treats the important theme on the boundary
condition for eigenstates of the associated Schr¨odinger equation. In subsection
3.1 the small tilting case is investigated in the context of a tunneling problem.
In subsection 3.2 the large tilting case is treated, where we shall apply the
complex-valued WKB method to the Schr¨odinger equation with the absorbing
boundary condition deduced from PBC for the original FP equation. Section
4 is devoted to a summary. Appendix is concerned with detailed calculations
in section 3.2.
2 Eigenvalues of the Fokker-Planck operator
The WKB treatment of FP equation has a long history since the monumental
work by van Kampen[15]. Most of the studies, however, were limited to the
system with a well-defined thermal equilibrium state[16]. This fact guaran-
teed a real-valued nature of eigenvalues for the non-Hermitian FP operator,
and one can reduce FP equation to the eigenvalue problem of the associated
Schr¨odinger equation with the same natural boundary condition as the FP
equation has. On the other hand, relaxation dynamics towards the nonequi-
librium steady state with the non-vanishing constant current demands a qual-
itatively different WKB approach. In fact FP operator in such a system has
complex eigenvalues in general.
Thissituationisrecognizedinoursystemwithatiltedsinusoidal potential[17],
2πx 2πW
U(x) = U cos x, (5)
0
L − L
(cid:18) (cid:19)
3
where U is the amplitude of the periodic part of the potential. Putting eigen-
0
values of the non-Hermitian FP operator as E and substituting P(x,t) =
−
e−EtP(x) into Eq.(2), we have the eigenvalue problem,
1 ∂ ∂U(x) ∂
EP(x) = +θ P(x). (6)
− η∂x ∂x ∂x!
While the steady state in Eq.(3) is periodic, Pst(x+L) = Pst(x), and is char-
acterized by the zero Bloch wavenumber (k = 0), it is difficult to analyze the
relaxation dynamics in general. Any initial distribution P(x,0) which include
k = 0 Fourier components can relax towards Pst(x): the k = 0 components
6 6
will decay out on the way of time evolution. In the present work, however, we
confine ourselves to the manifold of k = 0, which makes the analysis of the
decay phenomenon accessible. This implies that the time-dependent solution
of FP equation in Eq.(2) relaxing towards the steady state should satisfy the
spatial PBC, P(x+L,t) = P(x,t) and be constructed from solutions with a
natural boundary condition Q(x,t) (lim Q(x,t) = 0) as
x→±∞
∞
P(x,t) = Q(x+nL,t). (7)
n=−∞
X
Since the FP operator in Eq.(2) is invariant against space translation by L, it
hasBloch-typeeigenstates P (x) and(complex)eigenvalues E .Noting
n,k n,k
{ } { }
the normalization matrix 0LdxPk∗,n(x)Pk′,n′(x) = δk,k′Nn(k,n)′ due to the non-
Hermitian nature of FP operator, the solution P(x,t) in Eq.(7) is now written
R
explicitly as
P(x,t) = C P (x)e−En,k=0t, (8)
n,k=0 n,k=0
n
X
L
C = (N(0))−1 dxP∗ (x)P(x,0), C = 0, (9)
n,k=0 n,n′ n′,k=0 n,k6=0
n′ Z
X 0
which starts from an arbitrary k = 0 distribution P(x,0) at t = 0 and relaxes
towards Pst(x).
In the manifold with k = 0, the probability distribution can be expanded as
P(x) = ncne2πLixn. Then the eigenvalue problem (6) is reduced to
P
2π 2
A c = ηEc
n,m m n
L
(cid:18) (cid:19) m
X
nU nU
A = inW +θn2 δ 0δ + 0δ , (10)
n,m n,m n−1,m n+1,m
− − 2 2
(cid:16) (cid:17)
4
400 Im
200
Re
200 400 600 800
-200
-400
Fig. 1. Distribution of the complex eigenvalue with the second smallest real part
for full parameter range 3 U 20,0 W 20 (dots), and for the case of fixed
0
≤ ≤ ≤ ≤
potential barrier U = 10,0 W 20 (solid line). A pair of solid lines rapidly
0
≤ ≤
merge to the real axis as W decreases below U .
0
which can be solved numerically by truncating the infinite dimensional matrix
A to a large but finite dimendional matrix.
The zero-eigenvalue (E = 0) corresponds to the steady state. Under PBC, as
is obvious from Eq. (8), the probability distribution approaches exponentially
to the unique steady state, and the decay rate λ is obtained as the real part of
the second smallest eigenvalue. Among many complex eigenvalues of Eq.(10),
therefore, we have focused on the one with the second-smallest real part and
its variation against various parameter values U ,W (see Fig.1). From Fig.1
0
one can see that the eigenvalues with large positive real parts have nearly zero
imaginarypart.(SeetailsofsolidlinesinFig.1.)These eigenvalues correspond
to the small tilting case. In fact, a pair of solid lines rapidly merge to the real
axis as W decreases below U , ensuring the real-valued nature of eigenvalues
0
in the case of small tilting.
U and W dependence of the decay rate λ is shown in Fig. 2, which we have
0
confirmed bynumerically solving theFPequationinEq.(2). Thereisa smooth
but steep decrease of the decay rate around the crossover line of U W,
0
≃
where the original potential minima disappear.
In the next section, the Schr¨odinger equation associated with FP equation
will be solved with use of the WKB method. In the small tilting case in
section 3.1, we impose PBC for the wavefunction of the Schr¨odinger equa-
tion and apply the Bloch theorem. The decay rate for W < U is shown to
0
2 2
be 2π 1 W U . In the large tilting case in subsection 3.2, we have
L r − U0 0
reco(cid:16)urs(cid:17)e to the(cid:16)abs(cid:17)orbing boundary condition and apply the WKB method
which admits complex energies. The decay rate for W U is found to be
0
≫
(2π/L)2θ (with no W and U dependence). The crossover region where W is
0
comparable to U is numerically analyzed based on this approach.
0
5
800
600
decayrate
400
15
200
0
10
Uo
55
1100 5
ww 1155
2200
Fig. 2. W and U dependence of the decay rate λ.
0
120
100
80
60
40
W<U0
20 W>U0 W>>U0
W
5 10 15 20
Fig. 3. W dependence of the decay rate λ for U = 3,θ = 0.1. Dots show the decay
0
rates obtained from numerical diagonalization of FP operator. Three solid lines are
theoretical results (For small tilting case W < U the decay rate λ is obtained as
0
the first excited level ~ω of the associated Schro¨dinger equation (left part). For
a
crossover regime W U , λ is given as the numerical solution of (37)(middle part).
0
∼
For large tilting case W U , theperturbativesolution of (38) gives λ (right part).
0
≫
3 WKB analysis of probability distribution
The FP equation is transformed into the Schr¨odinger equation by the separa-
tion ansatz[15]
P(x,t) = e−U2(θx)φ(x)e−Et. (11)
Then one has
d2
θ φ(x)+(E V(x))φ(x) = 0, (12)
dx2 −
6
where we redefined eigenvalue as Eη E, and defined the effective potential
→
′ 2 ′′
U (x) U (x)
V(x) θ . (13)
≡ 2θ ! − 2θ
The Schr¨odinger operator (SO) in Eq.(12) looks Hermitian. Noting in Eq.(11)
U(x+L) = U(x) 2πW andP (x+L) = eikLP (x),however, thewave func-
n,k n,k
−
tion φ(x) of the Schr¨odinger equation should satisfy the absorbing boundary
condition,
φ(x+L) = eikL−2π2θWφ(x), (14)
which eventually breaks the Hermitian nature of SO.
In the following we shall analyze the Schr¨odinger equation in Eq.(12) first in
the small tilting case and then in the strong tilting case. In the small tilting
case, the numerical evidence in the previous section has shown the real-valued
nature of the second-lowest eigenvalue, indicating the recovery of Hermitian
nature of SO. Therefore we reduce Eq. (14) to
φ(x+L) = eikLφ(x), (15)
by assuming the absorbing term 2πW being vanishing. The validity of this
2θ
assumption will be verified by comparing between the analytic and numerical
decay rates. Below we shall confine to the k = 0 manifold.
3.1 Small tilting case
Substituting the tilted sinusoidal potential (5) into (13), one has the effective
potential,
π2 2πx 2 1 2π 2 2π
V(x) = U sin +W + U cos x. (16)
L2θ 0 L 2 L 0 L
(cid:18) (cid:19) (cid:18) (cid:19)
Since the first term is much larger than the second one in the low temperature
regime θ U , the second one is omitted hereafter.
0
≪
InsolvingEq.(12)thelowtemperatureconditionallowsustousetheWKBap-
proximation, since θ corresponds to ~2 in the standard Schr¨odinger equation.
2m
In the case that the external force is weak (W < U ), the effective potential
0
7
1500
1250
1000
V 750
500 X0X1 X2X3 X4X5
250 b d
a c e
0
-0.5 0x 0.5 1
-250
Fig. 4. Effective potential V(x) and the real part of the eigenvalue of FP operator
with second smallest real part E with the corresponding classical turning points
x ,....,x for W = 1,U = 3.
0 5 0
U
x
0.5 1 1.5 2 2.5 3
-20
a
-40
b
-60
c
-80
Fig. 5. Original tilted potential U(x) in real space. Small tilting case (a), crossover
region (b), and large tilting case (c).
V(x) has three wells. We consider the energy levels for the eigenstates which
experience all the three wells[18] (see Fig. 4).
In terms of the tilted potential U(x), the weak external force guarantees the
existence of potential minima in Fig. 5. Let us denote the wave functions in
each region a, ,e bordered by classical turning points x , ,x as φ , ,φ
0 5 a e
··· ··· ···
and introduce the action integrals
1 x1 1 x2
S = pdx, M = pdx
a ~ b ~
xZ0 xZ1
1 x3 1 x4
S = pdx, M = pdx, (17)
c ~ d ~
xZ2 xZ3
where ~2 = 2mθ and p = 2m E V(x) (p stands for the momentum in a
| − |
potential-well region). S anqd M are defined in each well and barrier, respec-
tively.Eachneighboringwavefunctionsaremutuallyrelatedbytheconnection
8
formula
1 1 i
e±i(S+π4) (E > V(x)) eS e−S (E < V(x)), (18)
√p ←→ √p ± 2
(cid:18) (cid:19)
where
y0
1
S = pdx , (19)
~
| |
Zx
with y being the classical turning point.
0
For example, φ and φ are related as
a b
1 −i x pdx i x pdx
φa(x)= c1e ~ x0 +c2e~ x0
√p
(cid:18) R R (cid:19)
φb(x)= 1 c1e−iSa e~1 xx1pdx + ie−~1 xx1pdx (20)
p (cid:26) (cid:18) R 2 R (cid:19)
| |
+qc2eiSa e~1 xx1pdx ie−~1 xx1pdx (21)
− 2
(cid:18) R R (cid:19)(cid:27)
In the same way, one obtains the expression for φ (x) as
e
1 −i x pdx i x pdx
φe(x) = (Ac1 +Bc2)e ~ x4 +(Cc1 +Dc2)e~ x4 , (22)
√p
(cid:26) R R (cid:27)
where A,B,C and D are given as
1 1 i
A=eiSa ieMb sinS + e−MbcosS e−Md +2 eMbcosS + e−MbsinS eMd ,
c c c c
2 4 4
(cid:26) (cid:18) (cid:19) (cid:18) (cid:19) (cid:27)
1 i 1
B=e−iSa eMbsinS + e−MbcosS e−Md 2 ieMb cosS + e−MbsinS eMd ,
c c c c
2 4 − 4
(cid:26) (cid:18) (cid:19) (cid:18) (cid:19) (cid:27)
1 i 1
C=eiSa eMbsinS e−MbcosS e−Md 2 ieMb cosS + e−MbsinS eMd
c c c c
2 − 4 − − 4
(cid:26) (cid:18) (cid:19) (cid:18) (cid:19) (cid:27)
1 1 i
D=e−iSa ieMb sinS + e−MbcosS e−Md +2 eMbcosS e−MbsinS eMd
c c c c
2 − 4 − 4
(cid:26) (cid:18) (cid:19) (cid:18) (cid:19) (cid:27)
(23)
On the other hand, with use of the Bloch theorem in Eq.(15) based on the
periodicity of the effective potential V, φ and φ should be related as
e a
φ = eikLφ , (24)
e a
9
where k is the wave number. This relation gives a constraint to guarantee the
nontrivial coefficients (c ,c ) of φ and φ :
1 2 a e
A eikL B
− = 0, (25)
(cid:12) (cid:12)
(cid:12) C D eikL(cid:12)
(cid:12) (cid:12)
(cid:12) − (cid:12)
(cid:12) (cid:12)
(cid:12) (cid:12)
(cid:12) (cid:12)
which can be rewritten as
1 1
eikL = cos(kL)+isin(kL) = (A+D) i 1 (A+D)2. (26)
2 ± s − 4
Here we have used the fact AD BC = 1. At low temperatures the tunneling
−
integrals are sufficiently large, M 1 and M 1. Then we can make the
b d
≫ ≫
approximation:
1
(A+D)2 = 4e2Mbe2Md cos2(S ) cos2(S ), (27)
∼ a c
4
Since (27) should satisfy (26), cosS and cosS must be sufficiently small
a c
| | | |
and be approximated as
1
cos(S ) ( 1)na,c S (n + ) π
a,c a,c a,c
≃ − − 2
(cid:26) (cid:27)
sin(S ) ( 1)na,c. (28)
a,c
≃ −
Furthermore, using the harmonic approximation in the neighborhood of the
a,c
potential minima x , the actions can be evaluated as
min
π
S = (E V ), (29)
a,c ~ω − 0a,c
a,c
where ω = 2θV′′(xa,c ) is the curvature of the potential minima, and V
a,c ~2 min 0a,c
stands for thqe local minima of V(x). Then Eq.(26) gives the quantization
condition for eigenvalues E
1
cos(kL) = (A+D)
2
1 1
= 4eMb+Md + e−Mb−Md cos(S )cos(S ) e−Mb+Md +eMb−Md sin(S )sin(S )
a c a c
2 4 −
(cid:26)(cid:18) (cid:19) (cid:16) (cid:17) (cid:27)
1 1 1 π 1 π
4eMb+Md + e−Mb−Md (n + )π (E V ) (n + )π (E V )
≃ 2 4 a 2 − ~ω − 0a c 2 − ~ω − 0c
(cid:26)(cid:18) (cid:19)(cid:18) a (cid:19)(cid:18) c (cid:19)
eMb−Md +e−Mb+Md ( 1)na+nc. (30)
− −
(cid:16) (cid:17)o
10
