ℏ 2 expansion of the transmission probability through a barrier

Ninety years ago, Wigner derived the leading order expansion term in h 2 for the tunneling rate through a symmetric barrier. His derivation included two contributions: one came from the parabolic barrier, but a second term involved the fourth-order derivative of the potential at the barrier top. He left us with a challenge, which is answered in this paper, to derive the same but for an asymmetric barrier. A crucial element of the derivation is obtaining the h 2 expansion term for the projection operator, which appears in the flux-side expression for the rate. It is also reassuring that an analytical calculation of semiclassical transition state theory (TST) reproduces the anharmonic corrections to the leading order of h 2 . The efficacy of the resulting expression is demonstrated for an Eckart barrier, leading to the conclusion that especially when considering heavy atom tunneling, one should use the expansion derived in this paper, rather than the parabolic barrier approximation. The rate expression derived here reveals how the classical TST limit is approached as a function of h and, thus, provides critical insights to understand the validity of popular approximate theories, such as the classical Wigner, centroid molecular dynamics, and ring polymer molecular dynamics methods.

where Vn is the nth derivative of the potential at the barrier top and ω is the (imaginary) barrier top frequency (V 2 = −Mω 2 ). The first term is well known, and it is the leading order expansion term for the tunneling probability through a parabolic barrier, as also derived by Wigner. The second term, which includes the leading order nonlinear term in the (symmetric) potential, is also not appreciated in the literature.
In his 1932 paper, Wigner noted that in what at the present time would be called the classical Wigner dynamics approximation, 2 the non-linear term is missing. In other words, when computing the rate using the correct thermal density of up to order h 23 but using a classical projection operator, 4,5 one regains only the parabolic barrier term. At present, any approximate thermal rate theory is considered to be valid, provided that in the high temperature limit, it reduces to the parabolic barrier estimate. However, do these theories reduce to the correct high temperature limit? As may be inferred from Eq. (1.1), the high temperature limit is linear in β, not quadratic as for the parabolic barrier.
Wigner also left us with a challenge: What is the correct leading order term in h 2 when the potential is not symmetric? The central theme of this paper is to first of all derive the leading order expansion term for any bell-shaped potential barrier. We will show that the recent semiclassical theory developed in Refs. 6-9 recovers the correct leading order term. However, all other popular approximations do not; these include classical Wigner dynamics, [10][11][12] centroid molecular dynamics, 13 and ring polymer molecular dynamics. 14 They approach the parabolic barrier limit differently from the correct dependence in the presence of nonlinearity.
Wigner himself after a remarkable and lengthy derivation notes that he does not expect that the term that goes as V 4 will be large The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp or important. Why, then, should one go at the present time to the considerable effort needed to rederive it and generalize it to the asymmetric case? The definition of the words "large" or "important" is vague and in the eyes of the beholder. At this point, suffice it to say that the theory has been used by various researchers in recent years to estimate rate constants, especially for heavy atom tunneling reactions 15,16 (Wigner usually considered only H atom transfer), where sometimes the parabolic barrier expression is employed. However, aside from practical issues, we note that in the symmetric case, as seen from Eq. (1.1), assuming as is usually the case that V 4 ≥ 0, the nonlinearity serves to increase the rate at high temperature, while based on our semiclassical intuition and the knowledge that the tunneling probability decreases exponentially with the action of the appropriate instanton trajectory, we know that typically nonlinearity increases the action and lowers the rate. Is the same true in the asymmetric case? We also note that in a recent paper 17 (which also motivated the present study), we showed how the leading order terms in h 2 and h 4 may be used to obtain an approximate expression for the thermal rate constant even for temperatures below the well-known crossover temperature between thermal activation and deep tunneling. 18,19 Development of a general high temperature theory for tunneling in multidimensional systems is, thus, also a challenge, whose solution has both theoretical and practical implications.
The central ingredient that also enables the derivation for asymmetric potentials is the formal h 2 expansion of the quantum expression for the projection operator onto products, derived in Ref. 4. It is combined with the h 2 expansion of the thermal density, obtained by Wigner in 1932, 3 and the well-understood formula for expressing the Wigner representation of a product of two operators in terms of the known Wigner representation of each operator separately. 20 The resulting expression is Noting that V 2 ≤ 0, we see that both the leading order asymmetric and symmetric anharmonic terms typically lead to an increase of the transmission factor, as compared to the parabolic barrier rate. The theory is presented in Sec. II. The demonstration that the semiclassical theory also known as the VPT2 semiclassical transition state theory (TST) 7 has the correct leading order term is presented in Sec. III. A brief analysis of the high temperature limits of other approximate theories is also presented in Sec. III. In Sec. IV, we present some numerical computations on the symmetric and asymmetric Eckart barrier 21 so as to give a perception for the magnitude of the effect of the nonlinearity on the transmission coefficient. We end with a discussion of the implications of the theory and further possible generalizations.

A. Background
We consider a one-dimensional Hamiltonian operator The potential is assumed to have a barrier at the coordinate location q = 0 with energy V(0) = 0, and (p,q) denote the momentum and coordinate operators, respectively. The flux operator, which measures the flux at the barrier top, is defined aŝ where δ(x) denotes the Dirac "delta" function. The transmission coefficient is defined as the ratio of the exact quantum thermal flux through the barrier to the classical thermal flux, 22,23 Here, the projection operator onto the product side is defined as 23 and θ(q) is the unit step function. The Wigner representation in phase space of an operatorÔ is defined as 3 We will use two well-known properties of Wigner representations. One is that the Trace operation for the product of two operators is just the phase space integration of the product of the separate Wigner representations of the two operators, The second is that the Wigner representation of a product of operatorsĈ may be expressed in terms of the representation of the separate operators through the relation 20,24 where the operatorΛ is defined aŝ The right and left arrows denote derivatives of the functions to the right and left of the operator, respectively. We will need to employ Eq. (2.8) since the expression for the transmission coefficient is the trace of a product of three operators.

B. h − 2 expansions
The exact Wigner representation of the flux operator has also its classical form, (2.10)

ARTICLE scitation.org/journal/jcp
We will also assume formal expansions for the thermal density operator, the projection operator, 12) and the product of the flux operator and the projection operator, With these notations, the expression for the transmission factor is given by (2.14) Wigner showed in 1932 3 that for the density operator, and ). (2.16) Using the Taylor expansion of the exponential function and the fact that the Wigner representation of the flux operator is known exactly [Eq. (2.10)] and is identical to the classical flux expression, we find from Eqs. (2.8) and (2.9) that the leading order terms of the Wigner representation of the product of the flux and projection operators are given by F P 0 (p, q) = F(p, q)P 0 (p, q), (2.17) We also know that the zeroth order projection operator is just the classical projection operator, The more challenging part is to write down the first order contribution P 1 (p, q) for the projection operator.
where primes denote derivatives with respect to the argument and δ(x) denotes the Dirac "delta" function. With these preliminaries, we can derive the leading order contribution to transmission factor.

C. Evaluation of the transmission factor
From Eq. (2.14), we find that the transmission factor may be written as One readily finds that and as also shown in Ref This is just the leading order correction for the parabolic barrier transmission coefficient. As detailed in Appendix A, some straightforward algebra shows that The more elaborate part is the remaining term As already noted, the flux F(p, q) localizes everything to the barrier location q = 0, which means that to estimate this factor, it is necessary to evaluate the functions g ± 21 (q), g ± 11 (q), and g ± 01 (q) in the limit that q → 0. For this purpose, we rewrite the potential as

29)
The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp and as a matter of notation, we have that the second order derivative is It is then a matter of some straightforward but tedious algebra to find that One, then, readily finds that so that putting it all together, we find the expression and this is the central result of this paper. Wigner's result for the symmetric barrier 1 has been now generalized to any (smooth) barrier.
Note that this result has been derived without using any knowledge about the asymptotic structure of the potential barrier. It is correct as long as the reactive flux expression for the rate is valid and includes in it only the local structure around the barrier top.

A. Semiclassical theory
The semiclassical transition state theory developed in Refs. 6-9 is based on the WKB approximation for the energy dependent transmission probability. 25,26 The classical action associated with quantum tunneling is The integral is along the constant energy trajectory between the turning points on the upside down potential. The corresponding semiclassical energy-dependent transmission coefficient is given as 25,26 TSC(E) = 1 , (3.2) and the thermal semiclassical transmission probability (normalized by the classical) is where V is the potential at the barrier top. If we take the classical limit, i.e., As an example, we consider a parabolic potential specified by where ω is the (imaginary) frequency at the barrier top. Then, Eq. (3.4) leads to In general, the energy associated with the action may be expanded for an anharmonic barrier as 7 where χ is expressed in terms of the anharmonicity coefficients of Eq. (2.29) as and E 0 is a constant energy defined by The corresponding transmission probability is given as where κ pb (b) is the transmission factor for the parabolic barrier, as given in Eq. (3.5). Here, ⟨⋅ ⋅ ⋅⟩ pb denotes an average using the parabolic barrier action as a reference such that to leading order in ̵ h 2 , the anharmonic corrections are obtained by the second moment.
Using κ pb (b) as a moment generating function, we find that When ̵ h → 0, this reduces to lim̵ h→0 ⟨( ϕ π ) 2 ⟩ pb = 1 12 so that the quantum correction becomes where the parabolic contribution from κ pb (b) is also truncated to leading order h 2 . Finally, we insert the values of χ and E 0 into the above expressions to arrive at  27 The relationship between the early work and the current analysis is discussed in Sec. V.

B. Classical Wigner dynamics
Classical Wigner dynamics as applied to rate theory is an approximation by which the thermal flux operator is evaluated exactly, but the projection operator is taken to be the classical projection operator. 4 This implies that in the analysis presented in Sec. II, the classical Wigner estimate is obtained by repeating the same computation but ignoring the quantum correction (P 1 (p, q)) to the projection operator. However, as we have already seen [Eq. (2.26)], this only leads to the parabolic barrier term, and the nonlinear terms arise from the quantum correction to the projection operator. Clearly, then, as also noted by Wigner, 1 the classical Wigner estimate does not give the correct answer for the leading order term in ̵ h 2 in the presence of nonlinearity. This also means that it does not give the correct high temperature leading order term.

C. Path integral based approaches 1. Variational centroid TST
Instead of using the flux side correlation function, one may also derive quantum reaction rates from what is known as the ImF formalism, 18,28 whereby one estimates the rates from the imaginary part of the partition function or equivalently the imaginary free energy. This concept originates from the analysis of the non-Hermitian Hamiltonian representation of the decay rate of a quantum state in a metastable well. At high temperature, centroid transition state theory 29 can be recovered within this framework. 19 Furthermore, anharmonicity corrections to the parabolic barrier rate can be incorporated by variational determination of an optimal free energy. Yet, the anharmonicity correction thus obtained is thermodynamic in nature, while the h 2 -expansion term derived in the previous section incorporates both thermodynamic and dynamic quantum corrections. Therefore, as we shall show, the anharmonic corrections derived from these two approaches are not necessarily equivalent and can be systematically examined by comparing the cubic and quartic contributions. Below, we outline the anharmonic corrections based on the variational ImF calculation and refer to the details in Appendix B.
In quantum TST, the centroid degree of freedom is identified as the unstable mode for the imaginary free energy. In the classical regime, the centroid-constrained quantum variance δ takes the value of a free particle, i.e., ⟨δ⟩ 2 c = ̵ h 2 β 12m . For the anharmonic potential defined in Eq. (2.29), the centroid density can be evaluated analytically using a cumulant expansion, where only the first order cumulant is kept. Evidently, the leading order correction associated with the barrier frequency recovers the parabolic barrier result, while any anharmonicity appears as higher order corrections. The TST calculation of Eq. (3.13) is based on a local expansion of the barrier potential, which sensitively depends on the location of the transition state. The free energy integrates all degrees of freedom and is, therefore, a global quantity. For a better representation of the free energy of an anharmonic potential, we adopt a reference quadratic potential V ref centered at the as yet arbitrary locationq with frequency ω. These are then optimized by a variational calculation of the free energy. Details can be found in Appendix B, which closely follows the theory presented in Refs. 30 and 31. Using the optimal quadratic potential, the centroid density is again evaluated via a cumulant expansion as in Eq. (3.13), giving (3.14) The anharmonic corrections in the above expression have the same functional form as Eq. (2.35) but different coefficients. The systematic expansion in Eq. (2.14) involves contributions from both the thermal part Ωn and the projection operator Pn. In comparison, the imaginary free energy approach is based on the thermal distribution and the flux operator in the imaginary time. The difference in the coefficients may be thus traced to the projection operator in the reactive flux.

Implications for path integral molecular dynamics
At high temperatures, centroid molecular dynamics 32,33 and ring polymer molecular dynamics 34,35 reduce to the same pathintegral TST estimate [Eq. (3.13)] so that these theories are not precise in this limit. They do reduce to the parabolic barrier expression without the nonlinear correction terms. Inspection of the various numerical results obtained with these theories shows that at high

IV. NUMERICAL OBSERVATIONS USING ECKART BARRIERS
The purpose of this section is to obtain a perception for the extent of the nonlinear correction. For this purpose, we study the asymmetric Eckart barrier potential, which has the form . (4.1) The potential is constructed such that when x → −∞, the potential vanishes and when x → ∞, the potential goes to V 1 − V 2 . The potential is symmetric when V 1 = V 2 . Here, V 1 is the barrier height, which is located at We use the reduced variables such that α is the reduced value of h and v is the asymmetry parameter of the potential. The barrier frequency ω is defined as (M is the mass of the scattered particle) The reduced force constant at the barrier is found to be (using Maple) and we note that it varies between − 1 2 and −2 as the asymmetry parameter goes from 1 to ∞. This implies that the barrier gets thinner as the asymmetry is increased.
With these preliminaries, one may write the (reduced) energy dependent transmission probability as , (4.6) and the thermal transmission probability normalized by the classical transmission probability is given by One, then, notes that the thermal transmission probability is a function of three variables only-α, v, and βV 1 .
With these variables, the leading order correction based on the parabolic barrier approximation is given by . (4.8) Using Maple, one readily finds that the fourth and third derivatives at the barrier are , (4.10) from which we deduce that the exact expansion for the transmission factor, as in Eq. (1.2), is while the path integral TST approximate expression given in Eq. (3.14) is It is, then, of interest to compare the numerically exact transmission coefficient with the exact leading order correction (κ 2 ), the parabolic barrier estimate (κ pb ), and the path integral TST estimate (κ 2,PI ). This is shown in Fig. 1 Figure 2 shows the values of ̵ hβω corresponding to the data presented in Fig. 1. One notes that for the symmetric barrier (v = 1), for all values shown, the analytical estimate is the most accurate, being essentially exact when α, the reduced value of h, is small. The path integral TST estimate is too large, however, more accurate than the parabolic barrier estimate. As one increases the asymmetry, the path integral TST estimate improves and is slightly better at larger values of α than the exact expansion. Of course, at low values of α, the exact expansion is superior. The values of hβω shown in Fig. 2 are in the range of 0 − 5, demonstrating that the leading order term in the expansion provides a good approximation to the exact transmission coefficient in the whole range and is significantly superior to the parabolic barrier estimate.
Increasing the temperature, as shown in Fig. 3, further improves the quality of the exact expansion, as might have been expected. In the symmetric case, it is essentially exact in the range 0 ≤ hβω ≤ 2, as seen in Fig. 4. However, also in the asymmetric cases, shown in the second and third panels of Fig. 3, the leading order expansion term is a much better approximation to the transmission coefficient over the whole range of values shown. The conclusion of all of this is that when considering heavy atom scattering, where typical tunneling factors are not greater than 4, one should use the second order expansion derived in this paper to estimate the rate.
On the other hand, it is also relevant to ask how much do the nonlinear terms change the transmission probability? From Eq. (4.11), we note that the nonlinear contribution term is a monotonically decreasing function of the asymmetry parameter, suggesting that the larger the asymmetry, the less important the nonlinear terms. This is readily seen from Figs. 1 and 3, where one notes that the error in the parabolic barrier estimate decreases as the asymmetry increases.

V. DISCUSSION
A central goal of this paper was to derive the general leading order expansion term of order h 2 for the thermal transmission coefficient when the potential barrier is non-quadratic. Since the flux side expression for the exact rate is a product of three operators, this exercise is rather lengthy, as shown in Sec. II. One may continue the theoretical developments presented in Sec. II to also derive the next order, except that the computation becomes increasingly involved as it includes many terms.
This paper is limited to the one-dimensional case. In principle, one may use the same formalism to derive correction terms in the multidimensional case, for example, by using the anharmonic expansions of molecular potentials, as presented in Refs. 6-9. This would also provide a check on the accuracy of the ̵ h 4 terms used in the semiclassical theories developed in Refs. 6, 8, and 9.
The theory derived in this paper presents a challenge to many present day approximate rate theories. We have shown that only the semiclassical based theory leads to the correct leading order term in ̵ h 2 . However, as far as other approximations are concerned, it is not enough to show that they reduce to the parabolic barrier result. The nonlinearity of the potential typically provides tunneling factors, which are substantially higher than the parabolic barrier rates, even when close to the classical limit.
This does not imply that the semiclassical estimate is always a "good" approximation. In Fig. 5, we plot the ratio of the numerically exact semiclassical approximation for the asymmetric Eckart barrier [Eq.  3)] and assuming a mass of a proton (∼2000 a.u.) and a barrier length scale d ∼ 2 a.u., one finds that α ∼ 3 implies that hω ∼ 10 cm −1 , which is an unrealistically broad barrier. Increasing the reduced barrier height would increase the value of the frequency roughly as the square root of the reduced barrier height, which implies that for most systems of chemical interest, unless the temperature is very low, the semiclassical theory would be valid. This is reassuring.
It is of special interest to analyze why path integral molecular dynamics (e.g., centroid molecular dynamics, ring-polymer molecular dynamics, and classical Wigner dynamics) fails in obtaining the precise leading order term in ̵ h 2 for the rate. The source of the discrepancy comes from the quantum correction to the projection operator. In a sense, a similar failure of path integral molecular dynamics has also been presented for the vibrational response of anharmonic oscillators. Centroid molecular dynamics was originally motivated by the observation that quantum and classical response functions (i.e., the Kubo-transformed correlation function) are identical for harmonic oscillators or coupled oscillators (i.e., normal modes). 36,37 The combination of this observation and the exact quantum thermal distribution provides the intuition underlying the path integral-based molecular dynamics. The fundamental frequency and beatings in harmonic systems are essentially classical coherences. A non-linear potential (i.e., anharmonic oscillator) exhibits quantum beatings, which signal a quantum coherence, which cannot be simply reproduced by continuous classical dynamics. The beating time defines the time scale at which the path integral-based molecular dynamics starts to deviate from the quantum result. A semiclassical analysis leads to an estimation of the beating time as where I is the classical action of the anharmonic oscillator and ω(I) is its frequency. 27 The beating time is proportional to the dispersion ∂ω(I) ∂I , which is determined by the anharmonicity and diverges for harmonic systems. Therefore, both in the region of the well and at the transition state, anharmonicity plays a crucial role in determining the reliability of the path integral-based molecular dynamics.
Does this invalidate the use of classical Wigner dynamics or path integral molecular dynamics methods? Not necessarily. In reality, quantum beatings or recurrences can rarely survive at long times, except when the reactants are isolated as in well-controlled gas phase chemistry. Usually, many-body interactions and dissipation will suppress the beatings and, thus, partly justify the application of the path integral molecular dynamics methods in condensed phase systems. 38 A similar argument is valid when considering the quantum expansion of the projection operator. The classical projection operator is a step function. The quantum corrections "smear" the edges of the projection operator, and it is this smearing that affects the rate. Introducing dissipation typically reduces the extent of such quantum coherences. In this context, it should be of interest to repeat the computation presented in this paper for the rate in the presence of dissipation to see whether these speculations can be given an analytical formulation.

ARTICLE scitation.org/journal/jcp
It is, then, necessary to evaluate three integrals. The first is which vanishes due to the localization at the barrier top. The second integral is given by The third integral is given by and one readily sees that the sum of the three integrals vanishes.
APPENDIX B: DERIVATION OF VARIATIONAL CENTROID TST

General considerations
Centroid transition state theory (c-TST) accounts for quantum corrections to TST above the crossover temperature so that the hightemperature expansion of c-TST provides an estimation of quantum and thermal effects. For illustration, we evaluate the leading order quantum correction. The potential surface is expanded as where Vn is the nth expansion coefficient. The centroid potential at the barrier top becomes ρc(0) = ⟨exp(−βV(δ)⟩c ≈ exp(−β⟨V(δ)⟩c), where the subscript "c" represents the centroid-constraint and δ is the displacement from the centroid. To leading order, the quadratic potential gives and this recovers the parabolic barrier rate expression to leading order in h 2 . In Eq. (B3), δ 2 is the centroid-constrained fluctuation, given as and the quadratic coefficient for the quadratic potential is Perturbative evaluation of the cubic and quartic potentials will yield h 4 or higher order corrections. Alternatively, one may estimate the effect of anharmonicity of the barrier non-perturbatively using the variational method presented below.

The variational method
The derivation here closely follows the discussion in Refs. 30 and 31. To obtain an analytic expression, one fits the exact centroid potential to a quadratic form determined by the "reference potential" The parameters (ω,q) of the reference potential are determined from the following conditions: To leading order, the first condition implies while the second one leads to At high temperature, using classical thermal fluctuations, we have that Thus, we obtainq which are, then, used to evaluate the leading order quantum correction to the centroid-TST rate.

Anharmonic corrections
To proceed, we use the variational result to evaluate the centroid correction. The optimal frequency can be directly applied to Eq. (B3), giving ρc(0) ≈ exp(−βV(0)) exp(− ̵ h 2ω 2 β 2 24 ) = exp(−βV(0)) exp(− ̵ h 2 ω 2 β 2 24 and this accounts for the correction coming from the quartic term. The optimal centroid locationq is associated with the cubic term so that one considers the centroid density at this location, The exponential part of the cubic correction can be evaluated explicitly, and it is always positive. To be consistent with the parabolic barrier notation, we rewrite the above expression as The corrections in Eqs. (B13) and (B15) are for a bounded potential, which can be converted to a reaction barrier via ω 2 → −ω 2 . Thus, these corrections have the same functional form as Eq. (2.35), but the coefficients are different, giving a larger estimate.