Decay 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