Thermal decay rates for an asymmetric cusped barrier at strong friction

The thermal decay of a metastable state over an asymmetric cusped barrier in the regime of overdamping (strong friction) is considered. This seems to be of importance for the nanoscale experiments on pulling polymeric molecules. The decay process is simulated numerically through computer solving of the Langevin equation. The quasistationary rates RD , as well as the mean lifetimes and transient times, are extracted from the numerical time-dependent rates Rn (t). The impact of the backscattering on the value of RD is discussed. The approximate analytical decay rate is derived for the asymmetric cusped barrier. The numerical results are confronted with this formula and with another analytical formula (the integral Kramers rate) available in the literature.


Introduction
One-dimensional Brownian motion, when a fictitious particle escapes from a potential well (metastable state) over a barrier, becomes presently a widely used model. Examples of its application range from fission of atomic nuclei [1][2][3] (the femtometer typical space scale) to cracks formation [4], singlemolecule pulling [5,6], and optical traps [7] (the nanometer scale). Often in these experiments, the shape of the barrier is known only approximately or not known [8,9]. Most of the theoretical approaches are based on the seminal Kramers paper [10] where two sorts of barriers had been considered: the parabolic and cusped (edged) ones. In the present work, we focus on the cusped barrier which seems to be relevant for the single-molecule pulling experiments [8] and for the electron transfer reactions [11,12]. These barriers are schematically shown in Fig 1. The physical meaning of the coordinate is, of course, different for different physical problems. For example, in the case of nuclear fission, it is proportional to the distance between the appearing fission fragments. When one considers the molecule dissociation, is again proportional to the distance between the centers of mass of the molecule parts. For the adhesive contact problem, can be proportional to the radius of the contact area [13]. To make our consideration as problem-independent as possible, we use the dimensionless coordinate .
Through the whole paper, we use two dimensionless scaling parameters: the governing parameter = (1) and the damping parameter (2) Here is the height of the barrier separating two wells; denotes the thermal energy of the particle motion; and denote the friction and inertia parameters, respectively; is the circular frequency of the metastable state.
The purpose of the paper is to investigate quantitatively to what extent the analytical formulas for the decay rate agree with the exact numerical rates for different asymmetries of the cusped potential (see Eq. (5) below). This study is performed at = 15 which corresponds to the overdamping regime [1,14,15]. We also analyze the behavior of the mean lifetimes and transient times as functions of .

Numerical modeling
The modeling is performed for the generalized dimensionless coordinate . The corresponding numerical Euler scheme [16] for the reduced Langevin equation reads: Here = − / is a regular (driving) force; = √ / stands for the amplitude of the random force; the superscripts and + 1 indicate two consequent time moments separated by the time-step ; ( ) are the random numbers possessing a normal distribution with zero average and the variance equal to 2.
The asymmetric cusped potential ( ) reads We keep the stiffness of the left parabola fixed and vary the ratio (the asymmetry of the barrier) and, consequently, the stiffness of the right parabola. For three different values of , the cusped potential energies are presented in Fig. 1. Only for the case of symmetric potential, the position of is shown in the figure. As becomes smaller, the sink moves further from the barrier and the potential at > becomes less steep. These two circumstances are expected to enhance rescattering [17] and thus to reduce the rate.
For each set of the parameters and , our numerical simulations result in ≅ 10 5 trajectories, each trajectory terminates at = provided it does not reach , or earlier otherwise ( is the input parameter). The time-dependent numerical decay rate is calculated as follows: Here ( ) is the number of Brownian particles reaching by the time moment ; ∆ denotes the number of the particles arriving at the sink during the time interval ∆ .
Examples of ( ) corresponding to our initial conditions can be found in many works: in [14] for six smooth potential profiles, in [18] for the symmetric cusped potential. The behavior of ( ) is similar for all the potentials: after a transient stage, the rate reaches a quasistationary value which is the principal target of our study. This value is found according to the algorithm described in [19,20]; the error of is within 3%.

Analytical Kramers rates
In [10], Kramers obtained several analytical formulas for the quasistationary decay rate relevant for different values of the damping parameter and two different barrier shapes. Moreover, he proposed a general flux-over-population method for finding the rate. For the case of overdamping, the method can be written as follows Here is the probability flux over the point and ( ) is the quasistationary probability density. In our case according to Kramers [10], the probability density approximately obeys the stationary Smoluchowski equation The stationary flux may be written as is proportional to the equilibrium probability density. Applying Eqs. (9) and (10) Considering the flux to be coordinate-independent, one integrates the left and right sides of Eq. (11) over the coordinate and obtains: The population of the well standing in the denominator of Eq. (7) can be written as Making use of Eqs. (7), (12), and (13), one arrives at what we call integral Kramers formula for the rate: This formula was not written in [10] explicitly, although in fact it was derived there. One finds Eq. (14) in an explicit form in the later publications (see, e.g., Eq. (5.109) in the Risken's book [15], Eq. (3) in [14], and Eq. (2) in [21]). With modern software and computers, it is easy to evaluate the integral Kramers rate . However, some authors still prefer approximate estimates of integrals [8,9] for the asymmetric cusped potential. Following this manner, we estimate the integrals in Eq. (14). Extending the limit of integration to the plus infinity, we reduce the left integral to the Poisson's one and obtain The second integral in Eq. (14) is convenient to be represented as the sum of two integrals (from up to and from up to ). Each of these integrals is evaluated using the Laplace method and keeping in mind that the integrands reach their maximum values at . Thus, the second integral becomes Clearly, this rate depends only upon the dimensionless universal parameters and and scales with as has been shown for the case of the two-parabolas harmonic barrier [22]. The first multiplier in Eq. (17) is just the Arrhenius transition state rate [23,24]. Next, the rate is inverse proportional to as for the two-parabolas harmonic barrier in the overdamped regime [10]. Finally, when = 1, Eq. (17) goes over to the formula obtained by Kramers in [10] for the case of the symmetric cusped barrier.

Numerical results and analysis
The modeling is performed in the space diffusion regime ( = 15) for three values of the governing parameter ( = 1.5, 2.4, 3.0). We know that the effect of the right well (in other words, of the slope beyond the barrier, see Fig. 1) is most pronounced at ≫ 1 and completely absents in the energy diffusion regime ( ≪ 1) [17].
To deal with the dimensionless quantities, we present in Fig. 2 the quasistationary rate divided by the frequency of the left well. The decrease of the rates ( ) with the potential asymmetry is due to the backscattering effect. The larger , the softer the right parabola. As the result, the Brownian particles have more chances to return to the left well (i) due to the larger coordinate of the sink and (ii) due to the less steep descent from the barrier. The decrease of with at fixed is purely statistical and follows qualitatively Eq. (17). In Fig. 3, the analytical rates (Eq. (14)) and 0 (Eq. (17)) are compared with the numerical quasistationary rate for the values of under consideration. From the left column of panels, one gets a qualitative impression that seems to provide a better approximation for at any value of the governing parameter and at any value of the barrier asymmetry.
The quantitative comparison is provided through the fractional differences (right column of panels). Clearly, the integral Kramers rate almost always agree with within 10%. Thus, results obtained earlier for smooth barriers [14,25] hold for the cusped asymmetric barriers, too. The character of ( )-dependence seems to be universal, i.e. -independent. Therefore, one can hope to invent a correction for enabling to restore a correct value of . The approximate Kramers rate 0 agrees with significantly worse: the absolute value of the fractional difference 0 sometimes exceeds 30%.
There are two characteristics of ( ): its quasistationary value , which has been already discussed, and the transient time , i.e. the time interval during which the rate ( ) reaches half of its quasistationary value. It is natural to consider now the ( )-dependence which is presented in Fig. 4. We find the value of in two stages. As long as is known, our code finds using the mesh with the timestep 0.01 . The error of the resulting is rather large: it is of the order of this timestep. In the second stage, we make the upper limit for the ( )-array equal to 2 found in the first stage and repeat the search procedure. Thus, the error, which is equal to the half of the timestep, over becomes equal to 1%.
Coming back to Fig. 4, one notices the increase of with . This seems to be explained by the elongation of the way the Brownian particles have to overcome while going to the sink because as decreases the sink moves further from the barrier.  In the nanoscale experiments, besides the decay rates, they often are interested in the mean lifetime of the metastable state. That is why we present in Fig. 5 these numerical mean lifetimes evaluated as described in [26,27]:  with the governing parameter is of statistical nature: every now and again one uses inverse Eq. (17) for a rough estimation of .

Conclusions
In the present work, the numerical modeling of the thermal decay via an asymmetric cusped barrier in the overdamped regime has been performed. This seems to be timely, for instance, for the reactions of electron or proton transfer. The process has been modeled through the stochastic (Langevin) equation.
Altogether about 30 cases corresponding to the three values of the governing parameter ( =1.5, 2.4, and 3.0) and to the barrier asymmetry ranging from 0.1 up to 30 have been considered. The following quantities characterizing the decay process have been obtained (the relative accuracy in the brackets): the quasistationary rates (3%), transient times (1%), mean lifetimes (3%). The dependence of these quantities upon resulting from the numerical modeling has been interpreted qualitatively. Two analytical rates have been compared with . This comparison has shown that the integral Kramers rate (Eq. (14)) agrees with typically within 5% for the cases of high barrier ( =2.4 and 3.0) whereas the approximate Kramers rate 0 (derived in the present work, Eq. (17)) overestimates typically by 20%. For the low barrier situation, neither nor 0 agrees with within 10% at > 10.