Numerical calculation of the rate of homogeneous gas-liquid nucleation in a Lennard-Jones system

J. Chem. Phys. 110, 2159 (1999); https://doi.org/10.1063/1.477826 110, 2159 © 1999 American Institute of Physics. Thermally driven escape over a barrier of arbitrary shape Cite as: J. Chem. Phys. 110, 2159 (1999); https://doi.org/10.1063/1.477826 Submitted: 13 March 1998 . Accepted: 21 October 1998 . Published Online: 12 January 1999 A. N. Drozdov, and J. J. Brey ARTICLES YOU MAY BE INTERESTED IN Theory of activated rate processes: Exact solution of the Kramers problem The Journal of Chemical Physics 85, 1018 (1986); https://doi.org/10.1063/1.451844 Two novel approaches to the Kramers rate problem in the spatial diffusion regime The Journal of Chemical Physics 111, 6481 (1999); https://doi.org/10.1063/1.479945 Activated rate processes: Generalization of the Kramers–Grote–Hynes and Langer theories The Journal of Chemical Physics 97, 2422 (1992); https://doi.org/10.1063/1.463081 Thermally driven escape over a barrier of arbitrary shape A. N. Drozdova) and J. J. Brey Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado de Correos 1065, Sevilla 41080, Spain ~Received 13 March 1998; accepted 21 October 1998! The Kramers theory for the thermally activated rate of escape of a Brownian particle from a potential well is extended to a barrier of arbitrary shape. The extension is based on an approximate solution of the underlying Fokker–Planck equation in the spatial diffusion regime. With the use of the Mel’nikov–Meshkov result for the underdamped Brownian motion an overall rate expression is constructed, which interpolates the correct limiting behavior for both weak and strong friction. It generalizes in a natural way various different rate expressions that are already available in the literature for parabolic, cusped, and quartic barriers. Applications to symmetric parabolic and cusped double-well potentials show good agreement between the theory and estimates of the rates from numerical calculations. © 1999 American Institute of Physics. @S0021-9606~99!01404-X# I. INTRODUCTION Ever since the pioneering contribution of Svante Arrhenius, the problem of thermally driven escape from a metastable state has become one of the most fundamental problems in physics and chemistry.1The modern theory of activated rate processes is essentially due to Kramers,2who provided a dynamical framework for the original concepts of Arrhenius. The underlying idea of the Kramers theory is to model the escape process by the motion of a Brownian particle with mass weighted coordinate xin a potential of mean force V(x). The dynamics is governed by the following Fokker–Planck equation for the probability density P(x,v,t) of finding the particle at time tat position xwith velocity v:2 ] tP~x,v,t!5@2v ] x1V8~x! ] v1 g ] v~v1 b 21 ] v!#P~x,v,t!. ~1.1! Here the prime denotes the derivative with respect to x, g is the friction coefficient, and b the inverse energy available from the thermal bath, b 215kBT. The potential is assumed to have a well with minimum at xw,0, separated from the continuum by a barrier at x50 of height E52V(xw). Hereby we set for convenience V(0)50. The quantity of interest is the escape rate Gof the particle from the well. The latter can always be written in the form G5 m GTST ,~1.2! where GTST is the transition state theory ~TST!result GTST5 H A 2 p b E 2` 0dx e2 b V~x! J 21,~1.3! and m is a transmission coefficient describing the deviation of the rate from GTST . Kramers studied the dependence of the escape rate on the frictional damping in two regimes, namely, for small and intermediate to large friction g . In the former regime, the coupling between the system and the bath is assumed to be vanishingly weak so that the rate limiting step is the transfer of energy from the bath to the particle. The transmission coefficient takes in this case the form m ~ g →0!5D~ g →0!52 gb E xp 0dx A 22V~x!,~1.4! where Dis the dimensionless loss of energy per oscillation of a particle with energy close to the barrier height, and xpthe left-hand side turning point of the asymptotic underdamped trajectory, V(xp)50. In the intermediate to large friction regime, when the transfer of energy becomes fast enough to maintain thermal equilibrium of escaping particles, the rate limiting step is spatial diffusion across the barrier region. One of the basic assumptions of the Kramers theory2in this regime is a parabolic barrier approximation. It consists in dividing the full potential into a parabolic barrier part U~x!52 1 2 v 2x2,~1.5! with v 252V9(0), and an anharmonic correction reading V~x!5U~x!1O~x3!.~1.6! In the immediate vicinity of the barrier top which dominates the dynamics, the nonlinearity of V(x) vanishes faster than the parabolic part 21 2 v 2x2and therefore can be neglected. This yields the following expression for the transmission coefficient: m pb5 A 11 g 2 4 v 22 g 2 v .~1.7! It should be noted that Eq. ~1.7!is valid for g * v /(2 p b E). Consequently, in the extreme high barrier ~low temperature!limit, b E→`, one will ultimately almost always be in the spatial diffusion regime. Kramers’ model, although simple, is of wide-ranging significance to a detailed understanding and evaluating the influence of the medium on reaction rates. It has found various generalizations to the full friction range,3non-Markovian activated rate processes,4,5 multidimensional systems,6and cases without detailed balance7~for a review see Ref. 1!.In a! Permanent address: Institute for High Temperatures, 13/19 Izhorskaya Street, 127412 Moscow, Russia. JOURNAL OF CHEMICAL PHYSICS VOLUME 110, NUMBER 4 22 JANUARY 1999 21590021-9606/99/110(4)/2159/5/$15.00 © 1999 American Institute of Physics all these investigations the barrier is assumed to be parabolic, though this assumption is not always met in real physical and chemical barrier crossing processes. For example, the barrier of charge transfer reactions is often of a cusp-shaped form.8 Kramers also derived the transmission coefficient for a cuspshaped barrier,2U(x)52a u x u , m cusp~ g →`!5a g A 1 2 p b ,~1.8! but his expression is valid only in the asymptotic limit of large friction where Eq. ~1.1!can be approximated by a Smoluchowski equation. There are various attempts in the literature to bridge the strong friction limit result for a cuspshaped barrier with the TST value, m TST51, at zero damping.9,10 An analogous interpolating formula is known for a quartic barrier.1Only very recently, Berezhkovskii et al.11 have extended this formula to an arbitrary nonparabolic barrier of the form U~x!52~a/ a ! u x u a .~1.9! Their generalization reads11 m a 5 H E 2` `dx exp@ b U~x!# J 21 3 E 2` `dx exp H b F U~x!21 2 g 2x2 G J .~1.10! There is, however, a certain irony here; the above formula agrees with the known escape rates for nonparabolic barriers, but fails to reproduce the exact result for a parabolic barrier. In the latter case, it yields instead of Eq. ~1.7!an approximate expression m 25~11 g 2/ v 2!21/2.~1.11! The aim of this paper is twofold. First, we want to present an approximate rate formula, which indeed is valid for arbitrarily shaped barriers and interpolates between the limits of small and large friction. And second, we wish to compare this formula with exact numerical rates in different types of potentials. II. INTERPOLATING FORMULA To begin with we consider the spatial diffusion regime. Our purpose is to derive an approximate solution of the Fokker–Planck equation which would allow one to recover Eqs. ~1.7!and ~1.10!. This goal can be achieved in many different ways.12 Here we employ the flux over population method developed by Kramers.2Within its scope, the escape rate is defined as the ratio of a stationary diffusion current at the top of the barrier to the population of the well. Accordingly, we have to look for a current carrying stationary probability density P(x,v), that smoothly matches the equilibrium distribution Peq~x,v!5exp@2 b V~x!21 2 b v2#~2.1! in the well and vanishes beyond the barrier. The two stationary densities are related by a form function j (x,v), P~x,v!5 j ~x,v!Peq~x,v!,~2.2! which is determined from $ 2v ] x1@V8~x!2 g v# ] v1 gb 21 ] vv 2 % j ~x,v!50. ~2.3! Once the form function is known, the reactive flux formula yields for the transmission coefficient m 5 b E 2` `dvv j ~ 0,v!exp S 21 2 b v2 D .~2.4! Following Kramers, we approximate the potential V(x) entering Eq. ~2.3!by its barrier part U(x). The latter is not necessarily parabolic, it may be a sum of arbitrary ~parabolic and nonparabolic!terms U~x!52 1 2 v 2x22a a u x u a 2¯.~2.5! Moreover, we assume that j (x,v) is a function of some linear combination of xand v, j ~x,v!5 j ~%!,%5cx1bv.~2.6! Then, it is not difficult to check by direct substitution that in leading order in %and ( b E)21an approximate solution to Eq. ~2.3!reads j ~x,v!5Z21 E % `dy e b U~y!,~2.7! with %5 A v /~ gm pb!@x2~ m pb / v !v#.~2.8! In the above m pb is given by Eq. ~1.7!, while the normalization constant Zis defined by the requirement that the form function j (x,v) approaches unity in the initial well and zero in the product side. This immediately yields Z5 E 2` `dy e b U~y!.~2.9! It will be recalled here that the barrier ~temperature!is assumed to be high ~low!enough so that the potential can be well approximated by its local behavior in the vicinity of the barrier top. Otherwise one can use in Eqs. ~2.7!and ~2.9! instead of the barrier part U(x) the full potential V(x) itself. In such a case, the integration has to be restricted to the barrier region with a lower limit at, say, xwand the upper limit at a value beyond the barrier from where the recrossing probability of a particle with zero initial velocity can safely be neglected. Inserting Eq. ~2.7!into Eq. ~2.4!, we obtain the following expression for the transmission coefficient: m ab5Z21 E 2` `dx exp H b F U~x!21 2~ g v / m pb!x2 G J . ~2.10! It is a simple matter to check that for a parabolic barrier the above formula coincides with the exact Kramers result, Eq. ~1.7!, while for a purely nonparabolic barrier ( v 50) it reproduces Eq. ~1.10!. One may also note that it agrees in the limiting case of high friction with the transmission factor for an arbitrarily shaped barrier following from the corresponding Smoluchowski equation13 2160 J. Chem. Phys., Vol. 110, No. 4, 22 January 1999 A. N. Drozdov and J. J. Brey m ~ g →`!5 H g A b 2 p E 2` `dx e b U~x! J 21 ,~2.11! and reduces to unity at zero damping. A rate expression valid in the full damping range can be obtained by making use of an elegant approach developed by Mel’nikov and Meshkov.3This gives in a straightforward way m 5 m abA~D!,~2.12! with A~D!5exp S 1 p E 0 `dx ln $ 12exp@2D~x211 4!# % x211 4 D , ~2.13! where Dis given by Eq. ~1.4!. It should be noted that the ansatz of writing a uniform formula for nonparabolic barriers as a product of a spatial diffusion expression and the depopulation factor Ais ad hoc. It follows neither from Mel’nikov and Meshkov nor from Pollak, Grabert, and Ha ¨nggi turnover theories. It is our aim here to prove the utility of Eq. ~2.12! by comparing with exact numerical rates. The latter is not so obvious as one might think. Specifically, Mel’nikov and Meshkov derived the depopulation factor ~2.13!under the assumption that the escape dynamics can be described by a probabilistic integral equation in energy-action variables, whose Green function corresponds to the barrier trajectory. For a smooth potential the trajectory that leaves the barrier with the entire energy close to zero returns to it after time T→`. This infinite time, however, is no longer true for a cusped barrier where the time is of the order of the period of particle oscillation in the well. Thus the interesting issue we shall address in our numerical applications is as follows: Does the finite period of the barrier trajectory spoil the applicability of Eq. ~2.12!? III. NUMERICAL RESULTS The aim of this section is to present exact numerical rates for different types of potential barriers that would allow one to test analytical predictions. One might, at first, believe that this issue should have been settled long ago, mainly because of its continuous importance in many problems of chemical physics. To the best of our knowledge, however, there are no numerical solutions of such a type, other than those obtained in Refs. 10 and 11 under the assumption that the potential consists only of a barrier part. This assumption results in a monotonic dependence of the transmission coefficient on g ; the coefficient increases with decreasing g and reaches its maximal value at zero damping, when there is no coupling between the system and the bath. It is clear that the data so obtained are not suited for testing analytical predictions in the most problematic intermediate and weak damping regimes. Here we deal with activated rate processes in a symmetric double-well potential of the form V~x!5E 112a@x424a u x u 22~12a!x2#,a.2 1 2. ~3.1! Its barrier part varies with the parameter afrom a purely parabolic (a50) to a purely cusped (a>1) barrier, see Fig. 1. Accordingly, the frequency v entering our rate expression reads v 25 H 4~12a!E/~112a!21 2,a<1, 0a.1. ~3.2! The method used to numerically solve Eq. ~1.1!will be described elsewhere.14,15 Table I shows a list of the first nonzero eigenvalue in the considered potential for b E510 and a50, 0.5, and 1. The calculation is performed over a large range of g which covers all regimes of chemical interest, from the underdamped Brownian motion to the spatial diffusion regime. Before testing the validity of the present rate expression, we note that Eq. ~2.12!gives the transmission coefficient for the escape from a metastable state. Using the approach suggested by Mel’nikov and Meshkov,3the coefficient for a symmetric double well can be written as m 5 m abA2~D!/A~2D!.~3.3! FIG. 1. Different shapes of the potential V(x), Eq. ~3.1!, for a50~the dashed line!and a51~the solid line!. TABLE I. First nonzero eigenvalue for symmetric double-well potentials, Eq. ~3.1!with b E510 and a50, 0.5, and 1. Exponential notation 2k means that the number preceding is to be multiplied by 102k. g a50a50.5 a51 0.05 0.17124 0.14624 0.14424 0.1 0.30424 0.25924 0.24724 0.25 0.59324 0.49424 0.45624 0.5 0.86824 0.71324 0.64024 1 0.10623 0.88924 0.79124 2 0.10623 0.95224 0.85624 5 0.85824 0.91424 0.85024 10 0.60724 0.78624 0.77024 20 0.36124 0.55724 0.56424 50 0.15424 0.26824 0.28624 100 0.78025 0.13324 0.14524 1000 0.78326a0.13425a0.14725a aExact estimate of the eigenvalue calculated from the respective Smoluchowski equation. 2161J. Chem. Phys., Vol. 110, No. 4, 22 January 1999 A. N. Drozdov and J. J. Brey The least nonvanishing eigenvalue of the corresponding Fokker–Planck operator is then given by twice the rate defined by Eq. ~1.2!. The numerical values of the transmission coefficient extracted in this way are exhibited in Fig. 2, together with the analytical predictions obtained in terms of Eq. ~3.3!. As evidenced by Fig. 2, the approximate rate expression gives an upper bound to the exact result for the rate in the parabolic double-well potential. For the cusped potentials the theory overestimates the rate in both limits of weak and strong friction and underestimates it in the intermediate friction region. It is also seen that for all values of athe best agreement is achieved in the strong damping limit ( g *100). With decreasing g the error made by the ansatz ~3.3! increases and reaches maximal values in the weak damping region ( g &0.1). The theoretical expression overestimates the rate in this region by 14% for a parabolic barrier (a 50) and by 18% for a purely cusped barrier. It should be pointed out that the same is true for the turnover theory of Pollak, Grabert, and Ha ¨nggi.5As we have shown in recent papers,15,18 their theory also considerably overestimates the rate in the weak friction regime. Finally, to conclude this section we note that the barrier frequency v appearing in Eq. ~2.10!may still be left even if the barrier is purely nonparabolic. In such a case, it should be treated as a variational parameter.16 Yet another way to improve the rate formula is to take into account finite-barrier corrections. These are obtainable systematically in both regimes of weak17 and intermediate to strong friction.12 A further improvement of the overall rate expression can be achieved by using in Eq. ~2.10!a properly determined energy loss of the deterministic particle dynamics. In contrast to the weak friction expression for Dproposed by Mel’nikov and Meshkov, Eq. ~1.4!,3as well that suggested by Pollak, Grabert, and Ha ¨nggi5in their turnover theory, the deterministic approach to this quantity yields an approximation which remains correct in the full damping range, regardless of the particular shape of the potential barrier.15,18 IV. CONCLUDING REMARKS In this paper, an approximate formula for the rate of escape over an arbitrarily shaped barrier has been constructed by means of the flux over population method and the approach by Mel’nikov and Meshkov. The resulting expression agrees in the limiting case of high friction with the rate following from the corresponding Smoluchowski equation and, in the extremely underdamped regime with the rate obtained by Kramers from a diffusion equation in energy ~action!variables. It generalizes in a natural way the known rate formulas for parabolic and nonparabolic barriers. Besides, we have presented for the first time numerically exact rate constants for potentials with different barrier shapes in all regimes of chemical interest, from underdamped to overdamped Brownian motion. These results proFIG. 2. Transmission coefficient and percentage error, 1003~approximate2exact!/exact, made in m by using Eq. ~3.3!. Exact numerical results are shown by circles. ~a!a50; ~b!a50.5; ~c!a51. 2162 J. Chem. Phys., Vol. 110, No. 4, 22 January 1999 A. N. Drozdov and J. J. Brey vides the necessary foundation for testing various different rate expressions that already exist in the literature. Comparison with the numerical data shows that the present overall rate expression is rather accurate in the strong damping limit, underestimates the rate by ;0%–18% in the intermediate friction region and overestimates the rate by ;14%–23% in the weak damping regime. ACKNOWLEDGMENTS One of us ~A.N.D.!is grateful to A. M. Berezhkovskii, P. Talkner, and V. Yu. Zitserman for many helpful discussions. We acknowledge the support of the Direccio ´n General de Investigacio ´n Cientı ´ficayTe ´ cnica of Spain for financial support ~A.N.D.!and for Project No. PB95-534 ~J.J.B.!. 1P. Ha ¨nggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 ~1990!. 2H. Kramers, Physica ~Utrecht!7, 284 ~1940!. 3V. I. Mel’nikov and S. V. Meshkov, J. Chem. Phys. 85, 1018 ~1986!. 4R. F. Grote and J. T. Hynes, J. Chem. Phys. 73, 2715 ~1980!;P.Ha ¨ nggi and F. Mojtabai, Phys. Rev. A 26, 1168 ~1982!. 5E. Pollak, H. Grabert, and P. Ha ¨nggi, J. Chem. Phys. 91, 4073 ~1989!. 6H. C. Brinkman, Physica ~Utrecht!22, 149 ~1956!; R. Landauer and J. A. Swanson, Phys. Rev. 121, 1668 ~1961!; J. S. Langer, Ann. Phys. ~N.Y.! 54, 258 ~1969!. 7P. Talkner, Z. Phys. B 68, 201 ~1987!; A. N. Drozdov, Physica A 187, 329 ~1992!. 8R. A. Marcus, Annu. Rev. Phys. Chem. 15, 155 ~1963!; L. D. Zusman, Chem. Phys. 49, 295 ~1980!. 9B. J. Matkowsky and Z. Schuss, SIAM ~Soc. Ind. Appl. Math.!J. Appl. Math. 33, 365 ~1977!; D. F. Calef and P. G. Wolynes, J. Phys. Chem. 87, 3387 ~1983!; H. Dekker, Physica A 136, 124 ~1986!; E. Pollak, J. Chem. Phys. 93, 1116 ~1990!. 10A. Starobinets, I. Rips, and E. Pollak, J. Chem. Phys. 104, 6547 ~1996!. 11A. M. Berezhkovskii, P. Talkner, J. Emmerich, and V. Yu. Zitserman, J. Chem. Phys. 105, 10890 ~1996!. 12P. Talkner, Chem. Phys. 180, 199 ~1994!. 13H. Risken, The Fokker-Planck Equation, Methods of Solution and Applications ~Springer, New York, 1989!. 14A. N. Drozdov and J. J. Brey, Phys. Rev. E 57, 1284 ~1998!. 15A. N. Drozdov and P. Talkner, J. Chem. Phys. 109, 2080 ~1998!. 16P. Talkner and E. Pollak, Phys. Rev. E 50, 2646 ~1994!. 17V. I. Mel’nikov, Phys. Rev. E 48, 3271 ~1993!. 18A. N. Drozdov and J. J. Brey, Chem. Phys. 235, 147 ~1998!. 2163J. Chem. Phys., Vol. 110, No. 4, 22 January 1999 A. N. Drozdov and J. J. Brey