A Simple Model for the Calculation of Diffusion Coefficient in a Periodic Potential

Chunhui Zhang, Zhanwu Bai^*

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

DOI: 10.4236/jmp.2021.122009                

Abstract

Motivated by developing a simple model to calculate the diffusion coefficient in moderate friction region, a simplified model is proposed to deal with the diffusion of Brownian particles in a periodic potential. Where the internal noise is a Gaussian white noise, and the basic cell of the periodic potential is composed of a parabolic potential linked with a harmonic potential. When the particles cross the joint point of the potential, a time coarse-graining scheme is used to obtain a simple analytical expression of the probability distribution. The particles drift and diffuse from the first barrier to the second barrier, the passing probability over the second barrier corresponding to the escape rate becomes decrease serves as the long-jump probability. The theoretical result is confirmed by numerical simulation results. The approach can be extended to color noise case.

Keywords

Diffusion Coefficient, Periodic Potential, Brownian Motion

1. Introduction

The diffusion of Brownian particles in a spatially periodic potential is a topic of great interest in many scientific areas of physics, chemistry, and biology [1] [2] [3]. Much effort has been devoted to the study of Brownian motion in periodic potentials. The diffusion coefficient has been investigated through numerical, simulation, and analytic approaches.

The matrix-continued-fraction method was employed to investigate the Brownian motion in one and two-dimensional periodic potentials. The diffusion coefficient was obtained through numerical calculation of the dynamic structure factor [4] [5] [6] [7]. Some characteristics of the diffusion coefficient were found, such as the resonant diffusion in one-dimensional periodic potential due to the interplay between two oscillatory motions [4], the anomalous dependence of the diffusion coefficient on the friction $D\propto {\gamma }^{-\sigma }$ with $\sigma <1$ in a coupled two-dimensional potential [5]. It is found that the coupling between two degrees of freedom always reduces the multiple-jump probability and then lowers the diffusion coefficient [6]. The diffusion-path approximation and quasi-2D approximation were examined by numerical results [7]. The former strongly overestimates the diffusion coefficient at large couplings, and the latter always gives rather good results.

Many diffusion features of Brownian particles in a periodic potential have been revealed by numerical simulations. For nonseparable and anisotropic potentials, molecular dynamics simulation found that the diffusion coefficient presents different dependence on friction in low friction regime as compared with separable potentials, which is directly related to the occurrence of long jumps [8]. For two-dimensional periodic or random potentials, superdiffusion, large-step diffusion, normal diffusion, and subdiffusion were observed through Langevin simulations [9]. These rich varieties of behaviors emerge naturally from an ordinary Langevin equation for a system described by ordinary canonical Maxwell-Boltzmann statistics, without injecting special assumptions such as Levy flights or special memory effects into models of surface diffusion. The Langevin simulation results show that the diffusion coefficient behaves as $D\propto {\gamma }^{-\sigma }$ with $0<\sigma <1/3$ in a two-dimensional periodic potential due to the coupling between the x and y degrees of freedom [10].

Some analytical approaches have been developed to study the diffusion of Brownian particles in a periodic potential. By expanding the distribution function into suitable eigenfunctions, a general method was given in Ref. [11] to calculate the distribution and correlation function of the diffusive motion of particles in a one-dimensional periodic potential. The one-dimensional diffusion in potentials which have a finite number of jumps in their value and in their derivative was investigated. The jump conditions of the eigenfunctions of the corresponding Fokker-Planck-operator were derived and applied to a periodic potential [12]. The modified PGH theory [13] is applied to the motion of a particle moving on a periodic potential influenced by friction and Gaussian thermal noise [14], a uniform expression for the diffusion coefficient valid for any friction value was derived, and the finite barrier corrections were also taken into account. A semiclassical theory for the diffusion of a particle moving on a periodic potential was presented in Ref. [15]. The analytical expressions for the diffusion coefficient and and hopping length distribution are valid for memory friction and any value of friction. Two kinds of approximate schemes, the quasi-2D approximation and the effective potential approach were employed to calculate the two-dimensional diffusion rate constant of a particle driven by a white or colored noise [16]. The theoretical result is qualitatively in agreement with the numerical result. Kramers theory was used to derive simple expressions for the hopping distribution in multidimensional activated surface diffusion [17]. The derived expressions are valid on condition that the average energy loss of the particle as it goes from one barrier to the next is of the order of $k_{B}T$ or more.

Although some analytical methods have been developed, simple and exact method is still deserve explored. Combine the physical picture of diffusion and the random walk model, a model to calculate the diffusion coefficient in the turnover region of damping is proposed in the present work. The proposed method provides a simple and exact approach to calculate the diffusion coefficient. Based on this approach, the diffusion of Brownian particles in the usual cosine periodic potential can be deal with by resort to the perturbation theory. The theoretical results for the applied periodic potential are confirmed by Langevin simulation results.

2. A Simplified Model for Calculation of Diffusion Coefficient

We consider a Brownian particle moving in a periodic potential with a basic cell composed of a parabolic potential barrier linked smoothly with a harmonic potential well, which is subjected to a Gaussian white noise. The equation of motion of the particle reads

$\stackrel{\dot{}}{x}=v,\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\stackrel{\dot{}}{v}=-\gamma v-\frac{1}{m}\frac{\partial V}{\partial x}+\frac{1}{m}\xi \left(t\right),$ (1)

where m is the mass of the Brownian particle, $\gamma$ is the damping coefficient, and $V\left(x\right)$ is the periodic potential, its basic cell is given by

$V\left(x\right)=\{\begin{array}{ll}{V}_b-\frac{1}{2}m{\omega }_b^2{x}^2,\hfill & \text{regionI};\hfill \\ \frac{1}{2}m{\omega }_0^2{\left(x-2\right)}^2,\hfill & \text{regionII};\hfill \\ V_b-\frac{1}{2}m{\omega }_b^2{\left(x-4\right)}^2,\hfill & \text{regionIII}.\hfill \end{array}$ (2)

The three parts in a basic cell of the piecewise potential are connected smoothly at $x=x_{c1},x_{c2}$ (Figure 1), $V_b$ is the height of the potential barrier, $x_{c1}=1$, $x_{c2}=3$, $V_b=1$, ${\omega }_0=1$, ${\omega }_b=1$, and $m=1$ are taken in the present work. Such a potential can serve as a zero-order approximation of a cosine periodic potential. The Gaussian white noise obeys the fluctuation-dissipation theorem: $\langle \xi \left(t\right)\xi \left(t^{\prime }\right)\rangle =2m\gamma k_{B}T\delta \left(t-t^{\prime }\right)$, $k_B$ is the Boltzmann constant and T the temperature.

The probability density function in every potential barrier or potential well can be obtained exactly. However, the exact form of probability density after several step jumps is a high dimension integration due to the particle passes through the joint points with stochastic times. The problem is in essence a complex nonlinear one. To simplify the calculation, we construct a time coarse-grain model: the particle passes the joint point $x_{c1}$ with mean passage time of the corresponding potential barrier region. The initial velocity of the particle starting diffusion from the barrier top is set to the average velocity calculated by the Kramers formula. We label the first potential barrier, the first potential well, and the second potential barrier as region I, II, and III, respectively, as shown in Figure 1. The transition probability density and probability density in the region I for an initial $\delta$ distribution of the probability density reads [18]

$W_1\left(x\mathrm{,}v\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{x}_0\mathrm{,}{v}_0\mathrm{,0}\right)=N_1\exp \left[-{\alpha }_1{x}^2-{\beta }_1{v}^2-{\gamma }_{1}xv\right]\mathrm{,}$ (3)

where

${\alpha }_1=1/2{\sigma }_{11}^{-1}\left(t\right),{\beta }_1=1/2{\sigma }_{22}^{-1}\left(t\right),{\gamma }_1={\sigma }_{12}^{-1}\left(t\right),N_1=\frac{1}{2\pi }{\left[4{\alpha }_1{\beta }_1-{\gamma }_1^2\right]}^{1/2}.$ (4)

The expressions of the second order moments ${\sigma }_{ij}$ can be found in Ref. [18]. The equivalent probability density of the particle at $x_{c1}$ is obtained by concentrating all probabilities the particle appearing in region I with $x>0$ on this point at the mean first passage time $t_1$ according to the coarse-grain approximation, that is

$W_{01}\left(x\mathrm{,}v\mathrm{,}{t}_1\right)=N_{20}\exp \left[-b_1{\left(v-{\bar{v}}_1\right)}^2\right]\delta \left(x-x_{c1}\right)\mathrm{,}$ (5)

with $b_1={\beta }_1,{\bar{v}}_1=-{\gamma }_1{x}_{c1}/\left(2{\beta }_1\right)$, and $N_{20}=1/2\sqrt{{\beta }_1/\pi }$. The mean first passage time $t_1$ is given by

$\begin{array}{c}{t}_1={\displaystyle {\int }_0^{\infty }\text{d}t}{\displaystyle {\int }_0^1\text{d}x}{\displaystyle {\int }_{\infty }\text{d}v}{W}_1\left(x,v,t;x_0,v_0,0\right)\\ ={\displaystyle {\int }_0^{\infty }\text{d}t}\frac{N_1\pi }{\sqrt{4{\alpha }_1{\beta }_1-{\gamma }_1^2}}\text{erf}\left(\sqrt{{\alpha }_1-\frac{{\gamma }_1^2}{4{\beta }_1}}{x}_{c1}\right)\end{array}$ (6)

The total probability for the particle appearing in region I for $x>0$ has been normalized as 1. The transition probability density [18] in region II is given by

$\begin{array}{l}{P}_2\left(x,v,t;x_{c1},v_1,t_1\right)\\ =N_2\exp \left(-{\alpha }_2{\left(x-{{\bar{x}}^{\prime }}_2\right)}^2-{\beta }_2{\left(v-{{\bar{v}}^{\prime }}_2\right)}^2-{\gamma }_2\left(x-{{\bar{x}}^{\prime }}_2\right)\left(v-{{\bar{v}}^{\prime }}_2\right)\right).\end{array}$ (7)

Figure 1. The potential profile. $x_{c1}\mathrm{,}{x}_{c2}\mathrm{,}\cdots$ are the joint points of parabolic barriers and harmonic potential wells.

where

$\begin{array}{l}{{\bar{x}}^{\prime }}_2=G_{11}{x}_{c1}+G_{12}{v}_1,{{\bar{v}}^{\prime }}_2=G_{21}{x}_{c1}+G_{22}{v}_1,\\ {\alpha }_2=1/2{\sigma }_{11}^{-1}\left(t\right)\mathrm{,}{\beta }_2=1/2{\sigma }_{22}^{-1}\left(t\right)\mathrm{,}{\gamma }_2={\sigma }_{12}^{-1}\left(t\right)\mathrm{,}\\ N_2=\frac{1}{2\pi }{\left[4{\alpha }_2{\beta }_2-{\gamma }_2^2\right]}^{1/2}\mathrm{,}\end{array}$ (8)

all of these quantities are defined in region II. The probability density in region II is then obtained

$\begin{array}{l}{W}_2\left(x\mathrm{,}v\mathrm{,}t\right)={\displaystyle {\int }_{-\infty }^{\infty }\text{d}{v}_1}{\displaystyle {\int }_{-\infty }^{\infty }\text{d}{x}_1}{W}_{01}\left(x_1\mathrm{,}{v}_1\mathrm{,}{t}_1\right)P\left(x\mathrm{,}v\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{x}_1\mathrm{,}{v}_1\mathrm{,}{t}_1\right)\\ =N_0{N}_2\sqrt{\frac{\pi }{D_2}}\exp \left(-a_2{\left(x-{\bar{x}}_2\right)}^2-b_2{\left(v-{\bar{v}}_2\right)}^2-c_2\left(x-{\bar{x}}_2\right)\left(v-{\bar{v}}_2\right)\right)\mathrm{,}\end{array}$ (9)

where

$\begin{array}{l}{\bar{x}}_2=G_{11}{\bar{x}}_1+G_{12}{\bar{v}}_1,{\bar{v}}_2=G_{21}{\bar{x}}_1+G_{22}{\bar{v}}_1,\\ D_2=b_1+{\alpha }_2{G}_{12}^2+{\beta }_2{G}_{22}^2+{\gamma }_2{G}_{12}{G}_{22},\\ a_2=\frac{4b_1{\alpha }_2+\left(4{\alpha }_2{\beta }_2-{\gamma }_2^2\right)G_{22}^2\left(t_2-t_1\right)}{4D_2},\\ b_2=\frac{4b_1{\beta }_2+\left(4{\alpha }_2{\beta }_2-{\gamma }_2^2\right)G_{12}^2\left(t_2-t_1\right)}{4D_2},\\ c_2=\frac{4b_1{\gamma }_2-2\left(4{\alpha }_2{\beta }_2-{\gamma }_2^2\right)G_{12}\left(t_2-t_1\right)G_{22}\left(t_2-t_1\right)}{4D_2}.\end{array}$ (10)

The transition probability density in region III is given by

$\begin{array}{l}{P}_3\left(x\mathrm{,}v\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{x}_{c2}\mathrm{,}{v}_2\mathrm{,}{t}_2\right)\\ =N_3\exp \left(-{\alpha }_3{\left(x-{{\bar{x}}^{\prime }}_3\right)}^2-{\beta }_3{\left(v-{{\bar{v}}^{\prime }}_3\right)}^2-{\gamma }_3\left(x-{{\bar{x}}^{\prime }}_3\right)\left(v-{{\bar{v}}^{\prime }}_3\right)\right)\mathrm{.}\end{array}$ (11)

Where

$\begin{array}{l}{{\bar{x}}^{\prime }}_3=G_{11}{x}_{c2}+G_{12}{v}_2,{{\bar{v}}^{\prime }}_3=G_{21}{x}_{c2}+G_{22}{v}_2,\\ {\alpha }_2=1/2{\sigma }_{11}^{-1}\left(t\right),{\beta }_2=1/2{\sigma }_{22}^{-1}\left(t\right),{\gamma }_2={\sigma }_{12}^{-1}\left(t\right),\\ N_2=1/\left(2\pi \sqrt{\det \left(\sigma \right)}\right)\mathrm{,}\end{array}$ (12)

All these quantities are defined in region III. The probability density in region III is expressed as

$W_3\left(x\mathrm{,}v\mathrm{,}t\right)={\displaystyle {\int }_{t_1}^t\text{d}{t}_2}{\displaystyle {\int }_{-\infty }^{\infty }\text{d}{v}_2}{v}_2{W}_2\left(x_{c2}\mathrm{,}{v}_2\mathrm{,}{t}_2\right)P_3\left(x\mathrm{,}v\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{x}_{c2}\mathrm{,}{v}_2\mathrm{,}{t}_2\right)\mathrm{.}$ (13)

The escape probability at time t crossing over the second barrier top is

$P_e\left(t\right)={\displaystyle {\int }_0^{\infty }\text{d}x}{\displaystyle {\int }_{-\infty }^{\infty }\text{d}v}{W}_3\left(x\mathrm{,}v\mathrm{,}t\right)\mathrm{,}$ (14)

Performing the Gaussian integrations over x and v, the escape probability can be expressed as

$\begin{array}{l}{P}_e\left(t\right)={\displaystyle {\int }_{t_1}^t\text{d}{t}_2}{\displaystyle {\int }_{-\infty }^{\infty }\text{d}{v}_2}{v}_2\frac{1}{2}\sqrt{\frac{\pi }{D_2}}{N}_0{N}_2\text{erfc}\left[-\sqrt{{\alpha }_3-\frac{{\gamma }_3^2}{4{\beta }_3}}{{\bar{x}}^{\prime }}_3\right]\\ \times \exp \left[-a_2{\left(x_{c2}-{\bar{x}}_2\right)}^2-b_2{\left(v_2-{\bar{v}}_2\right)}^2-c_2\left(x_{c2}-{\bar{x}}_2\right)\left(v-{\bar{v}}_2\right)\right].\end{array}$ (15)

In random walk model, the diffusion coefficient is expressed as [19] [20]

$D=k\langle l^2\rangle =k{\displaystyle \underset{n}{\sum }}{\left(nd\right)}^2{P}_n,$ (16)

where d is the spatial periodic of the potential, $P_n$ is the probability of n-step jumps, and k is the Kramers rate calculated by the Kramers formula in spatial diffusion regime, which is still valid for damping out of spatial diffusion regime due to the parabolic potential barrier. When the center of the probability packet moves toward to the second barrier, the passing probability over the barrier top

increases rapidly at almost a constant rate $k_e=\frac{\text{d}{P}_e}{\text{d}t}$, as shown by the numerical

results. When the center of the probability packet moves back to the potential well, the passing probability over the second potential barrier only increases due to diffusion and then becomes slow. We take the critical probability $P_c$ that the increase of $P_e\left(t\right)$ as a function of time from rapid to slow as the long jump probability (more than one step). Thereafter the process in the first basic cell is repeated periodically, which is a part of our simplified model. The probabilities for n-step jumps is then given by

$P_n=P_c^{n-1}\left(1-P_c\right)\left(n=1,2,\cdots \right)$ (17)

such a geometric progression jump probability distribution is a good approximation for several step jumps, as shown by simulation results (see Fig. 6 of Ref. [8] ).

3. Comparison with Langevin Simulation Results

To check the accuracy of the diffusion coefficient obtained by the simplified model,

Figure 2. The diffusion coefficient as a function of damping. Where the potential barrier height $V_b=1$, the spatial periodic $d=4$. (a) for $T=0.2$, (b) for $T=0.3$, and (c) for $T=0.4$.

we simulate the Langevin Equations (1) by second order Runge-Kutta algorithm. The number of test particles and the time step are taken as $N=3\times {10}^5$ and $\Delta t=2\times {10}^{-3}$, respectively. As shown in Figure 2, the theoretical results match the simulation results well in moderate friction region. Such a region is common in surface diffusion problem. The maximal error for applied parameters is less than 6% until lower reduced potential barrier height $V_b/k_{B}T=2.5$. The moderate friction region is called the turnover region in escape theory, which is not covered by the original Kramers escape theory. The calculation of escape rate (a factor of diffusion coefficient, see Eq. (16)) and diffusion coefficient in this region is lack of a simple method.

4. Conclusion

A simple model is proposed to calculate the diffusion coefficient for Brownian particles moving in a periodic potential. The basic cell of the periodic potential is composed of a parabolic potential barrier linked with a harmonic potential well smoothly, which can serve as a zero-order approximation of a cosine periodic potential. Further theoretical results for the common cosine periodic potential can be obtained by perturbation theory. The theoretical result for the applied potential is confirmed by the simulation result in moderate friction region, which is an often encountered region in surface diffusion problem. The proposed approach can be generalized conveniently to color noise case.

Conflicts of Interest

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

References

[1]	Ala-Nissila, T. and Ying, S.C. (1992) Progress in Surface Science, 39, 227-323. https://doi.org/10.1016/0079-6816(92)90017-C
[2]	Gardiner, C.W. (1990) Handbook of Stochastic Methods. Springer-Verlag, Berlin.
[3]	Oh, S.-M., Koh, S.J., Kyuno, K. and Ehrlich, G. (2002) Physical Review Letters, 88, Article ID: 236102.
[4]	Caratti, G., Ferrando, R., Spadacini, R., Tommei, G.E. and Zelenskaya, I. (1998) Chemical Physics Letters, 290, 509-513. https://doi.org/10.1016/S0009-2614(98)00543-0
[5]	Caratti, G., Ferrando, R., Spadacini, R. and Tommei, G.E. (1997) Physical Review E, 55, 4810. https://doi.org/10.1103/PhysRevE.55.4810
[6]	Caratti, G., Ferrando, R., Spadacini, R. and Tommei, G.E. (1998) Chemical Physics, 235, 157-170. https://doi.org/10.1016/S0301-0104(98)00126-8
[7]	Caratti, G., Ferrando, R., Spadacini, R. and Tommei, G.E. (1996) Physical Review E, 54, 4708. https://doi.org/10.1103/PhysRevE.54.4708
[8]	Chen, L.Y., Baldan, M.R. and Ying, S.C. (1996) Physical Review B, 54, 8856. https://doi.org/10.1103/PhysRevB.54.8856
[9]	Lacasta, A.M., Sancho, J.M., Romero, A.H., Sokolov, I.M. and Lindenberg, K. (2004) Physical Review E, 70, Article ID: 051104. https://doi.org/10.1103/PhysRevE.70.051104
[10]	Braun, O.M. and Ferrando, R. (2002) Physical Review E, 65, Article ID: 061107. https://doi.org/10.1103/PhysRevE.65.061107
[11]	Risken, H. and Vollmer, H.D. (1978) Zeitschrift für Physik B, 31, 209-216. https://doi.org/10.1007/BF01333894
[12]	Mrsch, M., Risken, H. and Vollmer, H.D. (1979) Zeitschrift für Physik B, 32, 245-252. https://doi.org/10.1007/BF01320120
[13]	Pollak, E. and Ianconescu, R. (2016) The Journal of Physical Chemistry A, 120, 3155-3164. https://doi.org/10.1021/acs.jpca.5b11502
[14]	Ianconescuab, R. and Pollak, E. (2016) Faraday Discuss, 195, 111-138. https://doi.org/10.1039/C6FD00105J
[15]	Georgievskii, Y. and Pollak, E. (1994) Physical Review E, 49, 5098. https://doi.org/10.1103/PhysRevE.49.5098
[16]	Bao, J.D. (2001) Physical Review E, 63, Article ID: 061112. https://doi.org/10.1103/PhysRevE.63.061112
[17]	Hershkovitz, E., Talkner, P., Pollak, E. and Georgievskii, Y. (1999) Surface Science, 421, 73-88. https://doi.org/10.1016/S0039-6028(98)00820-6
[18]	Risken, H. (1989) The Fokker-Planck Equation. Springer, Berlin. https://doi.org/10.1007/978-3-642-61544-3
[19]	Chudley, C.T. and Elliott, R.J. (1961) Proceedings of the Physical Society of London, 77, 353. https://doi.org/10.1088/0370-1328/77/2/319
[20]	Braun, O.M. and Sholl, C.A. (1998) Physical Review B, 58, 14870-14879. https://doi.org/10.1103/PhysRevB.58.14870