Diffusion of a Brownian Particle in a Periodic Potential with Memory Friction

Meng Hu, Zhanwu Bai^*

Department of Mathematics and Physics, North China Electric Power University, Baoding, China.

DOI: 10.4236/jmp.2022.138071                

Abstract

This study investigates the diffusive motion of a Brownian particle in a 1D periodic potential. The reactive flux theory for finite barriers and memory friction is developed to calculate the escape rate in the spatial diffusion regime. The diffusion coefficient is obtained in terms of the jump-model. The theoretical results agree well with the Langevin simulation results. The method can be generalized to other colored noises with Gaussian distribution.

Keywords

Escape Rate, Finite Barrier, Memory Friction

1. Introduction

The Kramers rate [1] theory and its successive development for all damping, the original Melnikov and Meshkov (MM) theory for white noise [2] , and Pollak, Grabert and Hanggi (PGH) theory for general noise [3] , are correct for high potential barriers. Nevertheless, the finite or the lower barrier is universal in practical stochastic processes. Two main approaches are employed for finite barrier correction. One is based on the generalized Langevin equation equivalent to Hamiltonian formulism, and the other is based on the Fokker-Planck equation. The perspective of the former is the finite barrier correction in the framework of PGH theory and improved PGH theory [4] . The improved PGH theory with finite barrier correction for Ohmic friction is quite accurate over the whole friction range up to the reduced barrier height 4. However, high accuracy agreement between theoretical and numerical exact results in the spatial diffusion regime cannot be expected, especially when the reduced barriers are low, because the perturbation parameter in the theory is no longer small. Within the framework of the Fokker-Planck processes in the spatial diffusion regime, a perturbation theory was developed to calculate the finite barrier escape rate by means of the flux over population expression, the Rayleigh quotient, and the mean first passage time to the stochastic separatrix [5] . Although a simple analytical perturbation solution can be obtained for a polynomial potential, obtaining a simple analytical solution for a nonpolynomial potential is difficult. The two approaches mentioned above have some limitations. The equivalent Hamiltonian approach is a perturbation theory for weak damping and extrapolation to spatial diffusion regime is questionable. In addition, obtaining a simple analytical expression from the finite barrier correction to the PGH theory is difficult. The Fokker-Planck equation approach relies on the existence of a Fokker-Planck equation. However, this situation is not always the case for a generalized Langevin equation with general noise and an arbitrary potential. This approach is also based on perturbation expansion, which is not suitable for low barriers.

To develop a simple, accurate and widely applicable approach to incorporate the finite barrier correction in analytical calculation of the escape rate, the reactive flux theory for finite barriers is proposed in [6] . The theoretical results match well with the simulation results until lower barriers for a Brownian particle moving in a cubic metastable potential and subjected to a Gaussian white noise.

The Brownian motion in a periodic potential can model many physical and chemical situations, such as the motion of ions in superionic conductors [7] , the Josephson supercurrent in tunneling junctions [8] [9] , the mass transport in solids [10] [11] , and the diffusion on surfaces [12] . Several analytical results with finite barrier correction were given in the spatial diffusion regime for the diffusion of a Brownian particle in a periodic potential. The escape rate for an internal Ornstein-Uhlenbeck noise has been investigated within the frameworks of Grote-Hynes [13] theory with high potential barriers and PGH theory [14] [15] for not too large damping and correlation time. In the present work, we extend our finite barrier correction scheme [6] to an internal Ornstein-Uhlenbeck noise to calculate analytically the diffusion coefficient in the spatial diffusion regime with finite barrier correction.

2. Reactive Flux Theory for Finite Barriers and Memory Friction

We consider a Brownian particle with unit mass diffusing in a 1D periodic potential under the influence of an internal Ornstein-Uhlenbeck noise. The particle is in contact with a heat bath at temperature T, which provides fluctuation and dissipation. The dynamics of the process is governed by the following generalized Langevin equation:

$\ddot{x}+{\displaystyle {\int }_0^t}\text{d}{t}^{\prime }\gamma \left(t-t^{\prime }\right)\stackrel{\dot{}}{x}\left(t^{\prime }\right)=-V^{\prime }\left(x\right)+\epsilon \left(t\right).$ (1)

where $\epsilon \left(t\right)$ is the Ornstein-Uhlenbeck noise, which is associated with the memory kernel function $\gamma \left(t\right)$ by the fluctuation-dissipation theorem:

$\langle \epsilon \left(t\right)\epsilon \left(t^{\prime }\right)\rangle =k_{B}T\gamma \left(t-t^{\prime }\right),$ (2)

where $k_B$ is the Boltzmann constant, $\gamma \left(t\right)$ is the memory kernel function, given by

$\gamma \left(t\right)=\frac{\gamma }{\tau }\exp \left(-t/\tau \right)\mathrm{,}$ (3)

where $\tau$ is the correlation time of the noise, and $\gamma$ is regarded as the effective friction coefficient due to $\gamma ={\displaystyle {\int }_0^{\infty }}\text{d}t\gamma \left(t\right)$ . In Equation (1), $V\left(x\right)$ is the 1D periodic potential, given by

$V\left(x\right)=V_0\cos x.$ (4)

What we want to investigate is the diffusion coefficient in the spatial diffusion regime, which can be attributed to the calculation of escape rate, because single jump is dominate in the spatial diffusion regime.

In the traditional reactive flux formulation, an initial equilibrium distribution is assumed for the trajectories of particles starting at the top of the barrier. The equilibrium state assumption at the top of the barrier is reasonable for higher potential barriers because the current is small. It is no longer a good approximation for lower potential barriers where the current is not small enough. For finite barriers, we remove the starting point of the trajectories to somewhere $x=x_0$ in the potential well where the probability distribution can be regarded as an equilibrium one, and replace the potential barrier from $x_0$ to $x_b$ (top of the barrier) with an equivalent parabolic barrier in the spatial diffusion regime. $x_0$ is given by $V\left(x_0\right)=V_b-k_{B}T$ , where $V_b=2V_0$ is the potential barrier height, and $k_{B}T$ is the average energy fluctuation of a quasi-equilibrium distribution. The equivalent potential $V_e\left(x\right)$ is given by

$V_e\left(x\right)=V_b-\frac{1}{2}m{\omega }_b^2{x}^2.$ (5)

The equivalent potential barrier frequency ${\omega }_b$ can be determined by the potential approach scheme in the barrier region, that is by minimization of the following average

$I={\displaystyle {\int }_{x_0}^{x_1}}\text{d}x{\left[V\left(x\right)-V_e\left(x\right)\right]}^2{P}_e\left(x\right),$ (6)

where $P_e\left(x\right)$ is the Boltzmann distribution normalized in barrier region, in which the potential $V_e\left(x\right)$ is used. $x_0\mathrm{,}{x}_1$ are two intersection points of the straight line $V=V_b-k_{B}T$ with the original potential. Similar to the derivation of the escape rate in original reactive flux theory, the expression of the escape rate for Gaussian white noise can be worked out [6] :

$k=\frac{k_{B}T}{nQ{\omega }_b}\left[\sqrt{{\omega }_b^2+{\gamma }^2/4}-\gamma /2\right]\exp \left(-\frac{V_b}{k_{B}T}\right).$ (7)

The population n in the expression can be approximated by

$n=\frac{1}{Q}{\displaystyle {\int }_{-\infty }^{x_0}}\text{d}x{\displaystyle {\int }_{-\infty }^{\infty }}\text{d}{v}_0\exp \left(-\frac{V\left(x\right)}{k_{B}T}\right)\exp \left(-\frac{v_0^2}{2m{k}_{B}T}\right),$ (8)

where Q is the partition function for the particles in the potential well.

When the theory is extended to an internal Ornstein-Uhlenbeck noise, a difficulty is encountered: the second order moments appearing in the probability density in barrier region is a oscillation function and the limits as time t tend to infinite do not exist, that is, the steady-state probability density does not exist. The steady-state probability density only exists for a small parameter region below the dashed line in Figure 1.

By using the method of characteristic function, the Fokker-Planck equation for transition probability density $P\left(x\mathrm{,}v\mathrm{,}t\right)$ in the barrier region is given as follows [16] :

$\frac{\partial P}{\partial t}=-v\frac{\partial P}{\partial x}-{\tilde{\omega }}_b^2\left(t\right)x\frac{\partial P}{\partial v}+{\tilde{\gamma }}_b\left(t\right)\frac{\partial (vP)}{\partial v}+{\psi }_b\left(t\right)\frac{{\partial }^{2}P}{\partial x\partial v}+{\varphi }_b\left(t\right)\frac{{\partial }^{2}P}{{\partial }^{2}v},$ (9)

where the subscript b signifies the dynamical quantities defined at the barrier top. The expressions of these quantities can be found in [16] . ${\tilde{\omega }}_b\left(t\right)$ and ${\tilde{\gamma }}_b\left(t\right)$ in the above equation can be called the renormalized (by the noise) potential barrier frequency and friction coefficient, and the original potential barrier frequency can be called bare frequency. For an internal Ornstein-Uhlenbeck noise, the steady-state transition probability density only exists in the region below the dashed curve in Figure 1. We use $s_i\left(i=1,2,3\right)$ to denote the three poles of the Laplace transformation of the coordinate Green function, i.e., the roots of the following characteristic equation:

$s^3+a{s}^2+bs+c=0,$ (10)

with $a=1/\tau ,b=\gamma /\tau -{\omega }_b^2,c=-{\omega }_b^2/\tau$ . For large time, i.e., $-\left(s_1+Re{s}_2\right)t\gg 1$ , ${\psi }_b\left(t\right)\mathrm{,}{\varphi }_b\left(t\right)$ in Equation (9) can be approximated as

${\psi }_b\left(t\right)=k_{B}T\left({\tilde{\omega }}_b^2\left(t\right)/{\omega }_b^2-1\right)\mathrm{,}{\varphi }_b\left(t\right)=k_{B}T{\tilde{\gamma }}_b\left(t\right)\mathrm{,}$ (11)

and ${\tilde{\omega }}_b^2\left(t\right)\mathrm{,}{\tilde{\gamma }}_b\left(t\right)$ are given by

Figure 1. The applicable parameter region of the theory is below the solid line. (a) for $T=0.4$ , (b) for $T=0.6$ .

${\tilde{\omega }}_b^2\left(t\right)=\frac{{\omega }_b^{2}Re{s}_2}{\tau {\left|s_2\right|}^2}+\frac{{\omega }_b^{2}Im{s}_2}{\tau {\left|s_2\right|}^2}tan\left(Im{s}_{2}t\right)\mathrm{,}$ (12)

${\tilde{\gamma }}_b\left(t\right)=-\left(s_1+Re{s}_2\right)+Im{s}_2\tan \left(Im{s}_{2}t\right)\mathrm{,}$

where $s_1$ and $s_2$ denote the roots of Equation (10) with the maximal and the secondary maximal real parts, respectively, and “Re” and “Im” denote the real and imaginary parts of the complex, respectively. If the root $s_2$ satisfies $\left|Im{s}_2\right|t\gg 1$ , we integrate two sides of the Fokker-Planck equation for large number periods, which results in an effective steady-state Fokker-Planck equation satisfied by an effective probability density, which is the long term average of the usual probability density and is just the quantity we need for the calculation of the mean escape rate for the time. If $\left|Im{s}_2\right|t\gg 1$ is not satisfied, the imaginary part of ${\tilde{\omega }}_b^2\left(t\right)$ , ${\tilde{\gamma }}_b\left(t\right)$ can be neglected because of the large barrier passing time t and the effective steady-state Fokker-Planck equation remains in the same form. Some details can be found in our previous work [17] . The barrier frequency ${\omega }_b$ in the above expressions is the renormalized one, which incorporates the finite barrier correction. The expression of the escape rate has the same form as (7), but ${\omega }_b^2\mathrm{,}\gamma$ are replaced with the renormalized long term average frequency and damping, given by

$\langle {\tilde{\omega }}_b^2\left(t\right)\rangle =-\frac{{\omega }_b^{2}Re{s}_2}{\tau {\left|s_2\right|}^2}\mathrm{,}$ (13)

$\langle {\tilde{\gamma }}_b\left(t\right)\rangle =-\left(s_1+Re{s}_2\right)\mathrm{,}$

and the expression of the escape rate is

$k=\frac{k_{B}T}{nQ\sqrt{\langle {\tilde{\omega }}_b^2\left(t\right)\rangle }}\left[\sqrt{\langle {\tilde{\omega }}_b^2\left(t\right)\rangle +{\langle {\tilde{\gamma }}_b\left(t\right)\rangle }^2/4}-\langle {\tilde{\gamma }}_b\left(t\right)\rangle /2\right]\exp \left(-\frac{V_b}{k_{B}T}\right)\mathrm{,}$ (14)

n is given by Equation (8). Because single-jump is dominant in the spatial diffusion regime, the diffusion coefficient is given by

$D=\frac{1}{2}\Gamma d^2=k{d}^2,$ (15)

where $\Gamma =2k$ is the total jump rate, and d is the spatial period of the potential, here $d=2\pi$ . For some parameters, $\langle {\tilde{\gamma }}_b\left(t\right)\rangle \le 0$ , which is physically unreasonable. The theory is applicable below the solid line (plotted with $\langle {\tilde{\gamma }}_b\left(t\right)\rangle =0$ ) in Figure 1. In contrast to the PGH theory [14] , the proposed method can also be applied to large damping and correlation time case provided that the parameter is below the solid line in Figure 1.

3. Diffusion Coefficient: Theory versus Simulation

To avoid a direct simulation of the generalized Langevin Equation (1), we introduce an auxiliary variable

$z=-{\displaystyle {\int }_0^t}\text{d}{t}^{\prime }\gamma \left(t-t^{\prime }\right)\stackrel{\dot{}}{x}\left(t^{\prime }\right)+\epsilon \left(t\right),$ (16)

which satisfies the following differential equation

$\stackrel{\dot{}}{z}+\frac{\gamma }{\tau }v+\frac{1}{\tau }z=0,$ (17)

in which the fluctuation-dissipation theorem (2) has been used, and v is the velocity. Hence, the generalized Langvin Equation (1) is equivalent to the following Markovian-type Langvin equations

$\stackrel{\dot{}}{x}=v$ (18)

$\stackrel{\dot{}}{v}=-V^{\prime }\left(x\right)+z$

$\stackrel{\dot{}}{z}=-\frac{\gamma }{\tau }v-\frac{1}{\tau }z$

We simulate the Langevin Equations (18) by the second-order Runge-Kutta algorithm. In the calculation, the natural unit ( $m=1$ , $k_B=1$ ), the dimensionless parameter $V_0=1$ , and the time steps are $\Delta t={10}^{-4}$ for $\tau =0.1$ and $\Delta t=5\times {10}^{-4}$ for $\tau =0.4$ . The test particles start from the a potential well and have zero velocity. The number of test particles $N=5\times {10}^5$ is used to describe the diffusion motion of a Brownian particle. The mean square displacement $\langle {\left(\Delta x\left(t\right)\right)}^2\rangle$ reveals a good linear relation at long times. The diffusion coefficient D is obtained by the long time behavior of the mean square displacement

$\langle {\left(\Delta x\left(t\right)\right)}^2\rangle =2Dt.$ (19)

Figure 2 shows that the theoretical results of the diffusion coefficient (Equation (14) and Equation (15)) match well with the Langevin simulation results. In

Figure 2. Diffusion coefficient: theory versus simulation. The solid lines are theoretical results, and the dashed lines are simulation results. Where $V_b=2$ , other parameters are: (a) $T=0.4,\tau =0.1$ ; (b) $T=0.4,\tau =0.4$ ; (c) $T=0.6,\tau =0.1$ ; (d) $T=0.6,\tau =0.4$ .

the spatial diffusion regime, specified by $\gamma \ge 2{\omega }_b$ , the maximal errors are about 4% for $T=0.4,\tau =0.4$ and 2% for other cases.

For higher temperatures or lower reduced potential barrier heights, the error increases, which can be interpreted as the jump-model being no longer a good approximation.

4. Summary

In the spatial diffusion regime, the calculation of diffusion coefficient is attributed to the calculation of the escape rate. The reactive flux theory for finite barriers is developed to incorporate finite barrier effect. The starting point of the Brownian particle is removed into the potential well where the probability density can be viewed as an equilibrium one, and the potential barrier is equivalent to a parabolic one. An equivalent steady-state Fokker-Planck equation is established to overcome the difficulty of the absence of the steady-state probability density. The theoretical results for diffusion coefficient indicate a good agreement with the Langevin simulation results in a certain range of parameters.

Fund

This work was supported by the National Natural Science Foundation of China under Grant No. 11905062 and the Fundamental Research Funds for the Central Universities under Grant No. 2019MS115.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Kramers, H. (1940) Physica, 7, 284-304. https://doi.org/10.1016/S0031-8914(40)90098-2
[2]	Melnikov, V.I. and Meshkov, S.V. (1986) The Journal of Chemical Physics, 85, 1018. https://doi.org/10.1063/1.451844
[3]	Pollak, E., Grabert, H. and Hanggi, P. (1989) The Journal of Chemical Physics, 91, 4073. https://doi.org/10.1063/1.456837
[4]	Pollak, E. and Ianconescu, R. (2016) The Journal of Physical Chemistry A, 120, 3155-3164. https://doi.org/10.1021/acs.jpca.5b11502
[5]	Talkner, P. (1994) Chemical Physics, 180, 199-213. https://doi.org/10.1016/0301-0104(93)E0426-V
[6]	Bai, Z.W. (2021) Physical Review E, 103, Article ID: 042144.
[7]	Dieterich, W., Fulde, P. and Peschel, I. (1980) Advances in Physics, 29, 527-605. https://doi.org/10.1080/00018738000101406
[8]	Bttiker, M., Harris, E.P. and Landauer, R. (1983) Physical Review B, 28, 1268. https://doi.org/10.1103/PhysRevB.28.1268
[9]	Ben-Jacob, E., Bergman, D.J., Matkowsky, B.J. and Schuss, Z. (1982) Physical Review A, 26, 2805. https://doi.org/10.1103/PhysRevA.26.2805
[10]	Munakata, T. and Kawakatsu, T. (1985) Progress of Theoretical Physics, 74, 262-271. https://doi.org/10.1143/PTP.74.262
[11]	Kawakatsu, T. (1986) Progress of Theoretical Physics, 75, 747. https://doi.org/10.1143/PTP.75.747
[12]	Agassi, D. and Eberly, J.H. (1985) Physical Review Letters, 54, 34. https://doi.org/10.1103/PhysRevLett.54.34
[13]	Grote, R.F. and Hynes, J.T. (1980) The Journal of Chemical Physics, 73, 2715. https://doi.org/10.1063/1.440485
[14]	Ianconescu, R. and Pollak, E. (2014) The Journal of Chemical Physics, 143, Article ID: 104104. https://doi.org/10.1063/1.4929709
[15]	Ianconescu, R. and Pollak, E. (2016) Faraday Discussions, 195, 111-138. https://doi.org/10.1039/C6FD00105J
[16]	RayChaudhuri, J., Chattopadhyay, S. and Banik, S.K. (2007) Physical Review E, 76, Article ID: 021125. https://doi.org/10.1103/PhysRevE.76.021125
[17]	Fan, W.Y., Bai, Z.W., Zhang, W. and Ding, L.P. (2020) Physica A, 547, Article ID: 123836. https://doi.org/10.1016/j.physa.2019.123836