# Quantum Instanton Evaluations of the Thermal Rate Constants for

Authors: Yi Zhao; Wenji Wang
Journal: Advances in Physical Chemistry
 (2012)

---

## Abstract

Quantum instanton (QI) approximation is recently proposed for the evaluations of the chemical reaction rate constants with use of full dimensional potential energy surfaces. Its strategy is to use the instanton mechanism and to approximate time-dependent quantum dynamics to the imaginary time propagation of the quantities of partition function. It thus incorporates the properties of the instanton idea and the quantum effect of partition function and can be applied to chemical reactions of complex systems. In this paper, we present the QI approach and its applications to several complex systems mainly done by us. The concrete systems include, (1) the reaction ofH+CH4→H2+CH3, (2) the reaction of  H+SiH4→H2+SiH3, (3) H diffusion on Ni(100) surface; and (4) surface-subsurface transport and interior migration for H/Ni. Available experimental and other theoretical data are also presented for the purpose of comparison.


---

## Body

## 1. Introduction

The accurate and efficient evaluation of chemical reaction rate constant is one of prime objectives of theoretical reaction dynamics. Since rigorous quantum mechanical approaches are limited to small molecular (several atoms) reactions, a variety of approximation approaches have been proposed. Benefited from the small recrossing dynamics at not-too-high temperatures, the transition state theories (TSTs), originally proposed by Eyring [1, 2] and Wigner [3],have become a possible and popular way to estimate rate constants. Due to their practical simplicity, they have been broadly applied to numerous reactions. The TST is inherently a classical theory and suitable at sufficiently high temperatures, where the classical description of nuclear motions may be adequate. At low temperatures, especially for the reactions involving the motions of light atoms (i.e., hydrogen), however, quantum effects become quite significant. To make the TST still valid for such low temperature reactions, many approaches have been proposed to quantize it [4–14]. However, there is no absolutely unambiguous way to do it.To develop a more accurate and lessad hoc quantum version of TST, with a specific focus on the tunneling regime, Miller et al. [15–17] have proposed a quantum instanton (QI) approach recently. The QI is based on an earlier semiclassical (SC) TST [18] that became known as the instanton [19, 20]. The similarity between the QI and SC instanton lies in using the steepest descent approximation to evaluate relevant integrals in the quantum rate formula, while the crucial difference is that the Boltzmann operator is evaluated by the quantum mechanics and semiclassical approximation, respectively. The QI theory thus incorporates the tunneling, corner cutting [21–24], and anharmonicity correctly and is expected to overcome the quantitative deficiency of the SC instanton model. In particular, the QI theory considers all tunneling paths and automatically gives each path its naturally weight factor from the quantum Boltzmann operator, instead of choosing a single optimal tunneling path, which is taken into account in the SC instanton and TSTs with SC tunneling corrections. Indeed, it has been numerically demonstrated that the QI predicts accurate quantum rates for one-dimensional and two-dimensional models within 20% error over a wide temperature range, from the deep tunneling to overbarrier regimes.A lot of developments and applications [25–39] have been made since the QI theory is proposed. The original QI [15] involves the second time derivative of the flux-flux correlation at time zero. It has been further improved [29] by taking into account the higher derivatives of the flux-flux function. For the 1D and collinear reactions, the improved model is considerable accuracy, giving the rates to within 5%–10% errors. For a practical purpose, a simple and general way for choosing dividing surfaces used in the QI is suggested [31], namely, using the family of (hyper)planes normal to the minimum energy path at various distances s. A “simplest” QI model [26] has also been suggested with one dividing surface, slightly less accurate than the original QI. To reveal the relationship with conventional TSTs, the classical limit of the QI has been derived [37]. It is found that the classical TST is just a special case of the QI in high temperature limit; moreover, the quantum correction of the prefactor is more important than that of the activation energy in the TST.Since the QI solely involves the Boltzmann operator and its relevant quantities, it can be applied to quite complex molecular systems (from gas phase [17, 37], liquid [28], to surface [38, 39]) via well-established imaginary time path integral techniques. The first implementation of QI with path integral Monte Carlo and adaptive umbrella sampling techniques is applied to the three-dimensional hydrogen exchange reaction D + H2 → HD + H [16]. Soon, the techniques are further extended to the reaction of H+CH4→H2+CH3 [17]. The thermodynamic integration with respect to the mass of the isotope and the inverse temperature is also proposed to compute the kinetic isotope effects [30] and rate constants [34] directly. To improve the convergence of the Monte Carlo simulation, the efficient “virial” estimators [32] have been derived from the logarithmic derivatives of the partition function and the delta-delta correlation functions, and it is found that their statistical errors do not increase with the number of discrete time slices in the path integral. Most recently, the QI has been compared with other conventional approaches [36] for an intramolecular proton transfer on a full-dimensional potential energy surface that incorporates high-levels ab initio calculations along the reaction path. The obtained kinetic isotope effects from the QI are in reasonable agreement with those from the path-integral quantum TST.In this paper, we firstly illustrate the QI formula and its path integral representation. Then, we display several applications, which are mainly done by ourselves. The systems include two gas phase reactionsH+CH4→H2+CH3 [17] and H+SiH4→H2+SiH3 [37], H diffusion on Ni(100) surface [38], and surface-subsurface transport and interior migration for H/Ni [39].


## 2. Method

In this section, we summarize the rate formula for the QI evaluation. The detailed derivation can be found in [15–17]. The QI model proposes the following thermal rate constant:(1)kQI=1Qrπℏ2Cff(0)ΔH(β).
Here, Qr is the reactant partition function per unit volume. Cff(0) is zero time value of the flux-flux correlation function(2)Cff(t)=tr⁡[e-βĤ/2F̂ae-βĤ/2eiĤt/ℏF̂be-iĤt/ℏ],
where β is the inverse temperature 1/(kBT), Ĥ is the Hamiltonian operator of the reaction system, and F̂a and F̂b are the flux operators given by(3)F̂γ=iℏ[Ĥ,h(sγ(r̂))],
with γ=a,b. In (3), h is the step-side function, r represents the Cartesian coordinates of the reaction system, and sa(r) and sb(r) define two separate dividing surfaces via the equations sa(r)=0 and sb(r)=0, both sa(r) and sb(r) being positive (negative) on the product (reactant) side of the dividing surfaces.ΔH(β) in (1) is a specific type of energy variance given by(4)ΔH(β)=ℏ-C̈dd(0)2Cdd(0).
In order to get the correct free particle (high temperature) limit (that would be 25% too large otherwise), an ad hoc term is added to ΔH(β),ΔHmod(β)=ΔH(β)+(π-2)/β, which has very little effect in the low temperature regime. Cdd(0) and C̈dd(0) are zero time value and its second time derivative, respectively, of the “delta-delta” correlation function(5)Cdd(t)=tr⁡[e-βĤ/2Δ(sa(r̂))e-βĤ/2eiĤt/ℏΔ(sb(r̂))e-iĤt/ℏ],
where the generalized delta-function operator is(6)Δ(sγ(r̂))=δ(sγ(r̂))∑i=1N1mi(∇isγ(r̂))2(γ=a,b).
Here, N is the total number of atoms, ∇i=∂/∂ri, ri denotes the Cartesian coordinates of the ith atom, and mi is its atomic mass.The dividing surfaces are determined by the stationary condition(7)∂∂ckCdd(0;{ck})=0,
where {ck} is a collection of parameters that is involved in the location of the dividing surfaces. This condition originates from the SC instanton model, and the resulting dividing surfaces correspond qualitatively to the turning points of the periodic orbit that runs on an upside down PES in imaginary time (see Appendix A in [15, 18]).Since all the relevant quantities in the QI expression (1) involved only the quantum Boltzmann operator, they can be readily evaluated using imaginary time path integral Monte Carlo (PIMC) [40] method.We begin with the simplest quantity,Cdd(0), which can be discretized according to the standard procedure [41–43] to give(8)Cdd(0)=C∫dr(1)∫dr(2)⋯∫dr(P)Δ(sa(r(0)))×Δ(sb(r(P/2)))exp⁡⁡[-βΦ({r(s)})],
where C is a multiplicative constant, P the number of imaginary time slices, and r(s)=(r1(s),r2(s),…,rN(s)) the Cartesian coordinates of the system associated with the sth time slice. Φ({r(s)}) is the discretized action given by(9)Φ({r(s)})=P2ℏ2β2∑s=1P∑i=1Nmi(ri(s)-ri(s-1))2+1P∑s=1PV(r(s)),
where r(0)=r(P) and {r(s)} represents {r(1),r(2),…,r(P)}. Path integral expressions for Cff(0) and C̈dd(0) are somewhat more complicated but can be obtained in a straightforward manner. The appropriate expressions are(10)Cff(0)=C∫dr(1)∫dr(2)⋯∫dr(P)Δ(sa(r(0)))×Δ(sb(r(P/2)))exp⁡⁡[-βΦ({r(s)})]fv({r(s)}),
with(11)fv({r(s)})=(iP2ℏβ)2∑i=1N∇isa(r(0))⋅(ri(1)-ri(P-1))∑i=1Nmi-1(∇isa(r(0)))2×∑i=1N∇isb(r(P/2))⋅(ri(P/2+1)-ri(P/2-1))∑i=1Nmi-1(∇isb(r(P/2)))2,C̈dd(0)=-1ℏ2C∫dr(1)∫dr(2)⋯∫dr(P)Δ(sa(r(0)))×Δ(sb(r(P/2)))exp⁡⁡[-βΦ({r(s)})]×{F({r(s)})2+G({r(s)})},
with(12)F({r(s)})=-Pℏ2β2{∑s=1P/2-∑s=P/2+1P}∑i=1Nmi(ri(s)-ri(s-1))2+2P{∑s=1P/2-1-∑s=P/2+1P-1}V(r(s)),G({r(s)})=2fPβ2-4Pℏ2β3∑s=1P∑i=1Nmi(ri(s)-ri(s-1))2,
with f being the total number of degrees of freedom (i.e., f=3N).In realistic calculations, we rewrite (1) as the product of several ratios(13)kQI=πℏ2Cdd(0)QrCff(0)Cdd(0)1ΔH(β).The terms ofCff(0)/Cdd(0) and ΔH are directly calculated as a constrained average over the same ensemble of paths [16, 17](14)Cff(0)Cdd(0)=〈fv〉,ΔH2=12〈F2+G〉,
with(15)〈⋯〉=∫dr(1)∫dr(2)⋯∫dr(P)Δ(sa(r(0)))Δ(sb(r(P/2)))exp⁡[-βΦ({r(s)})](⋯)∫dr(1)∫dr(2)⋯∫dr(P)Δ(sa(r(0)))Δ(sb(r(P/2)))exp⁡[-βΦ({r(s)})].The evaluation ofCdd(0)/Qr, however, meets a challenge because the Cdd(0) is the quantity associated with the transition state, while Qr with the asymptotic reactant domain, we evaluate it using adaptive umbrella sampling techniques [44].The QI approximation uses the short-time information of the flux-flux correlation function. Predescu and Miller [45] demonstrate that in the classical limit, Wigner’s variational principle and the quantum variational criterion based on the minimization of flux-flux correlation function produce the same optimal surface. Recently, Wang et al. [37] have shown that in the classical limit, the QI formula is the same as the classical transition state theory. These situations motivate us to write the QI formula (13) as the Arrhenius form(16)kQI=AQI⋅e-βΔF.
Here, the free energy ΔF is defined by
(17)ΔF=-1β⋅log⁡⁡(Cdd(0)Qr).
This two-dimensional free energy is related to Qr (the reactant partition function) and Cdd(0;ξa,ξb), and it is corresponding to the quality of probability density at (ξa,ξb).Cdd(0;ξa,ξb) has a similar property to the partition function at the transition state. The prefactor is given by(18)AQI=π2Cff(0)Cdd(0)ℏΔH.
The advantage of (16) allows us to investigate the respective quantum contributions to the rates from the quantum free energy and prefactor by comparing their quantum and classical values, since other factors such as the vibrational-rotational coupling and anharmonicity are automatically involved.


## 3. Results and Discussions

### 3.1. The Reaction ofH+CH4→H2+CH3

TheH+CH4 reaction is a prototype of polyatomic hydrogen abstraction reaction. Quantum dynamical studies of this reaction have become possible only recently, because it involves 12 internal degrees of freedom and thus poses difficulties to quantum dynamics calculations as well as construction of the potential surface. We apply the QI methodology to this reaction using the potential energy surface constructed by Espinosa-García [46]. All calculations are performed in terms of the Cartesian coordinates of all the atoms (i.e., 18 degrees of freedom).In the path integral simulations, the number of imaginary time slicesP is chosen to be 20 and 100 at temperatures T=1000 K and 200 K, respectively, while 3×107 Monte Carlo cycles are run to achieve <10% statistical convergence.


#### 3.1.1. Reaction Coordinate

For this reaction, we define a generalized reaction coordinates(r;ξ), where ξ is an adjustable parameter that shifts the location of the dividing surface (defined by s(r;ξ)=0). s(r;ξ) is defined by a linear interpolation between two constituent reaction coordinates s0(r) and s1(r) through the parameter ξ,(19)s(r;ξ)=ξ⋅s1(r)+(1-ξ)⋅s0(r).s1(r) is a reaction coordinate whose dividing surface is designed to pass through the top of the classical potential barrier, which is defined here as(20)s1(r)=max⁡⁡{sα(r),sβ(r),sγ(r),sδ(r)},
with sx(r)(x=α,β,γ,δ) being the reaction coordinate that describes the abstraction process of one of the methane hydrogens Hx by the incident one H(21)sx(r)=r(C-Hx)-r(Hx-H)-[r†(C-Hx)-r†(Hx-H)],
where r(X-Y) denotes the interatomic distance between atoms X and Y and r†(X-Y) is the value at the transition state geometry. s0(r), on the other hand, describes a dividing surface that is located far in the asymptotic reactant valley, which is given by(22)s0(r)=R∞-|R|.
Here, R is the scattering vector that connects the incident hydrogen and the center of mass of the methane. R∞ is an adjustable parameter which is chosen to be 9 Å in order to guarantee that the interaction potential energy between H and CH4 is negligible.Now, the termCdd(0) (5) becomes a function of two parameters, ξa and ξb, as follows:(23)Cdd(0;ξa,ξb)=tr⁡[e-βĤ/2Δ(s(r̂;ξa))×e-βĤ/2Δ(s(r̂;ξb))],
and thus, one seeks a stationary point of Cdd(0;ξa,ξb) in the two-dimensional (ξa,ξb) space to obtain the corresponding “optimal” values.


#### 3.1.2. Free Energy Surface

The quantityCdd(0;ξa,ξb) varies exponentially as a function of (ξa,ξb), it is convenient to define a quantum “free energy surface” as follows:(24)F(ξa,ξb)=-kBTlog⁡Cdd(0;ξa,ξb),
and locate the saddle point of F(ξa,ξb) by visual inspection.Figure1 shows the free energy surfaces for the H+CH4 reaction. It exhibits a barrier-like profile along the direction ξ̅=(ξa+ξb)/2, while it grows approximately quadratically with the increasing of the absolute value of Δξ=ξa-ξb. From this figure, it is seen that at a higher temperature T=1000K (Figure 1(a)), there appears only a single saddle point at (ξ̅,Δξ)=(1.02,0.0), while at a low temperature T=200K (Figure 1(b)), the saddle point bifurcates into a distinct pair at (ξ̅,Δξ)=(1.1,±0.35), which indicates the existence of nonnegligible tunneling effects in the rate constant [15, 16].Local topography of the quantum free energy surface defined by (24) near the top of the barrier. (a) T=1000K; (b) T=200K. The cross symbols show the location of the saddle points. The values of ξa and ξb at the saddle points are used as input for computing the quantum instanton rate.
(a)(b)


#### 3.1.3. Rate Constants

Having obtained the “optimal” values of the(ξa,ξb) at each temperature, one can now compute the quantum instanton rate by combining various quantities as in (8). The calculated QI rates are tabulated in Table 1 as well as other theoretical and experimental ones.Table 1
Theoretical and experimental rate constants obtained for theH+CH4→H2+CH3 reaction. (ξa,ξb) is the saddle point of the quantum free energy surface. kQI is the quantum instanton rate constant. kCVT/μOMT is the rate constant of the canonical variational theory with microcanonical optimized multidimensional tunneling (CVT/μOMT) [46]. kexpta and kexptb are the experimental Arrhenius fits, k(T)=2.18×10-20T3.0exp(-4045/T)[47] and k(T)=6.78×10-21T3.156exp(-4406/T) [48], respectively. Unit: cm3s-1 for rates.
T (K)(ξa,ξb)kQIkCVT/μOMTkexptakexptb200(0.93, 1.28)3.30 (−22)2.87 (−22)3.36 (−23)300(1.05, 1.05)7.80 (−19)5.03 (−19)8.20 (−19)1.86 (−19)400(1.05, 1.05)6.00 (−17)4.16 (−17)5.66 (−17)1.82 (−17)500(1.04, 1.04)9.00 (−16)6.96 (−16)8.35 (−16)3.33 (−16)600(1.03, 1.03)5.42 (−15)5.00 (−15)5.56 (−15)2.57 (−15)700(1.02, 1.02)2.40 (−14)2.31 (−14)1.19 (−14)800(1.02, 1.02)7.13 (−14)6.72 (−14)7.10 (−14)4.00 (−14)900(1.02, 1.02)1.85 (−13)1.78 (−13)1.07 (−13)1000(1.02, 1.02)3.97 (−13)3.58 (−13)3.82 (−13)2.43 (−13)2000(1.00, 1.00)1.52 (−11)1.70 (−11)2.31 (−11)1.96 (−11)Comparing the quantum instanton rates with others, we find thatkQI is in good agreement with the experimental data more specifically, it is closer to the rates obtained by Baulch et al. [47] than those by Sutherland et al. [48], and kQI agrees with kCVT/μOMT within 10% for the temperature range T = 600–2000 K, but it becomes somewhat larger than the latter as the temperature is decreased (the deviation becomes 30% and 55% for T=500 and 300 K, resp.). It should be noted that the differences between kQI and kCVT/μOMT are much smaller than the uncertainty of the experimental data.Figure2 displays the Arrhenius plots of the rate constants. In Figure 2, we also plot the accurate quantum dynamics results of MCTDH (the multiconfigurational time-dependent Hartree approach) [49]. Compared to the MCTDH ones, our QI rate constants are larger by factors of about 2 to 3 over the temperature range 300–400 K. This difference may partly be due to the recrossing effect which is not considered in QI theory and partly arise from the use of the J-shifting approximation and the neglect of the vibrational angular momenta Hamiltonian in MCTDH method [49].Figure 2
Arrhenius plots of the thermal rate constant for theH+CH4→H2+CH3 reaction. The solid line is the quantum instanton rate. Dotted and dashed lines are the Arrhenius fits of the experimental data from [47, 48], respectively. The dot-dashed line is the result of the multiconfigurational time-dependent Hartree approach (MCTDH) [49]. The open squares are the results of the canonical variational theory with microcanonical optimized multidimensional tunneling (CVT/μOMT) [46].


### 3.2. The Reaction ofH+SiH4→H2+SiH3

TheH+SiH4 reaction is an important step in the radical mechanism of thermal decomposition of monosilane. We calculate the rates and kinetic isotope effects (KIEs) of this reaction with the quantum instanton method in full Cartesian space, on the basis of analytical potential energy surface constructed by Espinosa-García et al. [50].The reaction coordinate ofH+SiH4 has the same form as that of H+CH4 (Section 3.1.1). In our QI calculations, the number of time slices, P, is set to be 20 (1000 K) and 120 (200 K) for the quantum evaluations. In our classical evaluations, the formula is the same as QI, but the number oftime slices in the path integral is set to be 1. The number of Monte Carlo is about (6-10)×106for computing a single ensemble average, and it converges the relevant quantities within 10% statistical errors.


#### 3.2.1. Free Energy and Prefactor

We have rewritten the QI formula (8) in the Arrhenius form (16), which consists of the free energy and prefactor. In this section, we calculate the corresponding quantum and classical quantities so as to investigate the respective quantum contributions to the rates from the quantum free energy and prefactor. Our calculated results are plotted in Figures 3 and 4.Figure 3
The temperature dependence of the free energy. The solid and dotted lines correspond to the quantum and classical calculations, respectively.Figure 4
The temperature dependence of the prefactor. The solid and dotted lines correspond to the quantum and classical calculations, respectively.Now, we look into the quantum effects from both the free energy and prefactor. Figure3 plots the quantum and classical free energies as a function of temperature at the optimized stationary point. One immediately observes that the quantum effect becomes significant at T < 600 K, whereas both quantum and classical results nearly coincide at T > 600 K. As expected, the quantum effect always decreases the classical free energy. The difference of the contributions to the rates is about 40% at 200 K.Figure4 shows the temperature dependence of the quantum and classical prefactors. Compared with Figure 3, one easily finds that this quantum contribution is much larger than that for the free energy correction. Even at 1000 K, the difference is observable. At 200 K, the difference is about several orders of magnitude.In the QI theory, it is not possible to explicitly distinguish the quantum contributions from the partition function and nuclear tunneling, but we can conclude that it is insufficient to estimate the accurate rate by only replacing the activation energy in the TST with its quantum analog, because the quantum prefactor plays more important rule in determining the quantum rate.


#### 3.2.2. Rate Constants

In Table2, we summarize the rate constants obtained from the QI theory and classical VTST (this kind of results are performed by setting the time slices in the path integral to be 1 in the QI calculations) as well as those from the canonical variational TST with the centrifugal-dominant small curvature SC adiabatic ground-state (CVT/CD-SCSAG) approach [50], the conventional TST with simple Wigner tunneling factor [51] and the experiment [52]. The corresponding Arrhenius plots are displayed in Figure 5.Table 2
Rate constants for theH+SiH4 reaction (cm3s-1). The CVT/CD-SCSAG results from [50]. DD is the conventional TST with simple Wigner tunneling factor [51]. Expt. is the experimental Arrhenius fits, k(T)=1.78×10-10exp(-3820/RT) (T = 290–660 K) [52].
T(K)Present QIVTSTCVT/CD-SCSAGDDExpt.2004.82 (−14)8.88 (−16)3.07 (−14)6.1 (−15)3002.74 (−13)1.38 (−14)2.09 (−13)1.7 (−13)2.81 (−13)4009.67 (−13)4.32 (−13)7.96 (−13)9.1 (−13)1.43 (−12)5002.37 (−12)1.23 (−12)2.11 (−12)2.7 (−12)3.78 (−12)6003.99 (−12)2.55 (−12)4.43 (−12)5.9 (−12)7.18 (−12)7007.60 (−12)5.18 (−12)8.00 (−12)8001.18 (−11)8.63 (−12)9001.76 (−11)1.29 (−11)10002.39 (−11)1.99 (−11)2.80 (−11)Figure 5
Arrhenius plot of the thermal rate constants for theH+SiH4 reaction. Solid line, the QI results; dashed line, the classical VTST results; dotted line, CVT/CD-SCSAG results from [50]; dot-dashed, the values from the conventional TST with a simple Wigner tunneling factor [51]; filled squares are the experimental values from [52].Table2 and Figure 5 display comparable results of the CVT/CD-SCSAG to the QI values, with 36% maximal errors at 200 K and a slightly different Arrhenius slope. Both results are in good agreement with experimental data in the tested temperature range. This manifests that the PES used is reasonable accurate. However, the values from the classical VTST are always smaller than the QI results especially in the deep tunneling regions. The conventional TST results have similar tendency to those from the classical VTST. Although the classical VTST can be much improved by evaluating the partition functions quantum mechanically, we do not focus on this improvement. However, it is noted that the anharmonicity, rotational-vibrational coupling are involved in the classical simulation. It is thus expected that these large errors come from the pure quantum effects.


#### 3.2.3. Kinetic Isotope Effects

Kinetic isotope effect (KIE) is the characteristic of chemical reactions which may reveal the quantum effect. We consider the following isotopic reactions:(25)(R1)H+SiH4⟶SiH3+H2(R2)D+SiH4⟶SiH3+HD(R3)H+SiD4⟶SiD3+HD.
The calculated values are tabulated in Table 3.Table 3
The kinetic isotope effects (KIEs) at different temperatures. CVT is the CVT/CD-SCSAG values from [50].
T (K)QI KIE (R1/R2)CVT KIE (R1/R2)QI KIE (R1/R3)CVT KIE (R1/R3)2000.560.6930.1033.603000.730.718.579.814000.770.764.395.075000.850.813.113.466000.900.842.412.727000.920.862.152.328000.941.789000.941.6810000.950.901.641.80The KIEs ofkQI(R1)/kQI(R2) and kQI(R1)/kQI(R3) in the temperature range of 200–1000 K are summarized in Table 3. Besides the QI values, this table also displays comparable results of the CVT/CD-SCSAG [50]. It is easily found that the k(R1)/k(R2) KIEs predicted by the QI are smaller than 1, in agreement with the CVT/CD-SCSAG values. The detailed comparison reveals that the CVT/CD-SCSAG predicts smaller values than those from QI theory. Again, the maximal error occurs at 200 K and is about 23%. For k(R1)/k(R3) KIEs, we find that although both approaches predict the values larger than 1 at 200 to 1000 K, the QI values are smaller than those from CVT/CD-SCSAG. Espinosa-García and coworkers [50] have pointed out that the CVT/CD-SCSAG rates may have been overestimated because of the high vibrational and tunneling contributions. This manifests that the QI approach indeed correctly accounts for the quantum effects.


### 3.3. H Diffusion on Ni(100) Surface

Diffusion plays a fundamental role in surface process. It reveals characteristics about the underlying surface potential and is intimately involved in determining the kinetics of surface catalyzed chemical reactions. The hydrogen atom and its isotopes are ideal candidates to exhibit quantum tunneling behavior due to their small masses. We explore the evaluation of the quantum instanton approximation to the process of H diffusion on Ni(100) surface using the EAM4 potential energy surface constructed by Truong and Truhlar [53].In the path integral calculations, the numbers of time slices,P and Pbath, for the degrees of freedom of H and quantized Ni atoms, respectively, are set to (P,Pbath)= (24–40, 6–8) in the temperature range of 40–300 K. The number of Monte Carlo cycle employed is (1–10) × 106, which converges most of the relevant quantities within 10% statistical error (some of the statistical errors are within 20% below 100 K).


#### 3.3.1. Model

Nickel crystallizes in a face-centered-cubic (fcc) lattice structure. The structural model used consists of 162 Ni atoms over four layers: 40 atoms are in each of the first and third layers, and 41 atoms are in each of the second and fourth layers.Figure6 shows the platform and profile chart of Ni(100) lattice. In the simulations, 89 Ni atoms (orange), in the four sides of the bulk metal and the bottom layer, are fixed. The 48 atoms (blue) closest to the fixed ones are considered to be movable and treated classically. The last 25 ones (red), surrounding the reactant and product sites and lying directly beneath the reactant and product sites, are treated quantum mechanically.Figure 6
The platform and profile chart of H diffusion on Ni(100) lattice. The gray circles represent H atom on different surface sites, and the orange circles represent the fixed Ni atoms. Blue circles are Ni atoms treated classically, while the red ones are treated quantum mechanically.In our calculations, the lattice (see Figure6) is rotated by 45 degrees counterclockwise, and we chose the x coordinate of H atom (along the reaction path) to be reaction coordinate.


#### 3.3.2. Probability Distribution of Paths

The path integral has an advantage that the paths of the particles can display the character of the diffusive process. We extract the probability distributions for the paths of H and Ni atoms (two reaction coordinate beads (x0 and xP/2) of H path are fixed at the two dividing surfaces) for the purpose of guaranteeing the instanton property (the instanton is a periodic orbit between the two dividing surfaces). The dividing surfaces can be obtained from the free energy surface (24), their values at different temperatures are shown in Table 4.Table 4
The positions of two dividing surfaces. (Unit in Å. The position 0.00 corresponds to the bridge site and the positions±1.24Å correspond to two hollow sites.)
T (K)30080605040xa†0.000.00−0.75−1.00−1.08xb†0.000.000.751.001.08Figure7 displays the probability distributions for both H and Ni atoms at several temperatures. It is found that the probability distribution of H is localized at 300 K, and it becomes more and more delocalized as the temperature decreases. Below 80 K, the distribution rapidly becomes obvious in both sides of the transition state, and the path begins to continuously distribute between the two hollow sites at lower temperatures. This phenomenon indicates that H tunneling becomes remarkable at the temperatures lower than 80 K. The probability distributions of Ni atoms have small changes in the whole temperature range, and they seem to be frozen at very low temperature.Figure 7
Distribution of H and Ni quantum paths on the surface with two reaction coordinate beads (x0 and xP/2) of H path fixed at the two dividing surfaces. Labels of <010> and <001> denote the crystal directions of nickel. The probabilities are normalized for both H and Ni.


#### 3.3.3. Diffusion Coefficients

According to the hopping model [54], the diffusion coefficient D is related to the rates through(26)D=kQI⋅b2,
where b is the hop length and is equal to 2.489 Å for the Ni(100) surface. We have calculated the diffusion coefficients at temperatures ranging from 40 to 300 K, and the results are plotted in Figure 8.Figure 8
Arrhenius plots of the diffusion coefficients in the range 40–300 K. Solid line: the QI results for the quantized lattice; dotted line: the QI results for the rigid lattice; dashed line: the CVT/SCSAG results for the rigid lattice [55]. Pluses, circles, squares, and crosses are experimental results of Lee et al. [58], George et al. [56], Mullins et al. [57], and Lin and Gomer [59], respectively.We plot two kinds of QI diffusion coefficients in Figure8, one is the result on a rigid lattice, the other is the result on a quantized lattice. At high temperatures, these two kinds of diffusion coefficients are nearly the same; however, at low temperatures, the ones on a rigid lattice are bigger than that on a quantized lattice, and this is mainly due to the fact that the free energy on a quantized lattice is higher than that on a rigid lattice [38]. Comparing the QI diffusion coefficients with the CVT/SCSAG ones [55] on a rigid lattice, we find that our QI results on a rigid lattice are very similar to those from the CVT/SCSAG. Now, we compare the present theoretical results with experimental ones. It is found that the QI results as well as the CVT/SCSAG results are in good agreement with the experiments of George et al. [56] and Mullins et al. [57] at 200 K and 300 K. However, at low temperatures, the theoretical results are significantly larger than Lee et al. [58] and Lin and Gomer [59]. The experimental transition temperature (100 K–160 K) is also different from the present calculations. The QI predicts it to be 70 K for the quantized lattice, while it is about 66 K for the rigid lattice obtained from both the QI and the CVT/SCSAG approaches. These discrepancies between theories and experiments may be attributed to the interaction potential. It is likely that the accuracy of the potential energy surfaces needs further improvement.


### 3.4. Surface-Subsurface Transport and Interior Migration for H/Ni

A full-dimensional potential energy surface (EAM6), constructed by Wonchoba and Truhlar [60], using the embedded diatomics-in-molecules (EDIMs) [61] potential energy function, is used in this QI calculations.In path integral calculations, the numbers of time slices,P and Pbath for the degrees of freedom of the H and quantized Ni atoms, respectively, are set to (P,Pbath)=(30-40,6-8) in the temperature range of 100–400 K. The number of Monte Carlo is about (2-6)×106 for computing a single ensemble average. It converges most of the values within 10% statistical errors (some of the statistical errors are within 20% at 100 K).


#### 3.4.1. Model

Figures9 and 10 show H diffusion processes in the nickel crystal with a face-centered-cubic (fcc) lattice structure. For a given rate process, we construct a lattice cell (each cell contains more than 200 Ni atoms) in which all atoms are treated to be movable to incorporate the effect of the crystal fluctuation on the rates. To be concrete, the Ni atoms in the four sides of the bulk metal and at the bottom layers are fixed, the Ni atoms surrounding the reactant and product sites and lying along the reaction path are treated quantum mechanically, and the others are treated classically. It should be noted that the lattice for the interior diffusion process has a structure of sphere and the outer layers are fixed.Figure 9
A lattice model with a few Ni atoms for the circle reaction processes of hydrogen.Ahcp,Afcc,Ocsub, and Tesub are the abbreviations for a hcp hollow site, a fcc hollow site, a subsurface octahedral vacancy, and a subsurface tetrahedral vacancy, respectively. The black lines stand for the directions of the reaction paths. The atoms of Ni1,Ni2, and Ni3 colored in orange are specially used to determine the reaction coordinate of the process Ocsub→Tesub.Figure 10
A lattice model with a few Ni atoms for the process of H diffusion in interior of bulk Ni. Ocinte and Teinte are the abbreviations for an interior octahedral vacancy and an interior tetrahedral vacancy, respectively. The black lines stand for the general directions of the reaction paths. The atoms of Ni4,Ni5, and Ni6 colored in orange are specially used to determine the reaction coordinate of the process Ocinte→Teinte.In the QI calculations, we need to define the reaction coordinate operatorsŝ  (6). For the systems considered in this section, the hydrogen coordinates are essentially good choices. We thus adopt the following reaction coordinates for different rate processes. In the H diffusion on Ni(111) surface from a hcp site to a fcc site (Ahcp→Afcc), the x coordinate of the H atom is chosen, whose direction is showed in Figure 9 as the black line connecting the Ahcp site to Afcc site. In the H resurfacing from a subsurface octahedral vacancy to the fcc site (Ocsub→Afcc) and from a subsurface tetrahedral vacancy to the hcp site (Tesub→Ahcp), the z coordinate of the H atom (vertical to Ni(111) surface) is taken. In the H diffusion between the adjacent subsurface octahedral and tetrahedral vacancies (Ocsub→Tesub) and between the adjacent interior octahedral vacancy and interior tetrahedral vacancy (Ocinte→Teinte), the reaction coordinates are along the directions that perpendicular to the planes of Ni1–Ni2–Ni3 and Ni4–Ni5–Ni6, respectively.


#### 3.4.2. Free Energy and Prefactor

In order to investigate the quantized Ni lattice effect on the rates, we have recast the QI formula in (1) into Arrhenius form, which consists of the free energy and the prefactor (17) and (18).In the free energy calculations, firstly, we calculate the free energy profile for each step of the H hopping paths,Ahcp→Afcc→Ocsub→Tesub→Ahcp and Ocinte→Teinte→Ocinte, with the reaction coordinate defined in Section 3.4.1. Then, we connect these free energy profiles one by one, the final free energy profiles for the whole processes are displayed in Figures 11 and 12.Figure 11
Free energy profiles with respect to the processes in Figure9 at 300 K. The green, the red, and the blue lines correspond to the rigid, classical, and quantized lattices, respectively.Figure 12
Free energy profiles with respect to the processes in Figure10 at 300 K. The green, the red, and the blue lines correspond to the rigid, classical, and quantized lattices, respectively.Figures11 and 12 display the calculated free energy profiles for the surface-subsurface and interior processes, respectively, with the quantized lattice, the classical lattice and the rigid one at room temperature (300 K). The corresponding free energy barriers, prefactors, and rates are tabulated in Table 5.Table 5
Free energies, prefactors, and rate constants.
Process  (300 K)ΔF  (kcal/mol)AQI (s-1)kQI (s-1)a2.191.58 (13)4.07 (11)Afcc→Ahcpb2.141.61 (13)4.49 (11)c2.141.27 (13)3.47 (11)a33.066.87 (13)6.22 (−11)Afcc→Ocsubb13.792.85 (13)2.64 (3)c13.752.29 (13)2.23 (3)a20.006.44 (13)1.83 (−1)Ocsub→Afccb5.402.10 (13)2.62 (9)c5.401.72 (13)2.15 (9)a34.924.29 (13)1.67 (−12)Ahcp→Tesubb14.92c14.87Tesub→Ahcpa3.255.13 (13)2.28 (11)a26.195.23 (13)4.31 (−6)Ocsub→Tesubb7.831.86 (13)3.51 (7)c7.861.45 (13)2.74 (7)a7.838.17 (13)1.51 (8)Tesub→Ocsubb1.751.91 (13)1.05 (12)c1.731.49 (13)8.22 (11)a27.825.12 (13)2.86 (−7)Ocinte→Teinteb11.142.50 (13)1.90 (5)c11.171.75 (13)1.32 (5)a7.077.80 (13)6.16 (8)Teinte→Ocinteb1.522.06 (13)1.60 (12)c1.471.53 (13)1.36 (12)a: The results for a rigid lattice.b: The results for a classical lattice.c: The results for a quantized lattice.Figures11 and 12 clearly show that the classical lattices always reduce the free energy barriers compared with the rigid lattices, however, the differences between the free energy barriers with the quantized lattices and the ones with the classical lattices are very small. For different processes; however, the relaxation effect on the free energies is very different. For hydrogen diffusion on Ni(111), the classical lattice only slightly lowers the free energy barrier, while it decreases the barriers by more than one-half in subsurface and interior processes. More careful analysis from Figure 11 reveals that the two preferred Ni(111) surface binding sites, that is, hcp and fcc hollows have symmetric wells, manifesting that the motions of Ni atoms in the layer beneath the surface have little influence on the surface free energies despite the fact that the Ni atoms beneath the hcp and fcc hollows have different arrangements. It is also found that the Ocsub site has a deeper well than the Tesub site and the well of the Tesub site nearly disappears as relaxed Ni atoms are considered. It manifests that the hydrogen at the Tesub site is very unstable and can easily move to the Ocsub site or resurface to the hcp site. In the interior processes, although the hydrogen in the Teinte site is much less stable than in the Ocinte site (see Figure 12), the well at the Teinte is obvious. One thus expects that a two-step reaction process can be used for the reaction from one Ocinte to the other.Next, we consider the prefactor. Table5 shows that the prefactors of the classical lattices are much smaller than those of rigid ones except for the process of H diffusion on Ni(111), and the prefactors of quantized lattices are always smaller than that of classical ones, but their amplitudes have the same order. This may be explained by the fact that the quantum motions of the lattice atoms can induce the dissipative effect on the tunneling degrees of freedom [62, 63], because the prefactor essentially incorporates the dynamical effect. It is well known that pure dissipation in the overdamping regime always hinders the reaction rates for a given reaction barrier. The present results thus are consistent with above analysis.The rate is determined by both the prefactor and the free energy barrier. In Table5, the rate with the rigid lattice is smaller than the one on the classical lattice, which is because the prefactor changes a little and the free energy determines the rate. However, the rate with the quantized lattice is smaller than that of the classical lattice, which is due to the fact that quantized lattice has a smaller prefactor while the free energies are similar. Generally speaking, for the quantized lattices, the rates are lower by 20%–40% when compared to the ones on the classical lattices.Another important feature of the free energy is its temperature dependence. Figure13 displays the free energy profiles with both H- and Ni-treated quantum mechanically at several temperatures. Generally speaking, the free energies have a slight difference at 300 K and 400 K, whereas this difference becomes pronounced for 100 K and 200 K, and the barrier positions move to the directions of shallow well for asymmetric reactions. These properties can be explained by the hydrogen tunneling effect. At lower temperatures, the tunneling plays a more important role. Indeed, the barrier heights decrease with decreasing temperature except for Ocsub→Tesub among 200 K to 400 K. This special case may be due to the special structure of the lattice. The thermal average displacements of Ni1 and Ni3 (in Figure 9) vertical to the Ni(111) surface increase with increasing temperature. Thus, H goes through reaction bottleneck easier at a higher temperature, which makes the barrier decrease. Figure 13 also displays that the free energy barriers of Ocinte→Teinte change little in the temperature range of 200–400 K, while those of the Teinte→Ocinte become smaller and smaller with decreasing temperature. The corresponding free energy barriers are 0.47, 1.09, 1.47, and 1.66 kcal/mol at 100, 200, 300, and 400 K, respectively. It manifests that the diffusing H atom may not equilibrate in the interior tetrahedral vacancy at very low temperatures.The temperature dependence of the free energies. The blue, the black, the red, and the green lines correspond to free energy profiles at 100 K, 200 K, 300 K, and 400 K, respectively.
(a)(b)(c)(d)


#### 3.4.3. Surface-Subsurface Transport

In the resurfacing process,Tesub→Ahcp does not show an obvious barrier (see Figure 11), as the lattice atoms are treated quantum mechanically. This step thus can be thought as a barrierless process. Table 6 tabulates the rate constants for the other resurfacing and subsurface processes in the temperature range of 100 to 600 K. We also list the available CVT/SCT results for Ocsub→Afcc and its reverse reaction [60]. Again, the CVT/SCT rates are close to the QI rates except at 100 K. We think our much bigger rate constant at 100 K is due to the contribution of tunneling.Table 6
Rate constants (Unit:s-1. Powers of 10 are in parentheses) for surface-subsurface processes. aThe QI results with a quantized lattice. bThe CVT/SCT results from [60].
100 K200 K300 K400 K600 KAfcc→Ocsuba3.81 (−14)1.14 (−1)2.23 (3)5.48 (5)2.43 (8)Afcc→Ocsubb1.09 (−15)1.23 (−1)9.72 (3)3.36 (6)1.16 (9)Ocsub→Afcca1.17 (6)4.59 (7)2.15 (9)1.23 (10)1.02 (11)Ocsub→Afccb1.62 (5)7.36 (7)2.93 (9)2.33 (10)2.15 (11)Ocsub→Tesuba2.78 (−4)1.99 (4)2.74 (7)1.18 (9)2.49 (10)Tesub→Ocsuba7.24 (8)4.68 (10)8.22 (11)2.00 (12)4.13 (12)


#### 3.4.4. Interior Migration

In the interior of bulk Ni, the two most stable sites to cage H are symmetric octahedral vacancies (see Figure12). H diffusion between them has been measured experimentally [64–67]. Several theoretical calculations have also been proposed to investigate this diffusion process. Wimmer et al. [68] calculate the diffusion coefficient via two-step reactions Ocinte→Teinte→Ocinte by using a transition state theory together with accurate ab initio energies, while Wonchoba and Truhlar [60] consider the kinetic step as a direct process with a double maximum barrier and calculate the diffusion coefficient using the CVT/SCT. In the present QI calculations, Figure 13 has explicitly shown that the free energies have a well at the tetrahedral site from 200 K to 400 K. It is thus reasonable to assume that the diffusing H atom temporarily equilibrates in the tetrahedral site before jumping forward or backward to a neighboring octahedral site. The free energy well, however, becomes very shallow at 100 K. In this case, the direct reaction from the octahedral site to the other one may be acceptable. Here, we only calculate the diffusion coefficients via the two kinetic steps at 200–400 K.The temperature dependence of the diffusion coefficients are commonly fitted to the Arrhenius equation(27)D(T)=D0⋅exp⁡⁡[-EaRT],
where R is the gas constant and D0 and Ea are the pre-exponential factor and the activation energy, respectively. The QI calculations predict D0=3.93×10-3cm2s-1 and Ea=10.26kcal/mol. In the calculations, the rates of Ocinte→Teinte are used to obtain the diffusion coefficients for Ocinte→Ocinte, because this process is much slower than that of Teinte→Ocinte and it determines the total reaction rates.Table7 tabulates the pre-exponential factors and activation energies coming from available experiments and theories. It is found that both D0 and Ea from the QI calculations are close to Ebisuzaki’s experimental data [65]. Further tracking down the comparisons with experiments is nontrivial, because the accuracy of the diffusion coefficients is much dependent of the potential energy surface. However, we may make a quantitative comparison for the QI and CVT/SCT results. With use of the diffusion coefficients from 295 K to 300 K obtained by the CVT/SCT, Wonchoba and Truhlar [60] predict 11.1 kcal/mol for Ea, and 1.3×10-3cm2s-1 for D0, respectively. These values are observably different from the QI calculations. Ea and D0 are 0.8 kcal/mol larger and 3 times smaller than those from the QI calculations, respectively. The origin of these differences can be explained by that Wonchoba and Truhlar have treated the processes Ocinte→Teinte and Teinte→Ocinte as a single kinetic step rather than as two kinetic steps, which results in a much longer tunneling path than that of the two steps. Compared with the results reported by Wimmer et al. [68], the activation energy is about 0.6 kcal/mol larger than the present one, and the pre-exponential factor is 10 times larger, which is also larger than all available experimental data. These differences may come from both the different potential energy surfaces and rate methods.Table 7
Arrhenius parameters for H diffusion in interior of bulk Ni.aThe QI diffusion coefficients for Ocinte→Ocinte. bFrom [60]. cFrom [68]. dFrom [64]. eFrom [65]. fFrom [66]. gFrom [67].
D0(cm2s-1)Ea(kcal/mol)Present QIa200–400 K3.93 (−3)10.26295–300 K1.3 (−3)11.1Wonchobab (theoretical data)300–627 K7.8 (−4)10.9627–1650 K4.4 (−4)10.2Wimmerc (theoretical data)273–1000 K3.84 (−2)10.92Yamakawad (experimental data)220–350 K1.9 (−3)8.89Ebisuzakie (experimental data)470–690 K5.22 (−3)9.56Eichenauerf (experimental data)660–930 K6.73 (−3)9.47Katzg  (experimental data)670–1270 K7.04 (−3)9.43


## 3.1. The Reaction ofH+CH4→H2+CH3

TheH+CH4 reaction is a prototype of polyatomic hydrogen abstraction reaction. Quantum dynamical studies of this reaction have become possible only recently, because it involves 12 internal degrees of freedom and thus poses difficulties to quantum dynamics calculations as well as construction of the potential surface. We apply the QI methodology to this reaction using the potential energy surface constructed by Espinosa-García [46]. All calculations are performed in terms of the Cartesian coordinates of all the atoms (i.e., 18 degrees of freedom).In the path integral simulations, the number of imaginary time slicesP is chosen to be 20 and 100 at temperatures T=1000 K and 200 K, respectively, while 3×107 Monte Carlo cycles are run to achieve <10% statistical convergence.


### 3.1.1. Reaction Coordinate

For this reaction, we define a generalized reaction coordinates(r;ξ), where ξ is an adjustable parameter that shifts the location of the dividing surface (defined by s(r;ξ)=0). s(r;ξ) is defined by a linear interpolation between two constituent reaction coordinates s0(r) and s1(r) through the parameter ξ,(19)s(r;ξ)=ξ⋅s1(r)+(1-ξ)⋅s0(r).s1(r) is a reaction coordinate whose dividing surface is designed to pass through the top of the classical potential barrier, which is defined here as(20)s1(r)=max⁡⁡{sα(r),sβ(r),sγ(r),sδ(r)},
with sx(r)(x=α,β,γ,δ) being the reaction coordinate that describes the abstraction process of one of the methane hydrogens Hx by the incident one H(21)sx(r)=r(C-Hx)-r(Hx-H)-[r†(C-Hx)-r†(Hx-H)],
where r(X-Y) denotes the interatomic distance between atoms X and Y and r†(X-Y) is the value at the transition state geometry. s0(r), on the other hand, describes a dividing surface that is located far in the asymptotic reactant valley, which is given by(22)s0(r)=R∞-|R|.
Here, R is the scattering vector that connects the incident hydrogen and the center of mass of the methane. R∞ is an adjustable parameter which is chosen to be 9 Å in order to guarantee that the interaction potential energy between H and CH4 is negligible.Now, the termCdd(0) (5) becomes a function of two parameters, ξa and ξb, as follows:(23)Cdd(0;ξa,ξb)=tr⁡[e-βĤ/2Δ(s(r̂;ξa))×e-βĤ/2Δ(s(r̂;ξb))],
and thus, one seeks a stationary point of Cdd(0;ξa,ξb) in the two-dimensional (ξa,ξb) space to obtain the corresponding “optimal” values.


### 3.1.2. Free Energy Surface

The quantityCdd(0;ξa,ξb) varies exponentially as a function of (ξa,ξb), it is convenient to define a quantum “free energy surface” as follows:(24)F(ξa,ξb)=-kBTlog⁡Cdd(0;ξa,ξb),
and locate the saddle point of F(ξa,ξb) by visual inspection.Figure1 shows the free energy surfaces for the H+CH4 reaction. It exhibits a barrier-like profile along the direction ξ̅=(ξa+ξb)/2, while it grows approximately quadratically with the increasing of the absolute value of Δξ=ξa-ξb. From this figure, it is seen that at a higher temperature T=1000K (Figure 1(a)), there appears only a single saddle point at (ξ̅,Δξ)=(1.02,0.0), while at a low temperature T=200K (Figure 1(b)), the saddle point bifurcates into a distinct pair at (ξ̅,Δξ)=(1.1,±0.35), which indicates the existence of nonnegligible tunneling effects in the rate constant [15, 16].Local topography of the quantum free energy surface defined by (24) near the top of the barrier. (a) T=1000K; (b) T=200K. The cross symbols show the location of the saddle points. The values of ξa and ξb at the saddle points are used as input for computing the quantum instanton rate.
(a)(b)


### 3.1.3. Rate Constants

Having obtained the “optimal” values of the(ξa,ξb) at each temperature, one can now compute the quantum instanton rate by combining various quantities as in (8). The calculated QI rates are tabulated in Table 1 as well as other theoretical and experimental ones.Table 1
Theoretical and experimental rate constants obtained for theH+CH4→H2+CH3 reaction. (ξa,ξb) is the saddle point of the quantum free energy surface. kQI is the quantum instanton rate constant. kCVT/μOMT is the rate constant of the canonical variational theory with microcanonical optimized multidimensional tunneling (CVT/μOMT) [46]. kexpta and kexptb are the experimental Arrhenius fits, k(T)=2.18×10-20T3.0exp(-4045/T)[47] and k(T)=6.78×10-21T3.156exp(-4406/T) [48], respectively. Unit: cm3s-1 for rates.
T (K)(ξa,ξb)kQIkCVT/μOMTkexptakexptb200(0.93, 1.28)3.30 (−22)2.87 (−22)3.36 (−23)300(1.05, 1.05)7.80 (−19)5.03 (−19)8.20 (−19)1.86 (−19)400(1.05, 1.05)6.00 (−17)4.16 (−17)5.66 (−17)1.82 (−17)500(1.04, 1.04)9.00 (−16)6.96 (−16)8.35 (−16)3.33 (−16)600(1.03, 1.03)5.42 (−15)5.00 (−15)5.56 (−15)2.57 (−15)700(1.02, 1.02)2.40 (−14)2.31 (−14)1.19 (−14)800(1.02, 1.02)7.13 (−14)6.72 (−14)7.10 (−14)4.00 (−14)900(1.02, 1.02)1.85 (−13)1.78 (−13)1.07 (−13)1000(1.02, 1.02)3.97 (−13)3.58 (−13)3.82 (−13)2.43 (−13)2000(1.00, 1.00)1.52 (−11)1.70 (−11)2.31 (−11)1.96 (−11)Comparing the quantum instanton rates with others, we find thatkQI is in good agreement with the experimental data more specifically, it is closer to the rates obtained by Baulch et al. [47] than those by Sutherland et al. [48], and kQI agrees with kCVT/μOMT within 10% for the temperature range T = 600–2000 K, but it becomes somewhat larger than the latter as the temperature is decreased (the deviation becomes 30% and 55% for T=500 and 300 K, resp.). It should be noted that the differences between kQI and kCVT/μOMT are much smaller than the uncertainty of the experimental data.Figure2 displays the Arrhenius plots of the rate constants. In Figure 2, we also plot the accurate quantum dynamics results of MCTDH (the multiconfigurational time-dependent Hartree approach) [49]. Compared to the MCTDH ones, our QI rate constants are larger by factors of about 2 to 3 over the temperature range 300–400 K. This difference may partly be due to the recrossing effect which is not considered in QI theory and partly arise from the use of the J-shifting approximation and the neglect of the vibrational angular momenta Hamiltonian in MCTDH method [49].Figure 2
Arrhenius plots of the thermal rate constant for theH+CH4→H2+CH3 reaction. The solid line is the quantum instanton rate. Dotted and dashed lines are the Arrhenius fits of the experimental data from [47, 48], respectively. The dot-dashed line is the result of the multiconfigurational time-dependent Hartree approach (MCTDH) [49]. The open squares are the results of the canonical variational theory with microcanonical optimized multidimensional tunneling (CVT/μOMT) [46].


## 3.1.1. Reaction Coordinate

For this reaction, we define a generalized reaction coordinates(r;ξ), where ξ is an adjustable parameter that shifts the location of the dividing surface (defined by s(r;ξ)=0). s(r;ξ) is defined by a linear interpolation between two constituent reaction coordinates s0(r) and s1(r) through the parameter ξ,(19)s(r;ξ)=ξ⋅s1(r)+(1-ξ)⋅s0(r).s1(r) is a reaction coordinate whose dividing surface is designed to pass through the top of the classical potential barrier, which is defined here as(20)s1(r)=max⁡⁡{sα(r),sβ(r),sγ(r),sδ(r)},
with sx(r)(x=α,β,γ,δ) being the reaction coordinate that describes the abstraction process of one of the methane hydrogens Hx by the incident one H(21)sx(r)=r(C-Hx)-r(Hx-H)-[r†(C-Hx)-r†(Hx-H)],
where r(X-Y) denotes the interatomic distance between atoms X and Y and r†(X-Y) is the value at the transition state geometry. s0(r), on the other hand, describes a dividing surface that is located far in the asymptotic reactant valley, which is given by(22)s0(r)=R∞-|R|.
Here, R is the scattering vector that connects the incident hydrogen and the center of mass of the methane. R∞ is an adjustable parameter which is chosen to be 9 Å in order to guarantee that the interaction potential energy between H and CH4 is negligible.Now, the termCdd(0) (5) becomes a function of two parameters, ξa and ξb, as follows:(23)Cdd(0;ξa,ξb)=tr⁡[e-βĤ/2Δ(s(r̂;ξa))×e-βĤ/2Δ(s(r̂;ξb))],
and thus, one seeks a stationary point of Cdd(0;ξa,ξb) in the two-dimensional (ξa,ξb) space to obtain the corresponding “optimal” values.


## 3.1.2. Free Energy Surface

The quantityCdd(0;ξa,ξb) varies exponentially as a function of (ξa,ξb), it is convenient to define a quantum “free energy surface” as follows:(24)F(ξa,ξb)=-kBTlog⁡Cdd(0;ξa,ξb),
and locate the saddle point of F(ξa,ξb) by visual inspection.Figure1 shows the free energy surfaces for the H+CH4 reaction. It exhibits a barrier-like profile along the direction ξ̅=(ξa+ξb)/2, while it grows approximately quadratically with the increasing of the absolute value of Δξ=ξa-ξb. From this figure, it is seen that at a higher temperature T=1000K (Figure 1(a)), there appears only a single saddle point at (ξ̅,Δξ)=(1.02,0.0), while at a low temperature T=200K (Figure 1(b)), the saddle point bifurcates into a distinct pair at (ξ̅,Δξ)=(1.1,±0.35), which indicates the existence of nonnegligible tunneling effects in the rate constant [15, 16].Local topography of the quantum free energy surface defined by (24) near the top of the barrier. (a) T=1000K; (b) T=200K. The cross symbols show the location of the saddle points. The values of ξa and ξb at the saddle points are used as input for computing the quantum instanton rate.
(a)(b)


## 3.1.3. Rate Constants

Having obtained the “optimal” values of the(ξa,ξb) at each temperature, one can now compute the quantum instanton rate by combining various quantities as in (8). The calculated QI rates are tabulated in Table 1 as well as other theoretical and experimental ones.Table 1
Theoretical and experimental rate constants obtained for theH+CH4→H2+CH3 reaction. (ξa,ξb) is the saddle point of the quantum free energy surface. kQI is the quantum instanton rate constant. kCVT/μOMT is the rate constant of the canonical variational theory with microcanonical optimized multidimensional tunneling (CVT/μOMT) [46]. kexpta and kexptb are the experimental Arrhenius fits, k(T)=2.18×10-20T3.0exp(-4045/T)[47] and k(T)=6.78×10-21T3.156exp(-4406/T) [48], respectively. Unit: cm3s-1 for rates.
T (K)(ξa,ξb)kQIkCVT/μOMTkexptakexptb200(0.93, 1.28)3.30 (−22)2.87 (−22)3.36 (−23)300(1.05, 1.05)7.80 (−19)5.03 (−19)8.20 (−19)1.86 (−19)400(1.05, 1.05)6.00 (−17)4.16 (−17)5.66 (−17)1.82 (−17)500(1.04, 1.04)9.00 (−16)6.96 (−16)8.35 (−16)3.33 (−16)600(1.03, 1.03)5.42 (−15)5.00 (−15)5.56 (−15)2.57 (−15)700(1.02, 1.02)2.40 (−14)2.31 (−14)1.19 (−14)800(1.02, 1.02)7.13 (−14)6.72 (−14)7.10 (−14)4.00 (−14)900(1.02, 1.02)1.85 (−13)1.78 (−13)1.07 (−13)1000(1.02, 1.02)3.97 (−13)3.58 (−13)3.82 (−13)2.43 (−13)2000(1.00, 1.00)1.52 (−11)1.70 (−11)2.31 (−11)1.96 (−11)Comparing the quantum instanton rates with others, we find thatkQI is in good agreement with the experimental data more specifically, it is closer to the rates obtained by Baulch et al. [47] than those by Sutherland et al. [48], and kQI agrees with kCVT/μOMT within 10% for the temperature range T = 600–2000 K, but it becomes somewhat larger than the latter as the temperature is decreased (the deviation becomes 30% and 55% for T=500 and 300 K, resp.). It should be noted that the differences between kQI and kCVT/μOMT are much smaller than the uncertainty of the experimental data.Figure2 displays the Arrhenius plots of the rate constants. In Figure 2, we also plot the accurate quantum dynamics results of MCTDH (the multiconfigurational time-dependent Hartree approach) [49]. Compared to the MCTDH ones, our QI rate constants are larger by factors of about 2 to 3 over the temperature range 300–400 K. This difference may partly be due to the recrossing effect which is not considered in QI theory and partly arise from the use of the J-shifting approximation and the neglect of the vibrational angular momenta Hamiltonian in MCTDH method [49].Figure 2
Arrhenius plots of the thermal rate constant for theH+CH4→H2+CH3 reaction. The solid line is the quantum instanton rate. Dotted and dashed lines are the Arrhenius fits of the experimental data from [47, 48], respectively. The dot-dashed line is the result of the multiconfigurational time-dependent Hartree approach (MCTDH) [49]. The open squares are the results of the canonical variational theory with microcanonical optimized multidimensional tunneling (CVT/μOMT) [46].


## 3.2. The Reaction ofH+SiH4→H2+SiH3

TheH+SiH4 reaction is an important step in the radical mechanism of thermal decomposition of monosilane. We calculate the rates and kinetic isotope effects (KIEs) of this reaction with the quantum instanton method in full Cartesian space, on the basis of analytical potential energy surface constructed by Espinosa-García et al. [50].The reaction coordinate ofH+SiH4 has the same form as that of H+CH4 (Section 3.1.1). In our QI calculations, the number of time slices, P, is set to be 20 (1000 K) and 120 (200 K) for the quantum evaluations. In our classical evaluations, the formula is the same as QI, but the number oftime slices in the path integral is set to be 1. The number of Monte Carlo is about (6-10)×106for computing a single ensemble average, and it converges the relevant quantities within 10% statistical errors.


### 3.2.1. Free Energy and Prefactor

We have rewritten the QI formula (8) in the Arrhenius form (16), which consists of the free energy and prefactor. In this section, we calculate the corresponding quantum and classical quantities so as to investigate the respective quantum contributions to the rates from the quantum free energy and prefactor. Our calculated results are plotted in Figures 3 and 4.Figure 3
The temperature dependence of the free energy. The solid and dotted lines correspond to the quantum and classical calculations, respectively.Figure 4
The temperature dependence of the prefactor. The solid and dotted lines correspond to the quantum and classical calculations, respectively.Now, we look into the quantum effects from both the free energy and prefactor. Figure3 plots the quantum and classical free energies as a function of temperature at the optimized stationary point. One immediately observes that the quantum effect becomes significant at T < 600 K, whereas both quantum and classical results nearly coincide at T > 600 K. As expected, the quantum effect always decreases the classical free energy. The difference of the contributions to the rates is about 40% at 200 K.Figure4 shows the temperature dependence of the quantum and classical prefactors. Compared with Figure 3, one easily finds that this quantum contribution is much larger than that for the free energy correction. Even at 1000 K, the difference is observable. At 200 K, the difference is about several orders of magnitude.In the QI theory, it is not possible to explicitly distinguish the quantum contributions from the partition function and nuclear tunneling, but we can conclude that it is insufficient to estimate the accurate rate by only replacing the activation energy in the TST with its quantum analog, because the quantum prefactor plays more important rule in determining the quantum rate.


### 3.2.2. Rate Constants

In Table2, we summarize the rate constants obtained from the QI theory and classical VTST (this kind of results are performed by setting the time slices in the path integral to be 1 in the QI calculations) as well as those from the canonical variational TST with the centrifugal-dominant small curvature SC adiabatic ground-state (CVT/CD-SCSAG) approach [50], the conventional TST with simple Wigner tunneling factor [51] and the experiment [52]. The corresponding Arrhenius plots are displayed in Figure 5.Table 2
Rate constants for theH+SiH4 reaction (cm3s-1). The CVT/CD-SCSAG results from [50]. DD is the conventional TST with simple Wigner tunneling factor [51]. Expt. is the experimental Arrhenius fits, k(T)=1.78×10-10exp(-3820/RT) (T = 290–660 K) [52].
T(K)Present QIVTSTCVT/CD-SCSAGDDExpt.2004.82 (−14)8.88 (−16)3.07 (−14)6.1 (−15)3002.74 (−13)1.38 (−14)2.09 (−13)1.7 (−13)2.81 (−13)4009.67 (−13)4.32 (−13)7.96 (−13)9.1 (−13)1.43 (−12)5002.37 (−12)1.23 (−12)2.11 (−12)2.7 (−12)3.78 (−12)6003.99 (−12)2.55 (−12)4.43 (−12)5.9 (−12)7.18 (−12)7007.60 (−12)5.18 (−12)8.00 (−12)8001.18 (−11)8.63 (−12)9001.76 (−11)1.29 (−11)10002.39 (−11)1.99 (−11)2.80 (−11)Figure 5
Arrhenius plot of the thermal rate constants for theH+SiH4 reaction. Solid line, the QI results; dashed line, the classical VTST results; dotted line, CVT/CD-SCSAG results from [50]; dot-dashed, the values from the conventional TST with a simple Wigner tunneling factor [51]; filled squares are the experimental values from [52].Table2 and Figure 5 display comparable results of the CVT/CD-SCSAG to the QI values, with 36% maximal errors at 200 K and a slightly different Arrhenius slope. Both results are in good agreement with experimental data in the tested temperature range. This manifests that the PES used is reasonable accurate. However, the values from the classical VTST are always smaller than the QI results especially in the deep tunneling regions. The conventional TST results have similar tendency to those from the classical VTST. Although the classical VTST can be much improved by evaluating the partition functions quantum mechanically, we do not focus on this improvement. However, it is noted that the anharmonicity, rotational-vibrational coupling are involved in the classical simulation. It is thus expected that these large errors come from the pure quantum effects.


### 3.2.3. Kinetic Isotope Effects

Kinetic isotope effect (KIE) is the characteristic of chemical reactions which may reveal the quantum effect. We consider the following isotopic reactions:(25)(R1)H+SiH4⟶SiH3+H2(R2)D+SiH4⟶SiH3+HD(R3)H+SiD4⟶SiD3+HD.
The calculated values are tabulated in Table 3.Table 3
The kinetic isotope effects (KIEs) at different temperatures. CVT is the CVT/CD-SCSAG values from [50].
T (K)QI KIE (R1/R2)CVT KIE (R1/R2)QI KIE (R1/R3)CVT KIE (R1/R3)2000.560.6930.1033.603000.730.718.579.814000.770.764.395.075000.850.813.113.466000.900.842.412.727000.920.862.152.328000.941.789000.941.6810000.950.901.641.80The KIEs ofkQI(R1)/kQI(R2) and kQI(R1)/kQI(R3) in the temperature range of 200–1000 K are summarized in Table 3. Besides the QI values, this table also displays comparable results of the CVT/CD-SCSAG [50]. It is easily found that the k(R1)/k(R2) KIEs predicted by the QI are smaller than 1, in agreement with the CVT/CD-SCSAG values. The detailed comparison reveals that the CVT/CD-SCSAG predicts smaller values than those from QI theory. Again, the maximal error occurs at 200 K and is about 23%. For k(R1)/k(R3) KIEs, we find that although both approaches predict the values larger than 1 at 200 to 1000 K, the QI values are smaller than those from CVT/CD-SCSAG. Espinosa-García and coworkers [50] have pointed out that the CVT/CD-SCSAG rates may have been overestimated because of the high vibrational and tunneling contributions. This manifests that the QI approach indeed correctly accounts for the quantum effects.


## 3.2.1. Free Energy and Prefactor

We have rewritten the QI formula (8) in the Arrhenius form (16), which consists of the free energy and prefactor. In this section, we calculate the corresponding quantum and classical quantities so as to investigate the respective quantum contributions to the rates from the quantum free energy and prefactor. Our calculated results are plotted in Figures 3 and 4.Figure 3
The temperature dependence of the free energy. The solid and dotted lines correspond to the quantum and classical calculations, respectively.Figure 4
The temperature dependence of the prefactor. The solid and dotted lines correspond to the quantum and classical calculations, respectively.Now, we look into the quantum effects from both the free energy and prefactor. Figure3 plots the quantum and classical free energies as a function of temperature at the optimized stationary point. One immediately observes that the quantum effect becomes significant at T < 600 K, whereas both quantum and classical results nearly coincide at T > 600 K. As expected, the quantum effect always decreases the classical free energy. The difference of the contributions to the rates is about 40% at 200 K.Figure4 shows the temperature dependence of the quantum and classical prefactors. Compared with Figure 3, one easily finds that this quantum contribution is much larger than that for the free energy correction. Even at 1000 K, the difference is observable. At 200 K, the difference is about several orders of magnitude.In the QI theory, it is not possible to explicitly distinguish the quantum contributions from the partition function and nuclear tunneling, but we can conclude that it is insufficient to estimate the accurate rate by only replacing the activation energy in the TST with its quantum analog, because the quantum prefactor plays more important rule in determining the quantum rate.


## 3.2.2. Rate Constants

In Table2, we summarize the rate constants obtained from the QI theory and classical VTST (this kind of results are performed by setting the time slices in the path integral to be 1 in the QI calculations) as well as those from the canonical variational TST with the centrifugal-dominant small curvature SC adiabatic ground-state (CVT/CD-SCSAG) approach [50], the conventional TST with simple Wigner tunneling factor [51] and the experiment [52]. The corresponding Arrhenius plots are displayed in Figure 5.Table 2
Rate constants for theH+SiH4 reaction (cm3s-1). The CVT/CD-SCSAG results from [50]. DD is the conventional TST with simple Wigner tunneling factor [51]. Expt. is the experimental Arrhenius fits, k(T)=1.78×10-10exp(-3820/RT) (T = 290–660 K) [52].
T(K)Present QIVTSTCVT/CD-SCSAGDDExpt.2004.82 (−14)8.88 (−16)3.07 (−14)6.1 (−15)3002.74 (−13)1.38 (−14)2.09 (−13)1.7 (−13)2.81 (−13)4009.67 (−13)4.32 (−13)7.96 (−13)9.1 (−13)1.43 (−12)5002.37 (−12)1.23 (−12)2.11 (−12)2.7 (−12)3.78 (−12)6003.99 (−12)2.55 (−12)4.43 (−12)5.9 (−12)7.18 (−12)7007.60 (−12)5.18 (−12)8.00 (−12)8001.18 (−11)8.63 (−12)9001.76 (−11)1.29 (−11)10002.39 (−11)1.99 (−11)2.80 (−11)Figure 5
Arrhenius plot of the thermal rate constants for theH+SiH4 reaction. Solid line, the QI results; dashed line, the classical VTST results; dotted line, CVT/CD-SCSAG results from [50]; dot-dashed, the values from the conventional TST with a simple Wigner tunneling factor [51]; filled squares are the experimental values from [52].Table2 and Figure 5 display comparable results of the CVT/CD-SCSAG to the QI values, with 36% maximal errors at 200 K and a slightly different Arrhenius slope. Both results are in good agreement with experimental data in the tested temperature range. This manifests that the PES used is reasonable accurate. However, the values from the classical VTST are always smaller than the QI results especially in the deep tunneling regions. The conventional TST results have similar tendency to those from the classical VTST. Although the classical VTST can be much improved by evaluating the partition functions quantum mechanically, we do not focus on this improvement. However, it is noted that the anharmonicity, rotational-vibrational coupling are involved in the classical simulation. It is thus expected that these large errors come from the pure quantum effects.


## 3.2.3. Kinetic Isotope Effects

Kinetic isotope effect (KIE) is the characteristic of chemical reactions which may reveal the quantum effect. We consider the following isotopic reactions:(25)(R1)H+SiH4⟶SiH3+H2(R2)D+SiH4⟶SiH3+HD(R3)H+SiD4⟶SiD3+HD.
The calculated values are tabulated in Table 3.Table 3
The kinetic isotope effects (KIEs) at different temperatures. CVT is the CVT/CD-SCSAG values from [50].
T (K)QI KIE (R1/R2)CVT KIE (R1/R2)QI KIE (R1/R3)CVT KIE (R1/R3)2000.560.6930.1033.603000.730.718.579.814000.770.764.395.075000.850.813.113.466000.900.842.412.727000.920.862.152.328000.941.789000.941.6810000.950.901.641.80The KIEs ofkQI(R1)/kQI(R2) and kQI(R1)/kQI(R3) in the temperature range of 200–1000 K are summarized in Table 3. Besides the QI values, this table also displays comparable results of the CVT/CD-SCSAG [50]. It is easily found that the k(R1)/k(R2) KIEs predicted by the QI are smaller than 1, in agreement with the CVT/CD-SCSAG values. The detailed comparison reveals that the CVT/CD-SCSAG predicts smaller values than those from QI theory. Again, the maximal error occurs at 200 K and is about 23%. For k(R1)/k(R3) KIEs, we find that although both approaches predict the values larger than 1 at 200 to 1000 K, the QI values are smaller than those from CVT/CD-SCSAG. Espinosa-García and coworkers [50] have pointed out that the CVT/CD-SCSAG rates may have been overestimated because of the high vibrational and tunneling contributions. This manifests that the QI approach indeed correctly accounts for the quantum effects.


## 3.3. H Diffusion on Ni(100) Surface

Diffusion plays a fundamental role in surface process. It reveals characteristics about the underlying surface potential and is intimately involved in determining the kinetics of surface catalyzed chemical reactions. The hydrogen atom and its isotopes are ideal candidates to exhibit quantum tunneling behavior due to their small masses. We explore the evaluation of the quantum instanton approximation to the process of H diffusion on Ni(100) surface using the EAM4 potential energy surface constructed by Truong and Truhlar [53].In the path integral calculations, the numbers of time slices,P and Pbath, for the degrees of freedom of H and quantized Ni atoms, respectively, are set to (P,Pbath)= (24–40, 6–8) in the temperature range of 40–300 K. The number of Monte Carlo cycle employed is (1–10) × 106, which converges most of the relevant quantities within 10% statistical error (some of the statistical errors are within 20% below 100 K).


### 3.3.1. Model

Nickel crystallizes in a face-centered-cubic (fcc) lattice structure. The structural model used consists of 162 Ni atoms over four layers: 40 atoms are in each of the first and third layers, and 41 atoms are in each of the second and fourth layers.Figure6 shows the platform and profile chart of Ni(100) lattice. In the simulations, 89 Ni atoms (orange), in the four sides of the bulk metal and the bottom layer, are fixed. The 48 atoms (blue) closest to the fixed ones are considered to be movable and treated classically. The last 25 ones (red), surrounding the reactant and product sites and lying directly beneath the reactant and product sites, are treated quantum mechanically.Figure 6
The platform and profile chart of H diffusion on Ni(100) lattice. The gray circles represent H atom on different surface sites, and the orange circles represent the fixed Ni atoms. Blue circles are Ni atoms treated classically, while the red ones are treated quantum mechanically.In our calculations, the lattice (see Figure6) is rotated by 45 degrees counterclockwise, and we chose the x coordinate of H atom (along the reaction path) to be reaction coordinate.


### 3.3.2. Probability Distribution of Paths

The path integral has an advantage that the paths of the particles can display the character of the diffusive process. We extract the probability distributions for the paths of H and Ni atoms (two reaction coordinate beads (x0 and xP/2) of H path are fixed at the two dividing surfaces) for the purpose of guaranteeing the instanton property (the instanton is a periodic orbit between the two dividing surfaces). The dividing surfaces can be obtained from the free energy surface (24), their values at different temperatures are shown in Table 4.Table 4
The positions of two dividing surfaces. (Unit in Å. The position 0.00 corresponds to the bridge site and the positions±1.24Å correspond to two hollow sites.)
T (K)30080605040xa†0.000.00−0.75−1.00−1.08xb†0.000.000.751.001.08Figure7 displays the probability distributions for both H and Ni atoms at several temperatures. It is found that the probability distribution of H is localized at 300 K, and it becomes more and more delocalized as the temperature decreases. Below 80 K, the distribution rapidly becomes obvious in both sides of the transition state, and the path begins to continuously distribute between the two hollow sites at lower temperatures. This phenomenon indicates that H tunneling becomes remarkable at the temperatures lower than 80 K. The probability distributions of Ni atoms have small changes in the whole temperature range, and they seem to be frozen at very low temperature.Figure 7
Distribution of H and Ni quantum paths on the surface with two reaction coordinate beads (x0 and xP/2) of H path fixed at the two dividing surfaces. Labels of <010> and <001> denote the crystal directions of nickel. The probabilities are normalized for both H and Ni.


### 3.3.3. Diffusion Coefficients

According to the hopping model [54], the diffusion coefficient D is related to the rates through(26)D=kQI⋅b2,
where b is the hop length and is equal to 2.489 Å for the Ni(100) surface. We have calculated the diffusion coefficients at temperatures ranging from 40 to 300 K, and the results are plotted in Figure 8.Figure 8
Arrhenius plots of the diffusion coefficients in the range 40–300 K. Solid line: the QI results for the quantized lattice; dotted line: the QI results for the rigid lattice; dashed line: the CVT/SCSAG results for the rigid lattice [55]. Pluses, circles, squares, and crosses are experimental results of Lee et al. [58], George et al. [56], Mullins et al. [57], and Lin and Gomer [59], respectively.We plot two kinds of QI diffusion coefficients in Figure8, one is the result on a rigid lattice, the other is the result on a quantized lattice. At high temperatures, these two kinds of diffusion coefficients are nearly the same; however, at low temperatures, the ones on a rigid lattice are bigger than that on a quantized lattice, and this is mainly due to the fact that the free energy on a quantized lattice is higher than that on a rigid lattice [38]. Comparing the QI diffusion coefficients with the CVT/SCSAG ones [55] on a rigid lattice, we find that our QI results on a rigid lattice are very similar to those from the CVT/SCSAG. Now, we compare the present theoretical results with experimental ones. It is found that the QI results as well as the CVT/SCSAG results are in good agreement with the experiments of George et al. [56] and Mullins et al. [57] at 200 K and 300 K. However, at low temperatures, the theoretical results are significantly larger than Lee et al. [58] and Lin and Gomer [59]. The experimental transition temperature (100 K–160 K) is also different from the present calculations. The QI predicts it to be 70 K for the quantized lattice, while it is about 66 K for the rigid lattice obtained from both the QI and the CVT/SCSAG approaches. These discrepancies between theories and experiments may be attributed to the interaction potential. It is likely that the accuracy of the potential energy surfaces needs further improvement.


## 3.3.1. Model

Nickel crystallizes in a face-centered-cubic (fcc) lattice structure. The structural model used consists of 162 Ni atoms over four layers: 40 atoms are in each of the first and third layers, and 41 atoms are in each of the second and fourth layers.Figure6 shows the platform and profile chart of Ni(100) lattice. In the simulations, 89 Ni atoms (orange), in the four sides of the bulk metal and the bottom layer, are fixed. The 48 atoms (blue) closest to the fixed ones are considered to be movable and treated classically. The last 25 ones (red), surrounding the reactant and product sites and lying directly beneath the reactant and product sites, are treated quantum mechanically.Figure 6
The platform and profile chart of H diffusion on Ni(100) lattice. The gray circles represent H atom on different surface sites, and the orange circles represent the fixed Ni atoms. Blue circles are Ni atoms treated classically, while the red ones are treated quantum mechanically.In our calculations, the lattice (see Figure6) is rotated by 45 degrees counterclockwise, and we chose the x coordinate of H atom (along the reaction path) to be reaction coordinate.


## 3.3.2. Probability Distribution of Paths

The path integral has an advantage that the paths of the particles can display the character of the diffusive process. We extract the probability distributions for the paths of H and Ni atoms (two reaction coordinate beads (x0 and xP/2) of H path are fixed at the two dividing surfaces) for the purpose of guaranteeing the instanton property (the instanton is a periodic orbit between the two dividing surfaces). The dividing surfaces can be obtained from the free energy surface (24), their values at different temperatures are shown in Table 4.Table 4
The positions of two dividing surfaces. (Unit in Å. The position 0.00 corresponds to the bridge site and the positions±1.24Å correspond to two hollow sites.)
T (K)30080605040xa†0.000.00−0.75−1.00−1.08xb†0.000.000.751.001.08Figure7 displays the probability distributions for both H and Ni atoms at several temperatures. It is found that the probability distribution of H is localized at 300 K, and it becomes more and more delocalized as the temperature decreases. Below 80 K, the distribution rapidly becomes obvious in both sides of the transition state, and the path begins to continuously distribute between the two hollow sites at lower temperatures. This phenomenon indicates that H tunneling becomes remarkable at the temperatures lower than 80 K. The probability distributions of Ni atoms have small changes in the whole temperature range, and they seem to be frozen at very low temperature.Figure 7
Distribution of H and Ni quantum paths on the surface with two reaction coordinate beads (x0 and xP/2) of H path fixed at the two dividing surfaces. Labels of <010> and <001> denote the crystal directions of nickel. The probabilities are normalized for both H and Ni.


## 3.3.3. Diffusion Coefficients

According to the hopping model [54], the diffusion coefficient D is related to the rates through(26)D=kQI⋅b2,
where b is the hop length and is equal to 2.489 Å for the Ni(100) surface. We have calculated the diffusion coefficients at temperatures ranging from 40 to 300 K, and the results are plotted in Figure 8.Figure 8
Arrhenius plots of the diffusion coefficients in the range 40–300 K. Solid line: the QI results for the quantized lattice; dotted line: the QI results for the rigid lattice; dashed line: the CVT/SCSAG results for the rigid lattice [55]. Pluses, circles, squares, and crosses are experimental results of Lee et al. [58], George et al. [56], Mullins et al. [57], and Lin and Gomer [59], respectively.We plot two kinds of QI diffusion coefficients in Figure8, one is the result on a rigid lattice, the other is the result on a quantized lattice. At high temperatures, these two kinds of diffusion coefficients are nearly the same; however, at low temperatures, the ones on a rigid lattice are bigger than that on a quantized lattice, and this is mainly due to the fact that the free energy on a quantized lattice is higher than that on a rigid lattice [38]. Comparing the QI diffusion coefficients with the CVT/SCSAG ones [55] on a rigid lattice, we find that our QI results on a rigid lattice are very similar to those from the CVT/SCSAG. Now, we compare the present theoretical results with experimental ones. It is found that the QI results as well as the CVT/SCSAG results are in good agreement with the experiments of George et al. [56] and Mullins et al. [57] at 200 K and 300 K. However, at low temperatures, the theoretical results are significantly larger than Lee et al. [58] and Lin and Gomer [59]. The experimental transition temperature (100 K–160 K) is also different from the present calculations. The QI predicts it to be 70 K for the quantized lattice, while it is about 66 K for the rigid lattice obtained from both the QI and the CVT/SCSAG approaches. These discrepancies between theories and experiments may be attributed to the interaction potential. It is likely that the accuracy of the potential energy surfaces needs further improvement.


## 3.4. Surface-Subsurface Transport and Interior Migration for H/Ni

A full-dimensional potential energy surface (EAM6), constructed by Wonchoba and Truhlar [60], using the embedded diatomics-in-molecules (EDIMs) [61] potential energy function, is used in this QI calculations.In path integral calculations, the numbers of time slices,P and Pbath for the degrees of freedom of the H and quantized Ni atoms, respectively, are set to (P,Pbath)=(30-40,6-8) in the temperature range of 100–400 K. The number of Monte Carlo is about (2-6)×106 for computing a single ensemble average. It converges most of the values within 10% statistical errors (some of the statistical errors are within 20% at 100 K).


### 3.4.1. Model

Figures9 and 10 show H diffusion processes in the nickel crystal with a face-centered-cubic (fcc) lattice structure. For a given rate process, we construct a lattice cell (each cell contains more than 200 Ni atoms) in which all atoms are treated to be movable to incorporate the effect of the crystal fluctuation on the rates. To be concrete, the Ni atoms in the four sides of the bulk metal and at the bottom layers are fixed, the Ni atoms surrounding the reactant and product sites and lying along the reaction path are treated quantum mechanically, and the others are treated classically. It should be noted that the lattice for the interior diffusion process has a structure of sphere and the outer layers are fixed.Figure 9
A lattice model with a few Ni atoms for the circle reaction processes of hydrogen.Ahcp,Afcc,Ocsub, and Tesub are the abbreviations for a hcp hollow site, a fcc hollow site, a subsurface octahedral vacancy, and a subsurface tetrahedral vacancy, respectively. The black lines stand for the directions of the reaction paths. The atoms of Ni1,Ni2, and Ni3 colored in orange are specially used to determine the reaction coordinate of the process Ocsub→Tesub.Figure 10
A lattice model with a few Ni atoms for the process of H diffusion in interior of bulk Ni. Ocinte and Teinte are the abbreviations for an interior octahedral vacancy and an interior tetrahedral vacancy, respectively. The black lines stand for the general directions of the reaction paths. The atoms of Ni4,Ni5, and Ni6 colored in orange are specially used to determine the reaction coordinate of the process Ocinte→Teinte.In the QI calculations, we need to define the reaction coordinate operatorsŝ  (6). For the systems considered in this section, the hydrogen coordinates are essentially good choices. We thus adopt the following reaction coordinates for different rate processes. In the H diffusion on Ni(111) surface from a hcp site to a fcc site (Ahcp→Afcc), the x coordinate of the H atom is chosen, whose direction is showed in Figure 9 as the black line connecting the Ahcp site to Afcc site. In the H resurfacing from a subsurface octahedral vacancy to the fcc site (Ocsub→Afcc) and from a subsurface tetrahedral vacancy to the hcp site (Tesub→Ahcp), the z coordinate of the H atom (vertical to Ni(111) surface) is taken. In the H diffusion between the adjacent subsurface octahedral and tetrahedral vacancies (Ocsub→Tesub) and between the adjacent interior octahedral vacancy and interior tetrahedral vacancy (Ocinte→Teinte), the reaction coordinates are along the directions that perpendicular to the planes of Ni1–Ni2–Ni3 and Ni4–Ni5–Ni6, respectively.


### 3.4.2. Free Energy and Prefactor

In order to investigate the quantized Ni lattice effect on the rates, we have recast the QI formula in (1) into Arrhenius form, which consists of the free energy and the prefactor (17) and (18).In the free energy calculations, firstly, we calculate the free energy profile for each step of the H hopping paths,Ahcp→Afcc→Ocsub→Tesub→Ahcp and Ocinte→Teinte→Ocinte, with the reaction coordinate defined in Section 3.4.1. Then, we connect these free energy profiles one by one, the final free energy profiles for the whole processes are displayed in Figures 11 and 12.Figure 11
Free energy profiles with respect to the processes in Figure9 at 300 K. The green, the red, and the blue lines correspond to the rigid, classical, and quantized lattices, respectively.Figure 12
Free energy profiles with respect to the processes in Figure10 at 300 K. The green, the red, and the blue lines correspond to the rigid, classical, and quantized lattices, respectively.Figures11 and 12 display the calculated free energy profiles for the surface-subsurface and interior processes, respectively, with the quantized lattice, the classical lattice and the rigid one at room temperature (300 K). The corresponding free energy barriers, prefactors, and rates are tabulated in Table 5.Table 5
Free energies, prefactors, and rate constants.
Process  (300 K)ΔF  (kcal/mol)AQI (s-1)kQI (s-1)a2.191.58 (13)4.07 (11)Afcc→Ahcpb2.141.61 (13)4.49 (11)c2.141.27 (13)3.47 (11)a33.066.87 (13)6.22 (−11)Afcc→Ocsubb13.792.85 (13)2.64 (3)c13.752.29 (13)2.23 (3)a20.006.44 (13)1.83 (−1)Ocsub→Afccb5.402.10 (13)2.62 (9)c5.401.72 (13)2.15 (9)a34.924.29 (13)1.67 (−12)Ahcp→Tesubb14.92c14.87Tesub→Ahcpa3.255.13 (13)2.28 (11)a26.195.23 (13)4.31 (−6)Ocsub→Tesubb7.831.86 (13)3.51 (7)c7.861.45 (13)2.74 (7)a7.838.17 (13)1.51 (8)Tesub→Ocsubb1.751.91 (13)1.05 (12)c1.731.49 (13)8.22 (11)a27.825.12 (13)2.86 (−7)Ocinte→Teinteb11.142.50 (13)1.90 (5)c11.171.75 (13)1.32 (5)a7.077.80 (13)6.16 (8)Teinte→Ocinteb1.522.06 (13)1.60 (12)c1.471.53 (13)1.36 (12)a: The results for a rigid lattice.b: The results for a classical lattice.c: The results for a quantized lattice.Figures11 and 12 clearly show that the classical lattices always reduce the free energy barriers compared with the rigid lattices, however, the differences between the free energy barriers with the quantized lattices and the ones with the classical lattices are very small. For different processes; however, the relaxation effect on the free energies is very different. For hydrogen diffusion on Ni(111), the classical lattice only slightly lowers the free energy barrier, while it decreases the barriers by more than one-half in subsurface and interior processes. More careful analysis from Figure 11 reveals that the two preferred Ni(111) surface binding sites, that is, hcp and fcc hollows have symmetric wells, manifesting that the motions of Ni atoms in the layer beneath the surface have little influence on the surface free energies despite the fact that the Ni atoms beneath the hcp and fcc hollows have different arrangements. It is also found that the Ocsub site has a deeper well than the Tesub site and the well of the Tesub site nearly disappears as relaxed Ni atoms are considered. It manifests that the hydrogen at the Tesub site is very unstable and can easily move to the Ocsub site or resurface to the hcp site. In the interior processes, although the hydrogen in the Teinte site is much less stable than in the Ocinte site (see Figure 12), the well at the Teinte is obvious. One thus expects that a two-step reaction process can be used for the reaction from one Ocinte to the other.Next, we consider the prefactor. Table5 shows that the prefactors of the classical lattices are much smaller than those of rigid ones except for the process of H diffusion on Ni(111), and the prefactors of quantized lattices are always smaller than that of classical ones, but their amplitudes have the same order. This may be explained by the fact that the quantum motions of the lattice atoms can induce the dissipative effect on the tunneling degrees of freedom [62, 63], because the prefactor essentially incorporates the dynamical effect. It is well known that pure dissipation in the overdamping regime always hinders the reaction rates for a given reaction barrier. The present results thus are consistent with above analysis.The rate is determined by both the prefactor and the free energy barrier. In Table5, the rate with the rigid lattice is smaller than the one on the classical lattice, which is because the prefactor changes a little and the free energy determines the rate. However, the rate with the quantized lattice is smaller than that of the classical lattice, which is due to the fact that quantized lattice has a smaller prefactor while the free energies are similar. Generally speaking, for the quantized lattices, the rates are lower by 20%–40% when compared to the ones on the classical lattices.Another important feature of the free energy is its temperature dependence. Figure13 displays the free energy profiles with both H- and Ni-treated quantum mechanically at several temperatures. Generally speaking, the free energies have a slight difference at 300 K and 400 K, whereas this difference becomes pronounced for 100 K and 200 K, and the barrier positions move to the directions of shallow well for asymmetric reactions. These properties can be explained by the hydrogen tunneling effect. At lower temperatures, the tunneling plays a more important role. Indeed, the barrier heights decrease with decreasing temperature except for Ocsub→Tesub among 200 K to 400 K. This special case may be due to the special structure of the lattice. The thermal average displacements of Ni1 and Ni3 (in Figure 9) vertical to the Ni(111) surface increase with increasing temperature. Thus, H goes through reaction bottleneck easier at a higher temperature, which makes the barrier decrease. Figure 13 also displays that the free energy barriers of Ocinte→Teinte change little in the temperature range of 200–400 K, while those of the Teinte→Ocinte become smaller and smaller with decreasing temperature. The corresponding free energy barriers are 0.47, 1.09, 1.47, and 1.66 kcal/mol at 100, 200, 300, and 400 K, respectively. It manifests that the diffusing H atom may not equilibrate in the interior tetrahedral vacancy at very low temperatures.The temperature dependence of the free energies. The blue, the black, the red, and the green lines correspond to free energy profiles at 100 K, 200 K, 300 K, and 400 K, respectively.
(a)(b)(c)(d)


### 3.4.3. Surface-Subsurface Transport

In the resurfacing process,Tesub→Ahcp does not show an obvious barrier (see Figure 11), as the lattice atoms are treated quantum mechanically. This step thus can be thought as a barrierless process. Table 6 tabulates the rate constants for the other resurfacing and subsurface processes in the temperature range of 100 to 600 K. We also list the available CVT/SCT results for Ocsub→Afcc and its reverse reaction [60]. Again, the CVT/SCT rates are close to the QI rates except at 100 K. We think our much bigger rate constant at 100 K is due to the contribution of tunneling.Table 6
Rate constants (Unit:s-1. Powers of 10 are in parentheses) for surface-subsurface processes. aThe QI results with a quantized lattice. bThe CVT/SCT results from [60].
100 K200 K300 K400 K600 KAfcc→Ocsuba3.81 (−14)1.14 (−1)2.23 (3)5.48 (5)2.43 (8)Afcc→Ocsubb1.09 (−15)1.23 (−1)9.72 (3)3.36 (6)1.16 (9)Ocsub→Afcca1.17 (6)4.59 (7)2.15 (9)1.23 (10)1.02 (11)Ocsub→Afccb1.62 (5)7.36 (7)2.93 (9)2.33 (10)2.15 (11)Ocsub→Tesuba2.78 (−4)1.99 (4)2.74 (7)1.18 (9)2.49 (10)Tesub→Ocsuba7.24 (8)4.68 (10)8.22 (11)2.00 (12)4.13 (12)


### 3.4.4. Interior Migration

In the interior of bulk Ni, the two most stable sites to cage H are symmetric octahedral vacancies (see Figure12). H diffusion between them has been measured experimentally [64–67]. Several theoretical calculations have also been proposed to investigate this diffusion process. Wimmer et al. [68] calculate the diffusion coefficient via two-step reactions Ocinte→Teinte→Ocinte by using a transition state theory together with accurate ab initio energies, while Wonchoba and Truhlar [60] consider the kinetic step as a direct process with a double maximum barrier and calculate the diffusion coefficient using the CVT/SCT. In the present QI calculations, Figure 13 has explicitly shown that the free energies have a well at the tetrahedral site from 200 K to 400 K. It is thus reasonable to assume that the diffusing H atom temporarily equilibrates in the tetrahedral site before jumping forward or backward to a neighboring octahedral site. The free energy well, however, becomes very shallow at 100 K. In this case, the direct reaction from the octahedral site to the other one may be acceptable. Here, we only calculate the diffusion coefficients via the two kinetic steps at 200–400 K.The temperature dependence of the diffusion coefficients are commonly fitted to the Arrhenius equation(27)D(T)=D0⋅exp⁡⁡[-EaRT],
where R is the gas constant and D0 and Ea are the pre-exponential factor and the activation energy, respectively. The QI calculations predict D0=3.93×10-3cm2s-1 and Ea=10.26kcal/mol. In the calculations, the rates of Ocinte→Teinte are used to obtain the diffusion coefficients for Ocinte→Ocinte, because this process is much slower than that of Teinte→Ocinte and it determines the total reaction rates.Table7 tabulates the pre-exponential factors and activation energies coming from available experiments and theories. It is found that both D0 and Ea from the QI calculations are close to Ebisuzaki’s experimental data [65]. Further tracking down the comparisons with experiments is nontrivial, because the accuracy of the diffusion coefficients is much dependent of the potential energy surface. However, we may make a quantitative comparison for the QI and CVT/SCT results. With use of the diffusion coefficients from 295 K to 300 K obtained by the CVT/SCT, Wonchoba and Truhlar [60] predict 11.1 kcal/mol for Ea, and 1.3×10-3cm2s-1 for D0, respectively. These values are observably different from the QI calculations. Ea and D0 are 0.8 kcal/mol larger and 3 times smaller than those from the QI calculations, respectively. The origin of these differences can be explained by that Wonchoba and Truhlar have treated the processes Ocinte→Teinte and Teinte→Ocinte as a single kinetic step rather than as two kinetic steps, which results in a much longer tunneling path than that of the two steps. Compared with the results reported by Wimmer et al. [68], the activation energy is about 0.6 kcal/mol larger than the present one, and the pre-exponential factor is 10 times larger, which is also larger than all available experimental data. These differences may come from both the different potential energy surfaces and rate methods.Table 7
Arrhenius parameters for H diffusion in interior of bulk Ni.aThe QI diffusion coefficients for Ocinte→Ocinte. bFrom [60]. cFrom [68]. dFrom [64]. eFrom [65]. fFrom [66]. gFrom [67].
D0(cm2s-1)Ea(kcal/mol)Present QIa200–400 K3.93 (−3)10.26295–300 K1.3 (−3)11.1Wonchobab (theoretical data)300–627 K7.8 (−4)10.9627–1650 K4.4 (−4)10.2Wimmerc (theoretical data)273–1000 K3.84 (−2)10.92Yamakawad (experimental data)220–350 K1.9 (−3)8.89Ebisuzakie (experimental data)470–690 K5.22 (−3)9.56Eichenauerf (experimental data)660–930 K6.73 (−3)9.47Katzg  (experimental data)670–1270 K7.04 (−3)9.43


## 3.4.1. Model

Figures9 and 10 show H diffusion processes in the nickel crystal with a face-centered-cubic (fcc) lattice structure. For a given rate process, we construct a lattice cell (each cell contains more than 200 Ni atoms) in which all atoms are treated to be movable to incorporate the effect of the crystal fluctuation on the rates. To be concrete, the Ni atoms in the four sides of the bulk metal and at the bottom layers are fixed, the Ni atoms surrounding the reactant and product sites and lying along the reaction path are treated quantum mechanically, and the others are treated classically. It should be noted that the lattice for the interior diffusion process has a structure of sphere and the outer layers are fixed.Figure 9
A lattice model with a few Ni atoms for the circle reaction processes of hydrogen.Ahcp,Afcc,Ocsub, and Tesub are the abbreviations for a hcp hollow site, a fcc hollow site, a subsurface octahedral vacancy, and a subsurface tetrahedral vacancy, respectively. The black lines stand for the directions of the reaction paths. The atoms of Ni1,Ni2, and Ni3 colored in orange are specially used to determine the reaction coordinate of the process Ocsub→Tesub.Figure 10
A lattice model with a few Ni atoms for the process of H diffusion in interior of bulk Ni. Ocinte and Teinte are the abbreviations for an interior octahedral vacancy and an interior tetrahedral vacancy, respectively. The black lines stand for the general directions of the reaction paths. The atoms of Ni4,Ni5, and Ni6 colored in orange are specially used to determine the reaction coordinate of the process Ocinte→Teinte.In the QI calculations, we need to define the reaction coordinate operatorsŝ  (6). For the systems considered in this section, the hydrogen coordinates are essentially good choices. We thus adopt the following reaction coordinates for different rate processes. In the H diffusion on Ni(111) surface from a hcp site to a fcc site (Ahcp→Afcc), the x coordinate of the H atom is chosen, whose direction is showed in Figure 9 as the black line connecting the Ahcp site to Afcc site. In the H resurfacing from a subsurface octahedral vacancy to the fcc site (Ocsub→Afcc) and from a subsurface tetrahedral vacancy to the hcp site (Tesub→Ahcp), the z coordinate of the H atom (vertical to Ni(111) surface) is taken. In the H diffusion between the adjacent subsurface octahedral and tetrahedral vacancies (Ocsub→Tesub) and between the adjacent interior octahedral vacancy and interior tetrahedral vacancy (Ocinte→Teinte), the reaction coordinates are along the directions that perpendicular to the planes of Ni1–Ni2–Ni3 and Ni4–Ni5–Ni6, respectively.


## 3.4.2. Free Energy and Prefactor

In order to investigate the quantized Ni lattice effect on the rates, we have recast the QI formula in (1) into Arrhenius form, which consists of the free energy and the prefactor (17) and (18).In the free energy calculations, firstly, we calculate the free energy profile for each step of the H hopping paths,Ahcp→Afcc→Ocsub→Tesub→Ahcp and Ocinte→Teinte→Ocinte, with the reaction coordinate defined in Section 3.4.1. Then, we connect these free energy profiles one by one, the final free energy profiles for the whole processes are displayed in Figures 11 and 12.Figure 11
Free energy profiles with respect to the processes in Figure9 at 300 K. The green, the red, and the blue lines correspond to the rigid, classical, and quantized lattices, respectively.Figure 12
Free energy profiles with respect to the processes in Figure10 at 300 K. The green, the red, and the blue lines correspond to the rigid, classical, and quantized lattices, respectively.Figures11 and 12 display the calculated free energy profiles for the surface-subsurface and interior processes, respectively, with the quantized lattice, the classical lattice and the rigid one at room temperature (300 K). The corresponding free energy barriers, prefactors, and rates are tabulated in Table 5.Table 5
Free energies, prefactors, and rate constants.
Process  (300 K)ΔF  (kcal/mol)AQI (s-1)kQI (s-1)a2.191.58 (13)4.07 (11)Afcc→Ahcpb2.141.61 (13)4.49 (11)c2.141.27 (13)3.47 (11)a33.066.87 (13)6.22 (−11)Afcc→Ocsubb13.792.85 (13)2.64 (3)c13.752.29 (13)2.23 (3)a20.006.44 (13)1.83 (−1)Ocsub→Afccb5.402.10 (13)2.62 (9)c5.401.72 (13)2.15 (9)a34.924.29 (13)1.67 (−12)Ahcp→Tesubb14.92c14.87Tesub→Ahcpa3.255.13 (13)2.28 (11)a26.195.23 (13)4.31 (−6)Ocsub→Tesubb7.831.86 (13)3.51 (7)c7.861.45 (13)2.74 (7)a7.838.17 (13)1.51 (8)Tesub→Ocsubb1.751.91 (13)1.05 (12)c1.731.49 (13)8.22 (11)a27.825.12 (13)2.86 (−7)Ocinte→Teinteb11.142.50 (13)1.90 (5)c11.171.75 (13)1.32 (5)a7.077.80 (13)6.16 (8)Teinte→Ocinteb1.522.06 (13)1.60 (12)c1.471.53 (13)1.36 (12)a: The results for a rigid lattice.b: The results for a classical lattice.c: The results for a quantized lattice.Figures11 and 12 clearly show that the classical lattices always reduce the free energy barriers compared with the rigid lattices, however, the differences between the free energy barriers with the quantized lattices and the ones with the classical lattices are very small. For different processes; however, the relaxation effect on the free energies is very different. For hydrogen diffusion on Ni(111), the classical lattice only slightly lowers the free energy barrier, while it decreases the barriers by more than one-half in subsurface and interior processes. More careful analysis from Figure 11 reveals that the two preferred Ni(111) surface binding sites, that is, hcp and fcc hollows have symmetric wells, manifesting that the motions of Ni atoms in the layer beneath the surface have little influence on the surface free energies despite the fact that the Ni atoms beneath the hcp and fcc hollows have different arrangements. It is also found that the Ocsub site has a deeper well than the Tesub site and the well of the Tesub site nearly disappears as relaxed Ni atoms are considered. It manifests that the hydrogen at the Tesub site is very unstable and can easily move to the Ocsub site or resurface to the hcp site. In the interior processes, although the hydrogen in the Teinte site is much less stable than in the Ocinte site (see Figure 12), the well at the Teinte is obvious. One thus expects that a two-step reaction process can be used for the reaction from one Ocinte to the other.Next, we consider the prefactor. Table5 shows that the prefactors of the classical lattices are much smaller than those of rigid ones except for the process of H diffusion on Ni(111), and the prefactors of quantized lattices are always smaller than that of classical ones, but their amplitudes have the same order. This may be explained by the fact that the quantum motions of the lattice atoms can induce the dissipative effect on the tunneling degrees of freedom [62, 63], because the prefactor essentially incorporates the dynamical effect. It is well known that pure dissipation in the overdamping regime always hinders the reaction rates for a given reaction barrier. The present results thus are consistent with above analysis.The rate is determined by both the prefactor and the free energy barrier. In Table5, the rate with the rigid lattice is smaller than the one on the classical lattice, which is because the prefactor changes a little and the free energy determines the rate. However, the rate with the quantized lattice is smaller than that of the classical lattice, which is due to the fact that quantized lattice has a smaller prefactor while the free energies are similar. Generally speaking, for the quantized lattices, the rates are lower by 20%–40% when compared to the ones on the classical lattices.Another important feature of the free energy is its temperature dependence. Figure13 displays the free energy profiles with both H- and Ni-treated quantum mechanically at several temperatures. Generally speaking, the free energies have a slight difference at 300 K and 400 K, whereas this difference becomes pronounced for 100 K and 200 K, and the barrier positions move to the directions of shallow well for asymmetric reactions. These properties can be explained by the hydrogen tunneling effect. At lower temperatures, the tunneling plays a more important role. Indeed, the barrier heights decrease with decreasing temperature except for Ocsub→Tesub among 200 K to 400 K. This special case may be due to the special structure of the lattice. The thermal average displacements of Ni1 and Ni3 (in Figure 9) vertical to the Ni(111) surface increase with increasing temperature. Thus, H goes through reaction bottleneck easier at a higher temperature, which makes the barrier decrease. Figure 13 also displays that the free energy barriers of Ocinte→Teinte change little in the temperature range of 200–400 K, while those of the Teinte→Ocinte become smaller and smaller with decreasing temperature. The corresponding free energy barriers are 0.47, 1.09, 1.47, and 1.66 kcal/mol at 100, 200, 300, and 400 K, respectively. It manifests that the diffusing H atom may not equilibrate in the interior tetrahedral vacancy at very low temperatures.The temperature dependence of the free energies. The blue, the black, the red, and the green lines correspond to free energy profiles at 100 K, 200 K, 300 K, and 400 K, respectively.
(a)(b)(c)(d)


## 3.4.3. Surface-Subsurface Transport

In the resurfacing process,Tesub→Ahcp does not show an obvious barrier (see Figure 11), as the lattice atoms are treated quantum mechanically. This step thus can be thought as a barrierless process. Table 6 tabulates the rate constants for the other resurfacing and subsurface processes in the temperature range of 100 to 600 K. We also list the available CVT/SCT results for Ocsub→Afcc and its reverse reaction [60]. Again, the CVT/SCT rates are close to the QI rates except at 100 K. We think our much bigger rate constant at 100 K is due to the contribution of tunneling.Table 6
Rate constants (Unit:s-1. Powers of 10 are in parentheses) for surface-subsurface processes. aThe QI results with a quantized lattice. bThe CVT/SCT results from [60].
100 K200 K300 K400 K600 KAfcc→Ocsuba3.81 (−14)1.14 (−1)2.23 (3)5.48 (5)2.43 (8)Afcc→Ocsubb1.09 (−15)1.23 (−1)9.72 (3)3.36 (6)1.16 (9)Ocsub→Afcca1.17 (6)4.59 (7)2.15 (9)1.23 (10)1.02 (11)Ocsub→Afccb1.62 (5)7.36 (7)2.93 (9)2.33 (10)2.15 (11)Ocsub→Tesuba2.78 (−4)1.99 (4)2.74 (7)1.18 (9)2.49 (10)Tesub→Ocsuba7.24 (8)4.68 (10)8.22 (11)2.00 (12)4.13 (12)


## 3.4.4. Interior Migration

In the interior of bulk Ni, the two most stable sites to cage H are symmetric octahedral vacancies (see Figure12). H diffusion between them has been measured experimentally [64–67]. Several theoretical calculations have also been proposed to investigate this diffusion process. Wimmer et al. [68] calculate the diffusion coefficient via two-step reactions Ocinte→Teinte→Ocinte by using a transition state theory together with accurate ab initio energies, while Wonchoba and Truhlar [60] consider the kinetic step as a direct process with a double maximum barrier and calculate the diffusion coefficient using the CVT/SCT. In the present QI calculations, Figure 13 has explicitly shown that the free energies have a well at the tetrahedral site from 200 K to 400 K. It is thus reasonable to assume that the diffusing H atom temporarily equilibrates in the tetrahedral site before jumping forward or backward to a neighboring octahedral site. The free energy well, however, becomes very shallow at 100 K. In this case, the direct reaction from the octahedral site to the other one may be acceptable. Here, we only calculate the diffusion coefficients via the two kinetic steps at 200–400 K.The temperature dependence of the diffusion coefficients are commonly fitted to the Arrhenius equation(27)D(T)=D0⋅exp⁡⁡[-EaRT],
where R is the gas constant and D0 and Ea are the pre-exponential factor and the activation energy, respectively. The QI calculations predict D0=3.93×10-3cm2s-1 and Ea=10.26kcal/mol. In the calculations, the rates of Ocinte→Teinte are used to obtain the diffusion coefficients for Ocinte→Ocinte, because this process is much slower than that of Teinte→Ocinte and it determines the total reaction rates.Table7 tabulates the pre-exponential factors and activation energies coming from available experiments and theories. It is found that both D0 and Ea from the QI calculations are close to Ebisuzaki’s experimental data [65]. Further tracking down the comparisons with experiments is nontrivial, because the accuracy of the diffusion coefficients is much dependent of the potential energy surface. However, we may make a quantitative comparison for the QI and CVT/SCT results. With use of the diffusion coefficients from 295 K to 300 K obtained by the CVT/SCT, Wonchoba and Truhlar [60] predict 11.1 kcal/mol for Ea, and 1.3×10-3cm2s-1 for D0, respectively. These values are observably different from the QI calculations. Ea and D0 are 0.8 kcal/mol larger and 3 times smaller than those from the QI calculations, respectively. The origin of these differences can be explained by that Wonchoba and Truhlar have treated the processes Ocinte→Teinte and Teinte→Ocinte as a single kinetic step rather than as two kinetic steps, which results in a much longer tunneling path than that of the two steps. Compared with the results reported by Wimmer et al. [68], the activation energy is about 0.6 kcal/mol larger than the present one, and the pre-exponential factor is 10 times larger, which is also larger than all available experimental data. These differences may come from both the different potential energy surfaces and rate methods.Table 7
Arrhenius parameters for H diffusion in interior of bulk Ni.aThe QI diffusion coefficients for Ocinte→Ocinte. bFrom [60]. cFrom [68]. dFrom [64]. eFrom [65]. fFrom [66]. gFrom [67].
D0(cm2s-1)Ea(kcal/mol)Present QIa200–400 K3.93 (−3)10.26295–300 K1.3 (−3)11.1Wonchobab (theoretical data)300–627 K7.8 (−4)10.9627–1650 K4.4 (−4)10.2Wimmerc (theoretical data)273–1000 K3.84 (−2)10.92Yamakawad (experimental data)220–350 K1.9 (−3)8.89Ebisuzakie (experimental data)470–690 K5.22 (−3)9.56Eichenauerf (experimental data)660–930 K6.73 (−3)9.47Katzg  (experimental data)670–1270 K7.04 (−3)9.43


## 4. Conclusion

We have presented the basic principle of the quantum instanton (QI) approximation and its applications to chemical reactions from gas phase to surface systems. The applications demonstrate that the QI method makes it possible to treat hundreds of atoms because of the well-established techniques of imaginary time path integrals. For instance, more than 200 atoms have been incorporated in the present study of the H diffusion processes.The QI approximation is a kind of ‘‘quantum transition state theory’’ in that there is no account of “recrossing” dynamics in the description. The recrossing effects on the quantum instanton rate constants have been quantified for several collinear reaction by Ceotto and Miller [25], and it is found especially evident for the collinear heavy-light-heavy reactions. Fortunately, the recrossing effects become generally less important in higher dimensions [69]. Therefore, the QI approach may become a suitable tool for the calculation of chemical reaction rate constants of complex systems.Compared to conventional TST theories, the QI approximation involves two dividing surfaces, which are quantum analogs of the two turning point surfaces of the imaginary time trajectory in the semiclassical instanton theory. At high-temperature limit, these dividing surfaces coalesce into the one, the same as the dividing surface from Wigner’s variational principle. In this case, the QI approximation becomes exactly the same as the classical TST. As the tunneling corrections are incorporated, the CVT rates are consistent with the QI rates except at deep tunneling regime, where QI rates are generally greater than the rates from CVT with various tunneling corrections, since the CVT method uses an optimized tunneling path, while the QI method considers all tunneling paths and automatically gives each path its natural weight by the quantum Boltzmann factor. For the reaction ofH+CH4→H2+CH3, accurate quantum dynamics rate constants are obtained with the MCTDH method. Compared to the MCTDH ones, the QI rate constants are larger by factors of about 2 to 3 over the temperature range 300–400 K. This difference may partly be due to the recrossing effect which is not considered in QI theory and partly arise from the use of the J-shifting approximation and the neglect of the vibrational angular momenta Hamiltonian in MCTDH method.



---

*Source: 483504-2011-12-11.xml*