Langevin Analysis of the Diffusion Model for Surface Chemical Reactions

An analysis is presented of the magnitude of some of the potential sources of error in a Iecently developed diffusion model of surface chemical reactions. Using single absorber Langevin simulations. comparisons are made between the diffusion equation model and the Fokker-PIanck equation for the rates of diffusion controlled surface chemical reactions. The diffusion equation is found to predict rates in good agreement with the Fokker-PIanck equation for physical values of the diffusion constant. For unphysica1Iy large diffusion constants, the rates predicted by the diffusion equation are found to be in error. By employing multiple absorber Langevin simulations errors in the single absorber approximation used in the diffusion model of surface reactions are examined. The single absorber model is found to be accurate for weakly bound adsorbates. For strongly bound adsorbates rate expressions derived from a two-dimensional model are found to be appropriate. The relative rates of Eley-Ridea1 and Langmuir-Hinshelwood mechanisms are also studied by multiple absorber Langevin simulations. The ratios of the Eley-Ridea1 to the Langmuir-Hinshelwood rate is found to be in good agreement with the predictions of the diffusion equation model for physica1Iy reasonable diffusion constants. The time dependent solution to the diffusion equation considered in a previous publication is given in an appendix.


I. INTRODUCTION
In recent years, diffusion constants for the migration of adsorbates on crystal surfaces have become available both from molecular dynamics calculations 1-4 and a variety of experimental techniques.5 The interest in surface diffusion is partially motivated by the fact that diffusion is a primary step in the dynamics of many surface processes.To connect diffusion information with kinetic rate constants in a recent publication 6  [hereafter referred to as I] we extended the theory for diffusion controlled reactions in solutions 7 to surface diffusion influenced reactions.In I, we analyzed the kinetic implications of the two-dimensional diffusion equation aW(r, t) _ DV2Wi( t) _ Wi(r, t) J at r, r which was first used for thin film nucleation studies.8 In Eq. (1), Wi(r, t) was the concentration of reactant species at coordinate r and time t, D was the sum of the diffusion constants for reactant molecules, r was the lifetime of a reactant molecule before desorption, and J was the number of reactant molecules per unit time per unit area adsorbing on the surface externally from the gas phase.Equation (1) differed from the usual diffusion equation by the addition of terms which allowed reactant molecules to enter or leave the reaction surface.
As pointed out in I, Eq. (1) has two parameters, Jr and where ')1.1 can be interpreted as one-half the average "'Visiting staff Member at the Los Alamos National Laboratory.Permanent address: Department of Chemistry, The University of Rhode Island, Kingston, Rhode Island 02881.
distance traveled by a nonreacting molecule before desorption.In terms of the key parameters from Eq.
(1), expressions were derived in I for diffusion influenced rate constants, the activation energies for diffusion influenced reactions, and the ratio of the rates of Eley-Rideal processes to Langmuir-Hinshelwood processes.In I, we showed within the model defined by Eq. ( 1) that the Eley-Rideal rate should be less important than the Langmuir-Hinshelwood rate for many reactions, and we showed the activation energy for Langmuir-HinsheLwood reactions to be bounded by the activation energy for diffusion and the activation energy for recombination.As pointed out in I, the expressions developed for diffusion influenced rate constants in terms of surface diffusion constants are particularly valuable for theoretical study because diffusion constants can be calculated from a portion of the potential energy surface required for full reaction dynamics calculations.
The model developed in I for diffusion influenced rate expressions for surface reactions contains a number of assumptions the effect of which require further analysis.The expressions developed in I are based on the behavior of a two-dimensional isotropic fluid which obeys the diffusion equation.Actual surface reactions occur on a lattice with appreciable activation barriers.The effect of such a lattice may be diffusion constants with spatial anisotropies.In I, isotropic diffusion constants were assumed.The rate constants within the diffusion model are evaluated from expreSSions for the concentration at an absorption boundary.It is well-known 9 that the diffusion equation gives inaccurate concentration profiles at such absorbing boundaries, and that a careful treatment requires solutions to the Fokker-Planck equation for the full phase space distribution function.In I by analogy with diffusion influenced reactions in solution the rate expressions were developed from the limiting form of a single absorber model.In contrast to solution kinetics the single absorber model may be inappropriate for surface reactions, because the influence of the absorption boundary conditions on the concentration profile is long range in two dimensions.
In the present work, we analyze the magnitude of the errors in I which arise from the use of the diffusion equation within the Single absorber approximation.We carry out this analysis by comparing the expressions derived in I with the results of simulations of diffusive processes using the Langevin equation.We show the diffusion equation model to be accurate when the magnitude of the diffusion constant is on the order of those measured and calculated for real physical systems.Only for nonphysically large diffusion constants do we find errors in the use of the diffusion equation to be appreciable.
The organization and contents of the remainder of this paper are as follows.In Sec.II, we compare the rate constants evaluated from the diffusion equation with the results of Langevin simulations within the single absorber model.To justify the parameters used in the simulations, we present the time dependent solution to Eq. (1) which extends the steady-state solution we gave previously.In Sec.II, we compare both rate constants and concentration profiles from the Langevin simulations with the diffusion equation.In Sec.III we compare the results obtained from the single absorber model with multiple absorber simulations.We determine the kinds of parameters for which the results of I are valid, and indicate expressions from I which can be used in other cases.We also calculate the ratio of Eley-Rideal rates to Langmuir-Hinshelwood rates and compare with the expression derived in I for this ratio.In Sec.IV, we summarize our conclusions.We derive the time dependent solution to Eq. (1) in the Appendix and develop the behavior of the solution in a number of limits.

II. SINGLE ABSORBER LANGEVIN SIMULATIONS
The migration of adsorbated on crystal surfaces is accurately described by the classical diffusion equation only on time scales which are long compared to the inverse of the effective friction constant. 9For time scales on the order of the inverse of the friction constant an accurate treatment of adsorbate migration requires solutions to the Fokker-Planck equation for the full phase space distribution function.For problems with absorbing boundary conditions the inaccuracies in the diffusion equation extend to distances close to the absorbing region at all times.For time scales on the order of the time for molecular motion the FOkker-Planck equation is inaccurate and surface migration must be described by the generalized Langevin equation. 10,l1   For diffusion controlled surface reactions it is imagined that the rate limiting step is the time necessary for reactants to diffuse to some critical reaction distance.Once reactants reach a critical distance the reaction is assumed to be very rapid.In I, we modeled such diffusion controlled reactions with the diffusion equation given in Eq. ( 1).The rate constant expres-sions derived in I required the evaluation of the .absorptionrate at the absorption boundary.This rate was evaluated at steady state.Although the time to reach steady state is long compared to the inverse of the friction constant, the evaluation of the rate at the absorption boundary may be inaccurate, because of deficiencies in the diffusion equation at short distances.To test the errors in the rate derived from the diffusion equation, in this section we present the results of numerical calculations of rates by Langevin simulations.The Langevin simulations give rates equivalent to rates which would be obtained from a solution to the Fokker-Planck equation.The Langevin simulations are carried out with geometries identical to the geometries used in I. Consequently, in this section, errors in the diffusion equation are examined within the single absorber model.Errors introduced by the use of the single absorber approximation will be examined in Sec.III.
To derive expressions for rate constants for diffusion controlled reactions in I, we solved Eq. ( 1) at steady state subject to the boundary conditions where WD(r) is the steady state solution to Eq. ( 1), RA is a critical distance within which we assume all diffusing absorbates to react with unit probability, and RB is an outer radius which provides a source of reactants at the initial concentration Co.The solution to Eq. ( 1) at steady state subject to the boundary conditions given in Eqs. ( 3) and ( 4) was found to be where In(x) and Kn(x) are modified Bessel functions of the first and second kind of order n, and A and B are coefficients whose detailed expressions are given in I [Eqs.
(49) and ( 50)].The rate constants were extracted from Eq. ( 5) from the relation In I, Eq. ( 6) was evaluated in a number of limits of which and are important to the present discussion.For finite RB and y the rate constant from Eqs. ( 5) and ( 6) is In Eqs. ( 6), ( 7), and (9) we have taken Co =JT for convenience.
Because a particle moving in a viscous medium will only display diffusive behavior at times long compared with the inverse of the friction constant it is well known that the diffusion equation provides a poor description of stochastic motion at short times and short distances. 9To describe the distribution of particles near the adsorber at radius RA for times long relative to the time scale of molecular motions it is necessary to solve the Fokker-Planck equation for the full phase space distribution function.When particles are allowed to flow into and out of the surface with the parameters of Eq. ( 1) the Fokker-Planck equation is In Eq. (10), f(r, u, t) is the phase space distribution function, m is the mass of the particles, u is the particle velocity, and f3 = (1lk B T), kB being the Boltzmann constant, and T being the absolute temperature.Direct solutions to Eq. (10) subject to absorbing boundary conditions are difficult to obtain.12 To compare the rate constants and concentrations profiles obtained from Eq. (1) with Eq. (10) we found it more convenient to solve the equivalent Langevin equation In Eq. ( 11), ~ is the friction constant related to the diffusion constant by 1 ~=f3mD and R(t) is a random force.To solve Eq. ( 11) under conditions appropriate for the present study we used the geometry shown in Fig. 1.At time t= 0, particles were scattered over the entire region at random so that the average concentration was Co and the velacities we thermalized to temperature T. For each integration, time step particles were added to region C to maintain a constant outer concentrations of Co.Reflecting boundary conditions were applied to region C at the square walls.Particles were allowed to pass freely between regions B and C. Particles entering region A were removed and counted as absorbed.To account for a desorption mechanism particles in region B were removed in each time step with probability AtIT, where At was the length of a time step.Particles were adsorbed onto region B according to a Poisson distribution with an average number of particles Jw(R~ -R~)At added in each time step.The Langevin equation was solved by methods implicit in Eq. ( 240) of Ref. 9 using Gaussian random noise for the integrated random force.
To understand the parameters used in the Langevin simulations it is useful to consider the solution to the time dependent diffusion equation [Eq.( 1 (13) The coefficients Crt and a", and the function Vo(a"r) used in Eq. ( 13 are defined at steady state we wish to solve the Langevin equation at steady state.We see that one parameter which will govern the rate of decay of the transient part of Eq. ( 13) is 'T"1.For small T the transient solution will decay quickly and a comparison with the steady state solutions will be meaningful.Consequently, in the Langevin simulations we chose very short lifetimes to desorption so that steady state comparisons could be made.is shown in Fig. 2. In Fig. 2, the solid line represents the concentration profile evaluated from the Langevin simulations and the dashed line gives the concentration evaluated with Eq. ( 5).The evaluated points are given for the Langevin case.The concentration from the diffusion equation underestimates the Fokker-Planck concentration profile in agreement with recent onedimensional studies of Harris.12 The concentration falls to zero far more abruptly in the Fokker-Planck case than the diffusion case.This does not necessarily imply a larger absorption rate because Eq. ( 6) is only valid for diffUSion equation solutions.Absorption rate constants were calculated from the Langevin simulations by evaluating the number of particles entering region A of Fig. 1 as a function of time.
2.25 X 10-2 3.42x10-3 2.03 x 10-3 1. 05 X 10-3 6.54 X 10-4 4.14x10-4 1. 54 X 10-4 2.17 X 10-2 3. 42x 10-3 2.03 X 10-3 1. 05 X 10-3 6.54 X 10-4 4.14x10-4 1. 54 X 10-4 8.99x10-4 7.11 X 10-4 7. 02X 10-4 5,74 X 10-4 4.94 X 10-4 3.04 X 10-4 1. 98 X 10-4 tained from the slope of Fig. 3 and the rate constant is defined as this rate per unit concentration.Rate constants from Langevin simulations kL are compared with kD [Eq.( 6)J and k Doo [Eq.( 7)] in Table I.The values used for R a , R B , Co, m, J, T, and 7 are the same in Table I as in Fig. 3.There is complete agreement between kD and k Doo for all values of D less than 10-2.For D = 10-2 , 1'-1 = 31.62 so that the concentration profile has not leveled off at RB = 50.Consequently, the behavior at RB = 50 is somewhat different than the limiting case of infinite R B • For D = 10-3 or less y-I:s 10 so that the steady-state solution to Eq. ( 1) at RB = 50 is the same as at infinite R B • The rate constants evaluated from the diffusion model are nearly linear in D as can be understood qualitatively from Eq. ( 7).This linear behavior is not observed from the Langevin simulations.For large values of D, kL is nearly D independent.The near D independence occurs for large D because the small friction constant leads to nearly free translational motion within the geometry of the calculation.Indeed, lim kL = const , values of D because of Eq. ( 17).As D becomes small the ratio appears to approach unity.Diffusion constants determined for physical systems tend to be on the order of 10-4 or 10-5 a. u.Consequently, the expressions for diffusion influenced rate constants developed in I can be expected to be sufficiently accurate for kinetic analysis of diffusion constant information.We can expect that basing a rate theory on the diffusion equation rather than the Fokker-Planck equation will be accurate for physical diffusion constants.

III. MULTIPLE ABSORBER SIMULATIONS
In solution kinetics, rate constants for diffusion controlled reactions are extacted from diffusion constants in the limit that RB is taken to infinity.The infinite RB limit is physically appropriate for reactions in solution, because such reactions are always dilute owing to the presence of a solvent.In addition, the influence of the outer boundary at RB is weak for solution reactions in contrast to the case of diffusion in two dimensions where the outer boundary has a strong influence on the concentration profile and calculated rate constants (see I, Sec.III A).The infinite RB limit is the essence of the single absorber model discussed in Sec.II.
To help assess the validity of the single absorber approximation (or infinite RB limit), we carried out a series of multiple absorber Langevin simulations to compare with the diffusion controlled reaction rate constants evaluated in I and in Sec.II of this work.At time t = 0, we dispersed 100 particles at random with thermalized velocities on a square two-dimensional surface of dimension 100x 100 a. u.We propagated the motions of the particles with the Langevin equation.In each time step, we remove and counted as reacted any two particles whose centers were separated by a distance R A • To simulate the absorption and desorption processes discussed in Sec.11 particles were desorbed with probability MIT in each time step of length t:.t, and an average of 10 4 J t:.t particles were added to the surface in each time step.To approximate an extended system periodic boundary conditions were imposed at the walls of the surface.
In each simulation the parameters were taken to be RA = 2, Co = O. 01, T= 300 K, m = 29166, J = 10-6 , and T = 10 5 , all in a. u.For each diffusion constant, we ran 100 trajectories except for D = 10-5 a. u., where ten trajectories were run.The decrease in the number of trajectories was necessary for small diffusion constants, because the computer time became prohibitive.For each calculation the number of reacted particles was measured as a function of time from t = 5 X 10 4 to t = 25 X 10 4 a. u.The data were then plotted an example of which is shown in Fig. 5 (D = 10-4 a. u.).Unlike the single absorber simulations (Fig. 3) no error bars are shown in Fig. 5, because the errors are smaller than the resolution of the graph.The slope of the graph of the number of reacted particles N as a function of time gives the reaction rate.To extract a rate constant from the rate it is necessary to determine the steady-state concentration of particles on the surface which in general will differ either from Co or JT.The steady state concentration was determined by mOnitoring the number of remaining surface particles at the end of each time step.
In Table II, we present the final concentration C, the multiabsorber rate constant kJlA along with kD and kL for each calculated value of D. The values of D used to evaluate kD and kL were twice that used for k JlA , because the total diffusion constant for both reactants are required in the single absorber model.It is important to recognize that a rigorous comparison between kJlA and kD or kL is not possible owing to the fact that the final concentrations are not identical.The meaning of J, T, and Co are not the same in the single and multiple absorber simulations.
in the single absorber Langevin simulations, and occurs because the particles do not exhibit any diffusive motion over distances on the order of the interparticle separations.As D decreases, kliA exhibits a dependence on D, and when D is on the order of 10-4 or ur 5 a. u. (the physical range of diffusion constants) the agreement between the rate constants is reasonably good.As the diffusion constant diminishes the rate constant decreases resulting in a higher final surface concentration.
The results presented in Table II are representative of the behavior of surface diffusion controlled reactions when T is very short.We were unable to carry out simulations for large T because the approach to steadystate behavior becomes prohibatively long for large T. The principal effect of making T large is to increase ,,-1.In a true multiparticle surface reaction the concentration profile should reach an asymptote over a distance on the order of the interparticle separation.For the multiple absorber simulations reported here that distance is given by Ci 1 / 2 • For the diffusion constants used in this study, this distance ranges between 12 and 14 a. u. which is always larger than ,,-1.Consequently, kliA and kD have been found to be in good agreement for physical values of the diffusion constant.However, for very large T, ,,-1 would become much larger than an interparticle spacing rendering the single absorber model inaccurate.A more accurate representation of the diffusion controlled rate constant would be a model with finite RB where RB is taken to be RB=Ci1/2/2.( 18) For large T [or small", see I, Eq. ( 55)] the finite RB rate constant reduces to kDo [Eq.( 8)].
The multiple absorber simulations also allowed us to calculate the ratio of the rate of Eley-Rideal reactions to the rate of Langmuir-Hinshelwood reactions.The E ley -Rideal rate was calculated by removing all reactant particles located within a distance RA of particles added by the external flux J.This calculated ratio is compared with the expression given in Eq. ( 79) of I, ( 19) where ( 20) Equation ( 19) is based on the diffusion model.The multiple absorber Langevin values of 0, OL, are given in Table III along with results evaluated from Eq. ( 19).The agreement between the multiple absorber simulation ratio and the diffusion equation ratio is good for small values of D. As expected the agreement is poor for large D, because the diffusion equation becomes inaccurate for the reasons discussed previously.

IV. CONCLUSIONS
In this work, we have carried out a critical analysis of some of the assumptions in our recent development of a two-dimensional isotropic model of diffusion influenced reactions on crystal surfaces.Using Langevin simulations we have found the single absorber model discussed in Sec.II to be well described by the diffusion equation for diffusion constants on the order of those expected for chemisorbed systems.
By analysis of multiabsorber Langevin simulations we have been able to infer that the single absorber model for surface reactions is inappropriate for strongly bound species with very long absorption lifetimes.For dilute weakly bound systems diffusion controlled rate constants can be best obtained from Eq. ( 7), whereas Eq. ( 8) is appropriate for strongly bound systems.When Eq. ( 8) is used RB should be taken from Eq. ( 18).
In I, we developed expressions for the activation energy and the ratio of the Eley-Rideal to Langmuir-Hinshelwood rates for diffusion controlled surface reactions at the infinite RB limit.If the finite RB limit is more appropriate the activation energy for the reaction is identical to the activation energy for diffusion in the diffusion controlled case [see Eq. ( 8)].In the infinite RB case, the activation energy for diffusion controlled reactions is only approximately given by the activation energy for diffusion [see I, Eqs.(87) to (89)].
For long adsorption lifetimes when RB is finite and" approaches zero the ratio of the Eley-Rideal to Langmuir-Hinshelwood rates is given by (see I, Sec.III C for details of such developments) Since T is large, 6 will be small for cases, where RB is taken to be finite, and the Langmuir-Hinshelwood mechanism will dominate.
We are presently using the results developed in I and in the present work to extract rate information from calculated diffusion constants.The results of these calculations will appear separately.
FIG.1.The geometry used in the single absorber simulations.
The spatial concentrations used initially in the Lange-Vin simulations are the same as the boundary conditions used to solve Eq. (1); namely, W(r,O)=C o , RA<rsR B , W(RA,t) =0 , t>O, W(RB,t)=C O , t>O.diffusion equation and Fokker-Planck equation concentration profiles will agree exactly at t = O.The principal discrepancies between the Fokker-Planck and diffusion equation concentration profiles will occur at long times for distances r near R A • To observe the discrepancies between the Fokker-Planck and diffusion equation concentration profiles, Langevin simulations were carried out with RA = 20, RB = 50, Co = 0.04, T = 300 K, m = 29166 [that of an oxygen atom], D = 10-3 , J = O. 04 X 10-5 , and T= lOS; all numbers expressed in a. u.The concentration profile evaluated from a 100 trajectory study at t = 500000 a. u.
FIG. 2. The concentration of reactants as a function of dis-tance from the absorber.The solid line is the Langevin result and the dashed line is the diffusion equation result.
FIG.3.The number of absorbed particles as a function of time from single absorber Langevin simulations.
FIG. 4.The ratio of the diffusion equation rate constant to the Fokker-Planck rate constant as a function of the diffusion constant.

FA
I for most phYSical systems, Co is the same order of magnitude as JT.When Co =JT, Eq. (22) becomes6 = R~ In(RB/RA ) • 2rrDT (23)

TABLE I .
Rate constants as a function of diffusion constant.

TABLE II .
Rate constants from single absorber and multiple absorber simulations as a function of diffusion constant.