Semiclassical instanton theory for reaction rates at any temperature: How a rigorous real-time derivation solves the crossover temperature problem
Joseph E. Lawrence


⟨#|1<#1|
|#⟩1|#1>
#1#2⟨#2| #1 |#2⟩
#1⟶^__X#1__X


























Department of Chemistry, New York University, New York, NY 10003, USA
joseph.lawrence@nyu.edu



§ ABSTRACT 
Instanton theory relates the rate constant for tunneling through a barrier to the periodic classical trajectory on the upturned potential energy surface whose period is τ=ħ/(k_ BT). 



Unfortunately, the standard theory is only applicable below the “crossover temperature”, where the periodic orbit first appears.





This paper presents a rigorous semiclassical (ħ→0) theory for the rate that is valid at any temperature.  The theory is derived by combining Bleistein's method for generating uniform asymptotic expansions with a real-time modification of Richardson's flux-correlation function derivation of instanton theory.


The resulting theory smoothly connects the instanton result at low temperature to the parabolic correction to Eyring transition state theory at high-temperature. 






Although the derivation involves real time, the final theory only involves imaginary-time (thermal) properties, consistent with the standard theory.
Therefore, it is no more difficult to compute than the standard theory.

The theory is illustrated with application to model systems, where it is shown to give excellent numerical results.

Finally, the first-principles approach taken here results in a number of advantages over previous attempts to extend the imaginary free-energy formulation of instanton theory. In addition to producing a theory that is a smooth (continuously differentiable) function of temperature, the derivation also naturally incorporates hyperasymptotic (i.e. multi-orbit) terms, and provides a framework for further extensions of the theory. 



















































  Semiclassical instanton theory for reaction rates at any temperature: How a rigorous real-time derivation solves the crossover temperature problem
    Joseph E. Lawrence
    September 9, 2024
=====================================================================================================================================================















§ INTRODUCTION

The cornerstone of chemical reaction rate theory is Eyring transition state theory (TST).<cit.> Developed in the 1930s,<cit.> TST is still widely used today to estimate the rate of chemical reactions. The advantage of the theory is its simplicity, allowing for the estimation of the rate constant from only local knowledge of the Born-Oppenheimer potential in the reactant minimum and at the saddle point separating the reactants and products. The TST approximation to the rate has the simple form

    k_ TST  =  κ/2πβħZ^/Z_r  e^-β V^

where β=1/k_ B T is the inverse temperature, V^ is the height of the potential barrier, κ is the transmission coefficient, Z_r is the reactant partition function and Z^ is the transition state partition function for the modes orthogonal to the unstable coordinate. Typically, both Z_r and Z^ are treated quantum mechanically, using a harmonic approximation for the vibrations and rigid-rotors for the rotations. In contrast, in the basic version of the theory, motion over the barrier is assumed to be classical, which in the absence of dynamical recrossing leads to a transmission coefficient of κ=1. 


Already in the 1930s, it was clear to the developers of TST that quantum tunneling along the reaction coordinate could result in a transmission coefficient, κ>1.
When the effect of tunneling is small, a reasonable approximation is given by Wigner's famous tunneling correction<cit.>

    κ_ W  = 1 + (βħω)^2/24,

where ω is the barrier frequency along the unstable coordinate. This is, in fact, just the first term in an expansion of the exact tunneling correction for the parabolic barrier<cit.>     


    κ_ PB = βħω /2  /sin(βħω/2).

Nevertheless, κ_ W is typically preferred for practical purposes over κ_ PB as the parabolic barrier approximation diverges at the  “crossover temperature” T_ c=ħω/2π k_ B.

At low temperatures, where tunneling can enhance the rate by many orders of magnitude, such simple corrections are insufficient. Therefore, a commonly used approach is to retain the separable approximation of TST and compute a one-dimensional tunneling correction, e.g. by fitting the barrier to a function for which the result is known analytically,<cit.> or by using a semiclassical tunneling probability based on Wentzel–Kramers–Brillouin (WKB) theory.<cit.> These one-dimensional approximations, however, fail to capture non-separable effects such as “corner cutting”, where the system tunnels through higher but narrower regions of the potential.<cit.> 

Instanton theory provides a rigorous way to go beyond these simple one-dimensional tunneling corrections. The instanton is the dominant tunneling path in a full-dimensional path-integral description of the reaction.<cit.>

Hence, it is not constrained to follow the minimum energy path and can, therefore, capture key nonseparable 

effects on the rate.
In fact, rather than being a minimum energy path, the instanton is a stationary action path. Following Lagrange's principle, the instanton can thus be interpreted as a classical trajectory. Specifically, it is a periodic imaginary-time trajectory (which is equivalent to a real-time trajectory on the upturned potential) whose period is the thermal time, τ=βħ. The instanton can be found practically by optimisation of the action for a discretised path,<cit.>

which is not significantly more computationally challenging than locating the transition state.<cit.> 
This approach has been used to apply instanton theory to a wide range of systems in full dimensions, both with pre-computed potential energy surfaces and also on-the-fly using high-level electronic structure theory.<cit.>



Having located the instanton path, the instanton approximation to the rate constant can be written in a similar form to Eyring-TST<cit.>

    k_ inst =  1/√(2πħ)(- d^2 S_ inst/dτ^2)^1/2Z_ inst/Z_r e^-S_ inst(τ)/ħ




where S_ inst(τ) is the instanton action and Z_ inst the instanton partition function. Here, Z_ inst generalises Z^ by capturing the effect of the changing vibrational frequencies and rotational constants along the instanton path. This can be seen explicitly in the formal expression for Z_ inst, which in the absence of rotations is<cit.>

    Z_ inst(τ) = ∏_j=1^f-11/2sinh[u_j(τ)/2],

here u_j(τ) is the j^  th stability parameter for the instanton orbit, which for a separable system reduces to u_j(τ)=τω_j^ vib.
















Part of the power of instanton theory is that it is a rigorous semiclassical theory.<cit.> 


In particular instanton theory is the first term in an asymptotic series expansion of the exact quantum rate as ħ→0. 

Asymptotic expansions can be thought of as generalising the idea of a perturbative expansion to problems where a simple power series may not be applicable.<cit.> An important example is obtaining the expansion in terms of ε of an integral of the form

    I(ε) = ∫_-∞^∞ A(x)   e^-f(x)/ε dx

for ε→0. The asymptotic expansion for this integral can be found by Laplace's method (or equivalently steepest descent integration), resulting in an asymptotic series of the form<cit.>^,[Where it is assumed A(x^⋆)≠0.]  

    I(ε)  ∼ A(x^⋆)√(2πε/f”(x^⋆))  e^-f(x^⋆)/ε ( 1+ a_1 ε + a_2 ε^2 + … ),

where x^⋆ is the global minimum of f(x), i.e. f'(x^⋆)=0. We see that, at leading order, this corresponds to approximating the integrand by expanding f(x) to second order about x^⋆ and treating A(x) as constant. 
The apparent simplicity of this procedure belies the power and rich complexity of asymptotic analysis. 
For example, although asymptotic series are generally not convergent, their first few terms typically give a very good approximation to the exact result, and when both are available are usually much more numerically efficient than a Taylor series representation.<cit.> Furthermore, despite the apparently approximate nature of asymptotics, in principle a thorough asymptotic analysis in combination with modern resummation methods allows exact results to be recovered.<cit.>  


One might at this stage reasonably ask, what one means by an expansion in ħ? In particular, there is an apparent ambiguity in whether one defines the inverse temperature in terms of the thermal time, β=τ/ħ, or vice-versa τ=βħ. 
The ambiguity is resolved by a more formal definition, such that, ħ→0 is actually shorthand for introducing a perturbation parameter ε next to the ħ in the path-integral exponent

    ⟨q_ f| e^-Ĥτ/ħ|q_ i⟩ = ∫^q(τ)=q_ f_q(0)=q_ i𝒟q  e^-S_τ[q(τ')]/(εħ)

and considering the behaviour as ε→0.
 This prescription is equivalent to a WKB analysis in one dimension. By analogy with Eq. (<ref>) one sees that the resulting asymptotic expansion of the path integral will be around a classical trajectory (where the action is stationary).   






There exist two qualitatively different approaches to the derivation of instanton theory: one based on reactive scattering theory,<cit.> and the other on the concept of imaginary free-energy (the “F premise”).<cit.> 
The F premise was first proposed by Langer in the context of droplet formation,<cit.> and then later by Coleman in the context of quantum field theory.<cit.>
Proposed heuristically, the F principle says that the rate of decay of a metastable state is related to the imaginary part of its free-energy as k= -2/ħ F. This is then evaluated using analytic continuation to define the integral over the unstable mode in an asymptotic (ħ→0) evaluation of the path-integral expression for the partition function. Separately, Miller arrived at instanton theory by starting from the flux-correlation formulation for the thermal rate, using arguments based on Weyl correspondence in combination with semiclassical results from Gutzwiller.<cit.>
Although the prefactor in Miller's theory appears qualitatively different to that of the F formulation, the two theories have been proven to be equivalent.<cit.>^,[We note that it has been argued recently that there should be an additional term that does not appear in the standard expressions for the instanton rate considered in these works, and that this additional term arises due to non-separability.<cit.> However, we note that these arguments were not based on a rigorous semiclassical (ħ→0) analysis. Differences to the standard expressions may, therefore, simply be explained as arising from subdominant terms.] 

More recently Richardson has derived instanton theory from first principles by evaluating the exact flux-correlation expression for the quantum rate asymptotically as ħ→0.<cit.> Unlike the F derivation, which is rather subtle and thus difficult to generalise or extend, the flux-correlation formalism provides a rigorous framework for extending the theory, as has been exploited in recent years in the study of electronically nonadiabatic chemical reactions.<cit.> For this reason this is the approach taken in the present study.

Despite its success in describing deep tunneling processes, instanton theory has a major issue: it breaks down at high temperature. 
This can  be understood by noting that the shortest period for an instanton orbit is determined by the barrier frequency, τ_ c=2π/ω.

Hence, Eq. (<ref>) cannot be applied for temperatures above the same “crossover temperature”, T_ c=ħ/(k_ Bτ_ c), that appears when considering the parabolic barrier rate. 
This is not a coincidence, as we will see later, it is because the parabolic barrier rate is closely related to the instanton result and is the rigorous semiclassical result for T>T_ c.
As one approaches the crossover temperature from below, instanton theory becomes less accurate (although does not diverge).

This is clearly undesirable, as we would like to be able to accurately describe the onset of tunneling in chemical reactions.
 The goal of the present paper is therefore to derive a rigorous uniform semiclassical theory that is valid for all values of the thermal time τ=βħ>0. Here, uniform is a technical term that refers to an expression that is valid for a range of values of an additional (non-asymptotic) parameter that spans two (or more) regions exhibiting qualitatively different asymptotic behaviour. 
 

That instanton theory breaks down at the crossover temperature was, of course, known by its initial proponents.<cit.> As such, there have been many previous attempts to both analyse and ameliorate this problem.<cit.> In particular, 
in 1981 Affleck argued based on WKB that the F principle should satisfy: k=-2/ħ F for T≤ T_c and k=-ωτ/πħ F for T≥ T_c.<cit.>

These ideas have been developed further by several different authors<cit.> resulting in a piecewise theory.<cit.>
However, there are a number of drawbacks to this approach.
Being based on a heuristic combination of the F principle and WKB it is difficult to see how to rigorously extend the theory, for example by incorporating higher order asymptotic (i.e. perturbative)  corrections.<cit.>
Furthermore, the derivation is based on the fundamental assumption that the instanton smoothly collapses to the transition state with increasing temperature, and hence cannot be generalised to more complex systems.













An alternative approach that has been suggested is to start from the uniform WKB approximation to the 1D microcanonical transmission probability,<cit.> and then calculate the thermal rate by computing the integral over energy numerically.<cit.> However, these approaches suffer from similar issues to those based on F, and require ad-hoc approximations to treat non-separable multidimensional systems.


Recently Upadhyayula and Pollak have proposed a theory they call “uniform semiclassical instanton theory” that replaces this numerical integration with an analytical approximation. However, despite the method's name (as is shown in the supplementary material) the resulting theory is not a rigorous uniform semiclassical approximation in the sense described above, and actually becomes less accurate near τ=τ_c as ħ→0.

















































































































    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    































§ THEORY





















In order to derive our general semiclassical rate theory and solve the crossover problem we return to a first-principles derivation of instanton theory. The basis for our approach is the derivation of Richardson from Ref. , but with a small change to the flux operators. After a preliminary recap of the quantum flux-flux formulation of rate theory, Sec. <ref> introduces key definitions and derives the asymptotic form of the flux-flux correlation function. Section <ref> then discusses how, below the crossover temperature, the change made to the flux operators results in a new real-time perspective on instanton theory. Using this new perspective, Sec. <ref> provides a unified understanding of the breakdown of both instanton theory and the parabolic barrier approximation at the crossover temperature. Section <ref> then discusses how this unified understanding can be combined with modern methods for deriving uniform asymptotics expansions. 
Completing the derivation, Sec. <ref> introduces the key result of the paper and discusses its behaviour.












The starting point for our derivation is the exact expression for the rate in terms of the time integral of a flux-flux correlation function<cit.>
 
    k Z_r = ∫_0^∞  c_ ff(t) dt,




where we make use of the most general form of the correlation function, in which the two flux operators can be chosen to be different<cit.>
 
    c_ ff(t) =  [e^-(τ/2 - it)Ĥ/ħF̂_p e^-(τ/2 + it)Ĥ/ħF̂_r].


The flux operators are formally defined as Heisenberg time derivatives of projection operators onto reactants (P̂_r) and products (P̂_p) as


    F̂_r    = -i/ħ[Ĥ,P̂_r] 
    F̂_p    = i/ħ[Ĥ,P̂_p]


where the projection operators are given by


    P̂_r    = 1-h(σ(q̂)-s_r) 
    P̂_p    = h(σ(q̂)-s_p)


here h(x) is the Heaviside step function, such that σ(q)<s_r defines the reactants and σ(q)>s_p the products.
Thus F̂_r measures the flux out of the reactant states (through the “dividing surface” σ(q)=s_r) and F̂_p measures the flux into the product states (through the “dividing surface” σ(q)=s_p).
 



















In Richardson's derivation<cit.> he took both flux operators to be the same and defined the dividing surface to pass through the saddle point of the potential.

Such a choice ensures that the flux correlation function has its maximum at t=0.
This is a natural choice when considering instanton theory as an extension of transition state theory, as integrating over time asymptotically (as ħ→0) then straightforwardly leads to a theory that involves no real-time quantities.


However, in the present work we choose instead to place the dividing surfaces well away from the barrier (out in the reactant and product asymptotes respectively) such that the correlation function reaches its maximum for finite real time. 

The reader might be worried that the presence of real time will mean that we must contend with the infamous sign problem.
This concern would be realised if we tried to evaluate the resulting path-integral expressions exactly or via a quantum instanton approximation<cit.> (an idea discussed in Ref. ).
However, as we will evaluate all integrals analytically using asymptotic analysis our semiclassical theory will have no such issue.
In fact, instead of making the problem more difficult, the presence of real time will actually be the key to solving the problem at the crossover temperature.











 §.§ Asymptotic expression for the flux-flux correlation function

Before we evaluate the integral over time, we begin by evaluating the flux-flux correlation function asymptotically. The steps we follow in this subsection follow closely those of Ref. , and are included here for completeness. However, for the sake of notational simplicity we restrict the discussion here to a system in one dimension.


To further simplify notation we define the forward and backward times, t_± =-iτ_± =± t-iτ/2, and introduce the imaginary time propagator K̂(θ)=e^-θĤ/ħ such that 

    c_ ff(t) = [K̂(τ_-)F̂_p K̂(τ_+) F̂_r].

Note that τ_± are therefore complex numbers, the imaginary part of which is (±) the “real time”.

In one dimension the flux operators can be written in the form

    F̂_α = p̂/2mδ(q̂-s_α)+δ(q̂-s_α)p̂/2m

for α=r or p, where δ(x) is the Dirac delta function.
By first inserting resolutions of the identity 1=∫dq |q⟩⟨q| and then making use of the relation ⟨q|p̂=-iħ∂/∂ q⟨q| followed by integration over the resulting Dirac delta functions, it is straightforward to show that 

    c_ ff(t) =  (-iħ/2m)^2  (    ∂ K(τ_-,s_p,s_r)/∂ s_r ∂ K(τ_+,s_r,s_p)/∂ s_p
    
         -    ∂^2 K(τ_-,s_p,s_r)/∂ s_r ∂ s_p K(τ_+,s_r,s_p) 
    
         -    K(τ_-,s_p,s_r)∂^2 K(τ_+,s_r,s_p)/∂ s_r ∂ s_p
    
         +   ∂ K(τ_-,s_p,s_r)/∂ s_p ∂ K(τ_+,s_r,s_p)/∂ s_r).

where the position representation of the propagator is defined as

    ⟨q”|K̂(θ)|q'⟩=K(θ,q',q”) .


























Now we turn to the asymptotic evaluation of this exact expression as ħ→0. This can be done numerically by first writing the propagator as a discretised path integral and then evaluating the integrals over position by steepest descent [i.e. using Eq. (<ref>)]. For the present purpose, we note that in the continuum limit this is formally equivalent to using the semiclassical (imaginary-time) van-Vleck propagator       

    ⟨q”|K̂(θ)|q'⟩=K(θ,q',q”) ∼∑_ paths(C/2πħ)^1/2 e^-S(θ,q',q”)/ħ

where S(θ,q',q”) is the Euclidean action calculated along the classical trajectory from q' to q” in total imaginary time θ, 
    S(θ,q',q”) = ∫_0^θ1/2m q̇^2(θ')+V(q(θ'))  dθ',

the prefactor is defined as

    C =  -∂^2 S/∂ q' ∂ q”,

and the sum is over all classical paths that go from q' to q” in time θ. Note the branch of the square root in the prefactor is chosen so that the prefactor is continuous along the trajectory.[The prefactor is of course a property of the entire trajectory. The prefactor being continuous along the trajectory should be understood as corresponding to considering the set of prefactors that are given by fixing q' and then varying θ and q” such that they follow the trajectory, and ensuring that this set is continuous for some parameterisation of the trajectory.]

In addition to the asymptotic expression for the propagator we will also need to make use of the asymptotic expressions for its derivatives


    ∂K/∂ q'   ∼∑_ paths -1/ħ∂ S/∂ q'(C/2πħ)^1/2 e^-S(θ,q',q”)/ħ
    ∂K/∂ q”   ∼∑_ paths-1/ħ∂ S/∂ q”(C/2πħ)^1/2 e^-S(θ,q',q”)/ħ
    ∂^2K/∂ q' ∂ q”   ∼∑_ paths1/ħ^2∂ S/∂ q'∂ S/∂ q”(C/2πħ)^1/2 e^-S(θ,q',q”)/ħ.





































Combining Eq. (<ref>) with Eqs. (<ref>)-(<ref>) and retaining just the dominant path we then obtain the following asymptotic expression for the correlation function (valid as ħ→0)

    c_ ff(t) ∼   (∂ S_ -/∂ s_r∂ S_ +/∂ s_p-∂ S_ -/∂ s_p∂ S_ -/∂ s_r
        -∂ S_ +/∂ s_p∂ S_ +/∂ s_r+∂ S_ -/∂ s_p∂ S_ +/∂ s_r)
       ×-1/4m^2√(C_+C_-/(2πħ)^2) e^-[S_ +(τ_+,s_r,s_p)+S_ -(τ_-,s_p,s_r)]/ħ,

where we have labelled the action (and associated quantities) ± to distinguish between the forward and backward paths. Hence, defining the total Euclidean action as the sum of the forward and backward parts

    S(t) = S_ +(τ/2+it,s_r,s_p) + S_ -(τ/2-it,s_p,s_r)

and the ħ independent part of the prefactor as

    A(t) =(C_+C_-)^1/2/8m^2π(∂ S_ -/∂ s_p∂ S_ -/∂ s_r
        -∂ S_ -/∂ s_r∂ S_ +/∂ s_p-∂ S_ -/∂ s_p∂ S_ +/∂ s_r+∂ S_ +/∂ s_p∂ S_ +/∂ s_r)

we see that the correlation function has the simple asymptotic form

    c_ ff(t) ∼A(t)/ħ e^-S(t)/ħ.
 
Note that, while the preceding analysis was specific to a one-dimensional system, 

the flux-flux correlation function for multidimensional systems also has the same form as Eq. (<ref>).



   














 §.§ A real-time derivation of instanton theory





Combining Eq. (<ref>) with Eq. (<ref>) we can recover the standard instanton theory below the crossover temperature by integrating over time asymptotically using Eq. (<ref>) to obtain

    k Z_r ∼A(t^⋆)/ħ√(2πħ/S”(t^⋆))e^-S(t^⋆)/ħ as ħ→ 0,




where t^⋆ satisfies the steepest descent condition S'(t^⋆)=0. We will now argue that, S(t^⋆)=S_ inst(τ), and that Eq. (<ref>) is equivalent to the standard instanton theory. 

First note that S'(t^⋆)=0 implies 

    ∂ S_ +/∂τ_+ - ∂ S_ -/∂τ_- = E_+ - E_- =0

i.e. the energy of the forward and backward paths, E_±, must be the same. Further, since K(θ,q',q”)=K^*(θ^*,q”,q') it follows that, on the real t axis, the action for the forward and backward paths are always complex conjugates S_ +(τ_+,s_r,s_p)=S^*_-(τ_-,s_p,s_r) making S(t) real. Differentiating this with respect to τ then gives E_+=E_-^* which combined with Eq. (<ref>) shows that the energy of the forward and backward trajectories at the stationary time are also real.

Before we discuss finding a trajectory that satisfies these conditions, note that there is a freedom we have not yet discussed: the contour of integration (the time path) in the definition of the action [Eq. (<ref>)]. This freedom corresponds to the order of the real- and imaginary-time propagators in the path-integral discretisation of Eq. (<ref>). Of course the exact expression is independent of this choice because real- and imaginary-time propagators commute. The semiclassical propagator must necessarily retain this property. This can be seen explicitly in Eq. (<ref>) as being a result of Cauchy's integral theorem. We are therefore free to choose the most convenient time path for our purposes (under the assumptions of the theorem). 




With a careful choice of time path it is 
 trivial to find a trajectory that satisfies S(t^⋆)=S_ inst(τ). Figure <ref> depicts the corresponding trajectory and time path. 
 The forward and backward time paths each consist of three segments; two pure real-time segments either side of a pure imaginary-time segment of length τ/2. 

Rather than defining the lengths of the real time segments (t_r and t_p) we instead define the trajectory along the imaginary-time segment to be half of the instanton orbit from one turning point to the other.

This then uniquely determines the real-time sections of the trajectories along with the times t_r and t_p. To see this note that, because the momentum is continuous along the trajectory it must be zero at the points connecting the real-time and imaginary-time segments. Hence, t_r (t_p) must be the time it takes for the system to roll from the reactant (product) end of the instanton to the reactant (product) dividing surface.  This, therefore, uniquely determines the stationary time as t^⋆=t_r+t_p. Since the forward and backward trajectories are entirely real and follow the same path it is clear that E_+=E_- and hence [cf. Eq. (<ref>)] the trajectory corresponds to a stationary time, t^⋆. 
 





















































The stationary trajectory gives an intuitive picture of 

the reaction process. 
The system starts at the reactants and moves in real time towards the barrier. Upon reaching the turning point of the trajectory at the barrier, instead of bouncing off the barrier in real time it switches to imaginary time. This effectively “turns the barrier upside-down” allowing the system to tunnel through to the product side. The system then switches back to real time and carries on to the products. Because of the cyclic nature of the trace in Eq. (<ref>) the system then retraces its steps, moving backwards in real time to the barrier before tunneling through the barrier in (positive) imaginary time and then back to the reactants again in negative real time. As the forward and backward real-time segments follow the same paths, their contributions to the action exactly cancel one another leaving only the imaginary time contribution, and hence S(t^⋆)=S_ inst(τ). 



Importantly, the uniqueness of asymptotic series means that the prefactor must also be equivalent to the usual instanton prefactor.  
Hence, we have that 

    A(t^⋆)/ħ√(2πħ/S”(t^⋆)) = Z_ inst(τ)/√(2πħ)(-d^2S_inst/dτ^2)^1/2.








For the inquisitive reader, a short aside: There are two obvious questions based on the preceding discussion. First, what would happen to the trajectory if we kept t^⋆ (and τ) fixed but deformed the time path? The answer is that the resulting trajectory must move into the complex position plane. This can be understood by noting that, if the momentum of the system is real then propagation in anything other than real time will lead to a change in the position that is complex. As the initial and final momenta are fixed changing the direction of the time path at any point (other than turning points) will, therefore, result in a complex trajectory. The second question is, when t is not t^⋆ can a time path still be found that keeps the trajectory real? The answer is that while this is possible for some t and τ in one dimension, it is generally not possible in multiple dimensions. Finally, the observant reader may note that there is another stationary trajectory that swaps the sign of the real time segments on the product side of the barrier to give t^⋆=t_r-t_p. One might, therefore, wonder why this trajectory doesn't also contribute. A careful consideration, however, shows that this corresponds to a different branch of the correlation function and hence does not need to be included.



















































































































 §.§ Diagnosing the problem at the crossover temperature
 
With this new real time perspective we can now obtain a simple intuitive picture of what happens as we approach the crossover temperature, and hence why the standard instanton theory stops working there. First consider approaching the crossover temperature from below.
Figure <ref> shows how the steepest descent time, t^⋆, varies with τ for a typical system close to the crossover. We see that as  τ→τ_ c the stationary time approaches infinity, t^⋆→∞. This can be understood intuitively by noting that as τ→τ_ c the turning points, where the real-time and imaginary-time parts of the path meet, get closer and closer to the top of the barrier. As this happens, the force at the turning points approaches zero. Hence, the  real-time segments have to become longer and longer to give time for the system roll away from (or come to a stop at) the barrier.






Now we consider the behaviour above the crossover temperature. Based on our preceding discussion it is clear that in this regime the correlation function must be dominated by the behaviour near t=∞. For a one-dimensional barrier we show in the supplementary material (following Ref. ) that for large values of real time the action behaves like

    S(t) ∼ S_∞ + a sin(ωτ/2)  e^-ω t as  t→∞

where S_∞=τ V^ and a is a constant with units of action. Furthermore, the pre-exponential factor varies as 

    A(t) ∼aω^2/4π e^-ω t as  t→∞.

Hence, the correlation function obeys

    c_ ff(t) ∼aω^2/4πħ e^-ω texp(-[S_∞ + a sin(ωτ/2)  e^-ω t]/ħ)

as ħ→0 and t→∞.











Upon making the substitution v=e^-ω t we can see that (above the crossover temperature) the rate is asymptotic to[Note that after making the substitution we also extend the integration range from 0 to 1 to be from 0 to ∞. This is the correct thing to do because as ħ→0 the integrand becomes more and more narrowly peaked at the boundary v=0.]

    k Z_r    ∼∫_0^∞aω/4πħexp(-[S_∞ + a sin(ωτ/2)  v]/ħ) dv 
       = ω/4πsin(ωτ/2)e^-S_∞/ħ≡ω/4πsin(βħω/2) e^-β V^

which is the well-known parabolic barrier rate.<cit.> Note that this expression differs qualitatively in its dependence on (the asymptotic parameter associated with) ħ from the instanton result evaluated below the crossover temperature, Eq. (<ref>). 


We can now gain a unified perspective on the behaviour both above and below crossover. To do so we begin by making the same variable transformation as above, v=e^-ω t, such that the rate can be expressed as

    k Z_r ∼∫_0^1 G(v)/ħ e^-S(v)/ħd v   as ħ→0

where we define G(v):=A(t(v))/(vω) and S(v):=S(t(v)).[Importantly, with this definition G(v) is a smooth function of v for all v≥0.] It is important to recognise that we will still recover the standard instanton result if we perform the integral over v asymptotically below the crossover temperature [i.e. using Eq. (<ref>)]. 

This can be seen explicitly by first integrating asymptotically over v to give

    k Z_r ∼G(v^⋆)/ħ√(2πħ/S”(v^⋆))e^-S(v^⋆)/ħ as ħ→ 0

and then using the chain rule  -vω S'(v)=S'(t) to show that S(t^⋆)=S(v^⋆) and (v^⋆ω)^2 S”(v^⋆)=S”(t^⋆).
The advantage of making this transformation is that the behaviour at the crossover temperature becomes very simple. It just corresponds to the temperature at which the stationary point, v^⋆, moves from being inside to outside of the integration range.   






Figure <ref> depicts e^-S(v)/ħ for three temperatures, one above crossover, one below crossover, and the other exactly at crossover. We see that, well below the crossover temperature, e^-S(v)/ħ is peaked far away from the boundary, and hence approximating the integrand as a Gaussian is a reasonable approximation. Equally, well above the crossover temperature we see that only the tail of the function e^-S(v)/ħ appears inside the integration region, and hence approximating the integrand with an exponential [as is done in Eq. (<ref>)] is also a reasonable approximation. However, exactly at the crossover temperature the stationary point S'(v^⋆)=0 occurs exactly on the integration boundary, v^⋆=0. 

This means that essentially half of the peak is outside of the integration range.
Hence, instanton theory (which assumes the full peak in inside the integration bounds) is approximately a factor of two too large. When approaching the crossover temperature from above things are much worse, as the first derivative at the boundary approaches zero and hence an exponential approximation will give a divergent integral.       




















 §.§ Uniform asymptotics
 
To obtain a rate theory that is valid both above and below the crossover temperature we need to make use of ideas from uniform asymptotics. A uniform asymptotic expression is one that is continuous in a parameter (distinct from the asymptotic parameter) that connects two or more regions of different asymptotic behaviour. We have seen in the previous section that, as τ varies the stationary point moves from inside to outside of the integration range, resulting in a qualitative change in the form of the asymptotic approximation for the rate. We are, therefore, interested in finding an asymptotic expression for the rate that is uniform in τ.

















A common approach is to construct uniform asymptotic approximations in a piecewise manner. Following this approach one might be tempted to suggest that, below the crossover temperature,  the integral over v should be approximated by the Gaussian integral 


    k Z_r ∼G(v^⋆)/ħ∫_0^∞ e^-[S(v^⋆)+(v-v^⋆)^2S”(v^⋆)/2]/ħdv,
 
and that this should then be combined with an expression valid above crossover that has the same value as τ→τ_ c.


This kind of uniform approximation is employed for example in Ref. .

The resulting expression has lots of the behaviour that one expects, it reduces the prediction by a factor of 1/2 at the crossover temperature and smoothly approaches the standard result as one lowers the temperature for fixed ħ (and also as one lowers ħ for fixed τ). 
However, this is not the correct approach. 

To see what is wrong with Eq. (<ref>), we can evaluate the integral and make use of the chain rules given earlier to show that it is equivalent to 

    kZ_r ∼A(t^⋆)/ħ√(2πħ/S”(t^⋆) ) e^-S(t^⋆)/
        ħerfc(- √(S”(t^⋆) /2ω ^2ħ))/2

where erfc(x) is the complementary error function. Crucially, we see that the difference between this and the standard instanton result is the factor,

    1/2erfc(- √(S”(t^⋆) /2ω ^2ħ)).

While this factor takes on values between 1/2 and 1, the precise value is dependent on the ratio S”(t^⋆) /2ω ^2ħ. This is a problem because the actual value of S”(t^⋆) is determined by the location of the dividing surfaces. This is clearly unphysical as the true quantum rate is independent of the choice of dividing surface. [As mentioned earlier, the standard instanton result is also independent of dividing surface.] Clearly we want our uniform asymptotic result to also have this property.

Fortunately, if one consults a textbook on uniform asymptotics<cit.> one will see that Eq. (<ref>) is not the recommended approach. 

Instead, 

one should use Bleistein's method<cit.> which has become the standard method for generating uniform asymptotic series.<cit.> Our final result will then be the first term in this series, just as the standard instanton result is 

the first term in a regular asymptotic series. This has the advantage that the theory will not only be continuous in τ, but it will also be smooth and continuously differentiable.


Following Bleistein's method we begin by defining a new variable transformation

    S(v) = u^2/2 - bu + c

and choosing the constants b and c so that u(v=0)=0 and u(v^⋆)=b. This then implies that


    b =  sgn(v^⋆)√(2S_∞-2S(v^⋆))


    c = S(v=0) = S_∞.


Note that we assume here that S(v) has been analytically continued to v<0 such that for τ<τ_ c we can find v^⋆<0.
Solving for u we then obtain

    u = b +sgn(v-v^⋆) √(2S(v)-2S(v^⋆))

where the sgn(v-v^⋆) ensures that the variable transform is single valued. After this variable transformation we have

    k Z_r ∼ e^-c/ħ∫_0^∞g(u)/ħ e^-(u^2/2-bu)/ħd u   as  ħ→0

where g(u) = G(v(u)) dv/du.



The final step in the Bleistein method is to write this pre-exponential term as a linear function that passes through the function at the boundary and the stationary point (u=b) plus a remainder term 




    g(u) = g(b) + (u-b)g(b)-g(0)/b + r(u).

As the remainder r(u) is zero at both the stationary point and the boundary it can be ignored at leading order in the uniform asymptotic expansion. With this we can then perform all integrals analytically to obtain








    kZ_r    ∼ e^-c/ħg(b)/ħ√(2πħ) e^|b^2/2ħ1/2erfc(-b/√(2ħ))
       +e^-c/ħg(b)-g(0)/b



To analyse this result we need to express it explicitly in terms of S(t) and A(t).
We begin by defining S_∞-S(t^⋆)=Δ S and noting that sgn(v^⋆) can be rewritten in terms of τ to give

    b =sgn(τ-τ_ c)√(2Δ S).

Second, we note that b^2/2 can be combined with c to give the instanton action

    -b^2/2 + c = S(t^⋆) = S_ inst(τ).

To complete our simplifications we need to determine explicit expressions for g(b) and g(0). This can be achieved by observing that 

    du/dv = sgn(v-v^⋆)S'(v)/√(2 S(v)-2S(v^⋆))

which can then be combined with the definition of g(u) to give

    g(u)    = G(v)dv/du = -A(t)dt/dvdv/du
       =sgn(v^⋆-v)A(t) √(2S(t)-2S(t^⋆))/S'(t).







Hence, evaluating this at u=0 and taking the limit as u→ b from above or below gives


    g(0)   = sgn(τ-τ_ c) √(2Δ S)lim_t→∞A(t)/S'(t)


and

    g(b) = A(t^⋆)/√(S”(t^⋆)).


Combining these results we obtain 


    k Z_r    ∼   e^-S_∞/ħ(   sgn(τ-τ_ c) A(t^⋆)/√(2 S”(t^⋆)Δ S )-lim_t→∞A(t)/S'(t))
       +e^-S(t^⋆)/ħA(t^⋆)/ħ√(2πħ/S”(t^⋆)) 1/2erfc(sgn(τ_ c-τ)√(Δ S/ħ))

which is a rigorous uniform asymptotic expression for the thermal rate constant that correctly bridges between the parabolic barrier result above the crossover temperature and the instanton result below crossover.

Although we motivated the derivation in terms of a one-dimensional system, this expression is valid for any system in which the instanton collapses smoothly to the transition state. Furthermore, while we have used the real-time formulation to derive the theory, it can immediately be rewritten in terms of quantities that involve only imaginary time. Hence, it is clearly independent of the choice of dividing surface.

Making use of the relations





    A(t^⋆)/√(S”(t^⋆)) = Z_ inst(τ)/2π(-d^2S_inst/dτ^2)^1/2

and 

    -lim_t→∞A(t)/S'(t) =  ω/4πsin(τω/2)Z^(τ)


allows us to write








    k Z_r    ∼   e^-τ V^/ħsgn(τ-τ_ c)Z_ inst(τ)/2π√(2Δ S)(-d^2S_inst/dτ^2)^1/2
       +e^-τ V^/ħω/4πsin(τω/2)Z^(τ) 
       +e^-S_inst(τ)/ħZ_ inst(τ)/√(2πħ)(-d^2S_inst/dτ^2)^1/21/2erfc(sgn(τ_ c-τ)√(Δ S/ħ))

which is equivalent to the more compact expression


    k     ∼  k_m-pb    +sgn(τ-τ_ c)k_ inste^-Δ S/ħ/√(4πΔ S/ħ)
       +k_ inst 1/2erfc(sgn(τ_ c-τ)√(Δ S/ħ))

where k_m-pb is the multidimensional parabolic barrier rate. 


We are very nearly done with our derivation, however, the keen-eyed reader will note that Eqs. (<ref>), (<ref>), and (<ref>) are still not valid at all temperatures. This can be seen by noting that (sin(τω/2))^-1 in the parabolic barrier rate, k_m-pb, not only diverges at τ_ c=2π/ω, but also at all τ_ c,n=2π n/ω for n=1,2,…,∞, i.e. half the crossover temperature, a third of the crossover temperature, and so on. Again the real-time formulation makes it easy to understand the cause of these divergences. Each one corresponds to a different stationary point of the action passing through the v=0 boundary.

We can give a physical interpretation to these stationary points by observing that τ_c,n corresponds to the shortest possible time for a trajectory which completes n orbits on the upturned potential. Hence, we can attribute these stationary points to the trajectories that involve multiple orbits. 
For example, when τ_c,3>τ>τ_c,2 (6π/ω>τ>4π/ω) there exists not one but two periodic trajectories on the upturned potential with a period τ. One is the usual instanton, that orbits just once, and the other is a trajectory that orbits twice in the same time, with each half of the trajectory following the same path as the standard instanton at twice the temperature. 

These multi-orbit instantons (typically referred to as periodic instantons) have a higher action and hence are exponentially suppressed compared to the one-orbit instanton. Within standard (Poincaré) asymptotics these terms are, therefore, not included as they are smaller than every term in the one-orbit instanton's asymptotic series (as ħ→0). However, we now see that in order to obtain a rigorous uniform theory valid at any temperature we are naturally led to include them. 





 §.§ Semiclassical instanton theory valid at any temperature








Given the preceding discussion it is clear that in order to cancel the divergences in k_m-pb we must modify Eq. (<ref>) by including multi-orbit instanton terms. The form of these terms can easily be determined by inspection. However, we do not have to rely on inspection alone. 

As shown in the Appendix, in one dimension we can rigorously derive the desired uniform thermal rate theory in an entirely different way. There, we start from the uniform energy dependent WKB transmission probability and then obtain the thermal rate by integrating over energy asymptotically using Bleistein's method.<cit.> 
The resulting theory contains exactly the multi-orbit terms we were expecting. Comparing this one-dimensional result  [Eq. (<ref>)] with  Eq. (<ref>) it is trivial to generalise to the multidimensional multi-orbit case. 
We thus arrive at the central result of the present paper, a uniform asymptotic expression for the thermal rate in multiple dimensions




    k     ∼   k_m-pb+ ∑_n=1^∞ (-1)^n+1sgn(τ-nτ_ c) k_n, inste^-Δ S_n/ħ/√(4πΔ S_n/ħ)
       +∑_n=1^∞ (-1)^n+1 k_n, inst 1/2erfc(sgn(nτ_ c-τ)√(Δ S_ n/ħ))

where the n-orbit instanton action is defined as

    S_n,inst(τ) = n S_ inst(τ/n)

and the difference to the collapsed/classical action is given by

    Δ S_n = τ V^ - S_n,inst(τ).

The only remaining terms to define are the effective n-orbit instanton rate constants that are given by

    k_n,instZ_r = e^-S_n,inst(τ)/ħZ_n,inst(τ)/√(2πħ)(-d^2S_n,inst/dτ^2)^1/2








with Z_n,inst(τ) the instanton partition function for the n-orbit instanton. 
In the absence of rotations this is given by

    Z_n, inst(τ) = ∏_j=1^f-11/2sinh[n u_j(τ/n)/2],

where again u_j(τ) is the j^ th stability parameter for the 1-orbit instanton of period τ. 
Note that, just as with the standard instanton theory, Eq. (<ref>) is rigorously independent of dividing surface and satisfies detailed balance.
We stress again that, although the theory is applicable to a wide range of multidimensional systems, our analysis, and hence Eq. (<ref>), assumes that the instanton collapses smoothly to the transition state as the temperature is increased. Notable exceptions, such as quartic barriers, where there are multiple (interacting) instantons with the same τ are important exceptions<cit.> and will be the subject of future work.



  §.§.§ Understanding the theory

Before demonstrating the accuracy of the theory numerically, we begin 
with a  qualitative discussion the terms that appear and how they interact with one another.

Perhaps the first thing one notices, is that the n-orbit terms appear with alternating signs. In the  derivation from WKB presented in the Appendix, these alternating signs occur as a direct consequence of the form of the uniform transmission probability. However, we note that the necessity of the sign alternation is also evident from Eqs. (<ref>) and (<ref>)
as it is required to match the alternating divergences of (sin(ωτ/2))^-1 in k_m-pb. To understand this cancellation explicitly, we consider the behaviour of k_m-pb about τ→ nτ_ c. Using the following expansion around τ=nτ_ c

    ω/4πsin(τω/2)∼(-1)^n/2π(τ-nτ_ c) +𝒪(τ-nτ_ c)

one can show that the parabolic barrier rate behaves like

    e^τ V^/ħ k_m-pbZ_r∼Z^(nτ_ c)(-1)^n/2π (τ-nτ_ c) +Z^'(nτ_ c)(-1)^n/2π ,

again with an error of 𝒪(τ-nτ_ c).
To see how this divergent behaviour is cancelled, we expand the terms appearing in the first sum about τ=nτ_ c, retaining terms up to 𝒪(τ-nτ_ c), to give

    sgn(τ-nτ_ c)/√(4πΔ S_n/ħ)∼√(nħ/2π/- S”_ inst(τ_ c))(1/τ-nτ_ c-S”'_ inst(τ_ c)/6n S”_ inst(τ_ c)).




Using this result, and noting that Z_n,inst(nτ_ c)=Z^(nτ_c) it can then be shown that

    sgn(τ-nτ_ c) k_n, instZ_re^+S_n,inst/ħ/√(4πΔ S_n/ħ)   ∼Z^(nτ_ c)/2π (τ-nτ_ c)
       + Z^(nτ_c)/2πS”'_ inst(τ_ c)/3n S”_ inst(τ_ c)
       + Z_n,inst'(nτ_c)/2π

again to 𝒪(τ-nτ_ c).
Comparing Eq. (<ref>) and Eq. (<ref>) we see that [once we combine  Eq. (<ref>) with the factor of (-1)^n+1] the divergent terms exactly cancel leaving just a constant as τ→ nτ_c. From this we can see that, exactly at the crossover temperature, our theory predicts a correction to the rough “factor of 2 error” of instanton theory discussed earlier. Specifically, (ignoring the hyperasymptotic multi-orbit terms) we find that 

    k(τ_c)- 1/2k_ inst(τ_ c)  ∼1/2πe^-τ_ c V^/ħ(   Z^(τ_c)/Z_r(τ_c)S”'_ inst(τ_ c)/3S”_ inst(τ_ c)
    
    +   Z_ inst'(τ_c)/Z_r(τ_c)-Z^'(τ_c)/Z_r(τ_c)).


Having considered the behaviour of the theory close the crossover temperature(s) where Δ S_n=0 let us know consider the behaviour of the new terms for Δ S_n/ħ≫ 1, both when τ>nτ_ c and τ<nτ_ c. The simpler of these two cases is τ>nτ_ c, where (half) the complementary error function approaches one and hence we have

    k_n, inste^-Δ S_n/ħ/√(4πΔ S_n/ħ)  + k_n, inst 1/2erfc(-√(Δ S_n/ħ))  ∼  k_n, inst.

The slightly more complicated case is τ<nτ_ c. Here we can make use of the standard asymptotic result

    erfc(x)∼e^-x^2/√(π x^2)(1-1/2x^2+…) as x→∞

from which we observe that

    -k_n, inste^-Δ S_n/ħ/√(4πΔ S_n/ħ)  + k_n, inst 1/2erfc(√(Δ S_n/ħ))  ∼ -k_n, inste^-Δ S_n/ħ/4√(π(Δ S_n/ħ)^3)

which as we would expect is subdominant to k_m-pb even for n=1.  

 






















































The derivation from the WKB transmission probability given in the Appendix suggests a natural separation of the theory into three parts: classical above barrier transmission, quantum above barrier reflection, and quantum tunnelling. Making this separation we can write Eq. (<ref>) as

    k ∼ k_ TST - k_ reflect + k_ tunnel

where

    k_ TST = 1/2πτZ^(τ)/Z_r(τ)e^-τ V^/ħ

is the above barrier transmission contribution (equivalent to the TST rate with κ=1). The reflection rate can then be expressed as  

    k_ reflect = 1/2πZ^(τ)/Z_r(τ)e^-τ V^/ħ∑_λ=1^∞ (-1)^λ+11/λτ_ c+τ

where to avoid potential confusion in the following sections we have used λ rather than n as the dummy index in the sum.
Note that as defined k_ reflect is always positive, and can be expressed as k_ reflect=ϕ(τ/τ_ c)k_ TST where ϕ(x) is a system independent function, with ϕ(0)=0, ϕ(1)=1-ln(2), and ϕ(∞)=1/2. Finally the tunneling contribution, which is system dependent can be expressed as

    k_ tunnel   = 1/2πZ^(τ)/Z_r(τ)e^-τ V^/ħ∑_λ=1^∞ (-1)^λ+11/λτ_ c-τ
       + ∑_n=1^∞ (-1)^n+1sgn(τ-nτ_ c) k_n, inste^-Δ S_n/ħ/√(4πΔ S_n/ħ)
       +∑_n=1^∞ (-1)^n+1 k_n, inst 1/2erfc(sgn(nτ_ c-τ)√(Δ S_ n/ħ)).







  §.§.§ Numerical Considerations

Having given a qualitative interpretation of the theory we now turn to some practical considerations about the numerical implementation of the theory. First, the presence of infinite sums in Eq. (<ref>) might appear daunting, and raises the obvious question: How many terms are needed to reach numerical conversion? Clearly, one requires  enough terms to avoid the divergence of k_m-pb at the temperature of interest. We will see in the next section that this is the dominant consideration and that for 0<τ<2τ_ c excellent convergence is obtained with only n=2. It should be stressed that if implementing the theory using Eqs. (<ref>) and (<ref>) one must include a sufficient number of terms in the sum over λ to recover the parabolic barrier rate above the crossover temperature, and hence the maximum value of λ may be higher than n. Note that including an arbitrarily large number of terms in this sum is trivial. 





One remaining aspect that we have not yet addressed is the meaning of S_ inst(τ) [and Z_ inst(τ)] when τ<τ_ c.
Because there exist no real-position periodic orbits for these values of imaginary time, the resulting trajectories must move into the complex position plane. Finding such trajectories is clearly impractical for realistic chemical applications. However, since the terms containing S_ inst(τ<τ_ c) are always subdominant we can develop numerical approximations of the action in this region that require only real positions without affecting the key asymptotic behaviour of the theory.

In the following section we will give an example of how this can be for one-dimensional systems, and compare to the results obtained using the exact action.





§ NUMERICAL RESULTS




 §.§ Symmetric Eckart Barrier

To illustrate the accuracy of the new theory we consider the prototypical one-dimensional model of reactive scattering: the symmetric Eckart barrier. The potential for the symmetric Eckart barrier is defined as

    V(q) = V^^2(q/L).

For this simple model system the instanton action can be evaluated analytically as




    S_ inst(τ) = 2  τ_ cV^ - τ_ c^2 V^/τ

where the barrier frequency is given by





    ω = √(2 V^/m L^2).





To aid comparison with previous work<cit.> we consider the following parameters L=0.66047  a_0, m=1836 m_e, V^=72ħ^2/(m π^2  L^2). The exact result was calculated for comparison by numerical integration of 

    k Z_r = 1/2πħ∫_0^∞  P(E)e^-τ E/ħdE

using the exact analytical result for the transmission probability, P(E).<cit.>






Figure <ref> compares the rate calculated with the new theory against the exact, classical, parabolic barrier, and standard (old) instanton results for the symmetric Eckart barrier.  
The top panel shows the logarithm of the rate for each theory (multiplied by Z_rτ_ c to ensure it is dimensionless), and the bottom panel shows the ratio of the approximate theories to the exact result (plotted on a logarithmic scale so that overestimation and underestimation are treated on an equal footing).
One immediately notices the divergence of the parabolic barrier rate at the crossover temperature (β/β_ c=1) and the (approximate) factor of two error of the standard instanton result.


In contrast to both the standard theory and parabolic barrier rate, the new theory smoothly connects the low and high temperature limits. The accuracy of the new theory is as good as one could have hoped, with the error no greater in the vicinity of the crossover temperature than it is at low temperature. It is important to stress that quantum effects are not insignificant in this regime. The quantum rate more than a factor of 6 larger than the classical rate at T_ c and still more than factor 2 larger than the classical rate even at  T=1.5× T_ c. This highlights the importance of the new theory in accurately describing the onset of quantum tunneling.



At lower temperatures, we see that the error of both the new theory and the standard instanton have approximately plateaued at τ=2τ_ c (half the crossover temperature). The close agreement between the standard theory and the new theory at τ=2τ_ c indicates that in this regime the multi-orbit terms do not significantly affect the rate. However, it is important to remember that if one were to neglect these terms and try to use Eq. (<ref>) instead of Eq. (<ref>) the result would diverge at τ=2τ_ c. One could of course consider approximating the full theory [Eq. (<ref>)] by removing both the multi-orbit terms and making appropriate modifications to k_m-pb to remove the divergences. This would break the formal asymptotic properties of the theory, but may nevertheless provide a useful approximation to the rate. 






 §.§ Asymmetric Eckart Barrier



Next we consider the asymmetric version of the Eckart barrier, for which the potential can be written as




    V(q) = (√(V_1)+√(V_2))^2/4cosh^2(q/L) -V_2-V_1/1+exp(-2q/L),

note that V_1=V^.
For this system evaluating the action analytically is slightly more involved. Rather than give the expression for S_ inst(τ) explicitly we  instead note it can be computed from knowledge of the reduced action, W(E). Note the reduced action is the Legendre transform of the action, S_ inst(τ)=W(E)+τ E. 
The reduced action can be expressed as<cit.>





    W(E) = A  ( 1+√(V_2/V_1) - √(E/V_1) - √(E/V_1-1+V_2/V_1))




where

    A = 4π/1+√(V_2/V_1)√(V_1 V_2/ω^2).

One can then obtain S_ inst(τ), by noting that the period of the orbit, τ, is related to the energy by τ(E)=W'(E). Inverting this to find E(τ) then allows one to calculate S_ inst(τ).






While for these simple models one can obtain analytical expressions for the action, in general it must be found numerically e.g. using the ring-polymer instanton method.<cit.> For τ>τ_ c this is precisely what is already done in standard instanton calculations. 
However, for τ<τ_ c there does not exist a real periodic orbit. One might reasonably be concerned that finding a trajectory in complex positions would be impractical, however, we will now argue that this is unnecessary. 
Note first that the terms in Eq. (<ref>) containing S_ inst(τ) for τ<τ_ c are subdominant.
This allows us to approximate the action in this region without changing the key asymptotic properties of the theory.  




For example, when S”'_ inst(τ_ c)>0, we can use 

    S_ inst(τ) ≈[                                    S_ inst(τ)                                     if τ≥τ_ c; S_ inst(τ_ c) + 1/τ∑_j=1^3   s_j/j!(τ-τ_ c)^j                                     if τ<τ_ c ]

where

    s_1    = S'_ inst(τ_ c)τ_ c≡ V^τ_ c≡ S_ inst(τ_ c) 
    
        s_2    = 2S'_ inst(τ_ c)+S”_ inst(τ_ c)τ_ c
    
        s_3    = 3S”_ inst(τ_ c) + S”'_ inst(τ_ c)τ_ c.

Under the condition that S”'_ inst(τ_ c)>0, this approximation satisfies a number of key properties. First, it matches the first three derivatives of the action at τ_ c. Second, it correctly predicts that S_ inst(τ)→-∞ as τ→0. Finally, it also satisfies S_ inst(τ)<τ V^, and hence it is guaranteed that Δ S_ n>0. Importantly, as this approximation only involves information at τ≥τ_ c it is straightforward to evaluate using standard techniques.




As with the symmetric Eckart barrier the parameters for the model are chosen to correspond to previous studies<cit.> and are given in reduced units as
m=1, L=8/√(3π), 
V_1=(6/π)^1/4, 

V_2=(24/π)^1/4, 

and ħ=1. The exact results are again calculated using Eq. (<ref>), with the exact analytical P(E).<cit.>  



Figure <ref> compares the rate calculated with the new theory against the exact, classical, parabolic barrier, and standard (old) instanton results for the asymmetric Eckart barrier. We see that the new theory performs just as well for the asymmetric barrier as it does for the symmetric one. This is an important property of the standard instanton result that the present theory also retains, and stands in contrast to centroid based path-integral methods that are known to breakdown at low temperatures for asymmetric systems.<cit.> Pleasingly, we find that the results of the new theory using either the exact action, or the approximate action given by Eq. (<ref>) are graphically indistinguishable at all values of τ for this system, differing by at most 0.02%. This illustrates that the new theory can be accurately applied using just information that is available from standard ring-polymer instanton calculations. 

 The lower panel of Fig. <ref> also shows how the results converge with increasing number of instanton orbits included in the sums in Eq. (<ref>). One can immediately see that, for this range of temperatures, including terms up to n=2 already agrees almost perfectly with the full sum: the largest deviation being just over 0.2%. In fact, even including only a single orbit (n=1) is practically sufficient for nearly all temperatures considered, with the only significant deviation occurring in a very narrow range near the divergence at β/β_ c=2 [almost hidden by the frame of the graph]. Away from this divergence the largest deviation occurs close to the crossover temperature, where the full sum is just over 1% larger than the single orbit result. This is a significant result as it indicates that Eq. (<ref>), which only involves information available from a standard instanton calculation, is already sufficient for studying the behaviour of the rate for a wide range of temperatures near crossover. Finally, the narrow range of temperatures affected by the divergence at β/β_ c=2 may seem surprising when compared to the broad divergence of the parabolic barrier rate. However, this difference can be understood by noting that in the case of the parabolic barrier rate the diverging term is asymptotically dominant, whereas at β/β_ c=2 the term that is divergent is formally subdominant and has to overcome an exponential suppression.  

 




















































































 





§ DISCUSSION

Despite the limited numerical significance of the multi-orbit terms for the systems studied here, their natural appearance in the theory is theoretically interesting. For a fixed value of τ these multi-orbit terms are formally subdominant (i.e. negligible in comparison) to every term in the asymptotic series for the single orbit. Therefore, they do not contribute within a basic asymptotic analysis, be it either Poincaré asymptotics, where a fixed number of terms in the asymptotic series are included, or superasymptotics, where the number of terms is chosen to minimise the error.<cit.> The multi-orbit terms, instead, correspond to what are often referred to as hyperasymptotic contributions.<cit.>

These terms are significant because
the error made by the original asymptotic series can be written recursively in terms of them.



Hyperasymptotic series and their resummation are a central part of the closely related areas of resurgent asymptotics, exponential asymptotics, and transseries.<cit.> Such approaches allow one to go beyond the accuracy of superasymptotics and in principle to even obtain exact results. While these techniques have found a wide range of uses in modern physics, with applications in string theory, quantum field theory, and cosmology,<cit.> they have as yet not found wide use in chemistry or chemical physics. 


It will therefore be interesting in the future to make connection to these approaches, for example by using the method proposed by Berry and Howls<cit.> to derive the multi-orbit contributions to the theory.

























It should be noted that Miller's original derivation of the standard thermal instanton theory proceeded via a microcanonical expression involving multi-orbit instantons.<cit.> This raises the obvious question, how is the present theory related to Miller's original result?
While Miller's microcanonical result is equivalent in one dimension to the uniform WKB transmission probability, Miller noted immediately that it did not recover the correct result in a separable system. Despite this, once integrated asymptotically over energy it recovers the thermal instanton, which does correctly describe both separable and non-separable systems.
One might, therefore, wonder, can Eq. (<ref>) be derived from Miller's microncanonical expression following the method used in the Appendix in one dimension. 
Interestingly, the answer is no. One finds that, the resulting theory is incorrect, failing to recover the separable result except at low temperature.



Miller's original expression is not the only microcanonical instanton theory. Shortly after his first paper on instanton theory, Miller, along with Chapman and Garrett,  suggested an ad-hoc correction to the original microcanonical formula designed to correctly recover the separable limit.<cit.> 
More recently an alternative (easier to implement) method, known as the density of states microcanonical instanton, has been proposed.<cit.>

Notably, thermalising either of these microcanonical theories  by integrating asymptotically (ħ→0) over energy does recover Eq. (<ref>).[This is trivially true starting from the density of states instanton. The connection to Chapman, Garrett, and Miller's expression can be made by noting that it differs from the density of states instanton by a term of order ħ that can, therefore, be discarded.] 
It is important to stress, however, that the present theory should not be considered as an approximation to these microcanonical results. This can be seen most clearly for escape from a metastable well, for which the thermal instanton correctly describes the plateau in the rate at low temperature, whereas exact integration of the semiclassical transmission probability approaches zero.  
It should also be emphasised that, for the calculation of thermal rates the present approach is also more practical than thermalising a microcanonical theory.<cit.>


This is because integrating over energy requires information from instanton calculations at a wide range of imaginary times, including instantons whose period is greater than the thermal time. In contrast, Eq. (<ref>) only requires calculations at the temperature of interest, T, and a small number of integer multiples, nT.












































Instantons are, of course, not the only way to incorporate the effects of tunneling into the calculation of reaction rate constants. In particular path-integral theories, such as ring-polymer molecular dynamics (RPMD)<cit.> and related quantum transition state theories,<cit.> are capable of describing the transition between the high and low temperature regimes.

However, it is important to recognise that, as these approaches involve path-integral sampling, they are practically very different methods.

In particular, numerical determination of free-energy differences, and the associated need to compute the potential energy at a large number of configurations typically makes sampling methods more expensive than instanton theory. The present theory and RPMD rate theory thus have different use cases. RPMD is most useful in liquid systems where instanton theory cannot be applied,<cit.> and instanton theory is most useful for gas phase and surface reactions in combination with high-level ab-initio electronic structure theory.<cit.> 

One might reasonably suggest that making a harmonic approximation to RPMD-TST would avoid the need to perform path-integral sampling, and hence give a competitor to the present theory. However, to evaluate harmonic RPMD-TST in the vicinity of the crossover temperature one would need to derive a uniform theory that would likely be very similar to Eq. (<ref>). This is because, above the crossover temperature a basic harmonic approximation to RPMD-TST recovers the parabolic barrier rate,<cit.> and below the crossover temperature it is closely related to the F formulation of the thermal instanton.<cit.> 






§ CONCLUSIONS AND FUTURE WORK

The breakdown of instanton theory at the crossover temperature has been a long standing problem in reaction rate theory. Although suggestions had been made to overcome this problem,<cit.> none were entirely satisfactory. Here the problem has been rigorously solved for the general class of multidimensional problems in which the instanton collapses smoothly to the transition state. To derive this result we have combined a new real-time version of Richardson's flux-correlation function derivation of instanton theory<cit.> with the modern method for developing uniform asymptotic expansions due to Bleistein.<cit.>

The resulting theory is a rigorous semiclassical theory for the rate that is uniformly valid as a function of τ=βħ. Unlike previous approaches<cit.> the present result is a smooth function of temperature.
The new theory is also rigorously independent of dividing surface, and obeys detailed balance. 
Although the derivation was motivated by considering real-time trajectories, the final result only involves imaginary-time quantities (just like the original instanton theory).
Importantly, close to the crossover temperature the new theory only requires information already available from a standard instanton calculation, meaning that the theory can immediately be applied within existing instanton codes. 

There are a number of exciting avenues for further the theoretical development. First, the present theory is an important step towards the development of more accurate microcanonical instanton theories. In particular, as the present theory is a globally valid (and smooth) function of temperature, it is now possible to use the steepest descent approximation to the inverse Laplace transform to obtain the microcanonical rate at any energy.<cit.> [Note in this approach the asymptotic parameter is not associated with ħ.] 



The real-time formalism developed here also opens up exciting opportunities for future research. Although the real-time components of the trajectories do not contribute to the current theory, they could be used to extract extra information in other contexts, such as vibrationally state-resolved or electronically nonadiabatic reaction rates.





Another interesting direction for further development lies in connecting the present results to the F formulation. While the derivation used here was based on reactive scattering, the final result should be equally applicable to escape from a metastable well—the basis of the F derivation of instanton theory. Exploring this connection further would nevertheless be useful to link the current work more closely to the high-energy physics literature.











A key advantage of the first-principles derivation of the new theory is that it provides a rigorous framework for future work.


One such area is the generalisation of the theory to more complex systems in which the instanton does not collapse smoothly to the transition state with increasing temperature. 

While these systems have not received much attention so far in the chemical literature, this is likely because there has not existed a theory that can treat them.

Another advantage of the present theory is that it can be systematically improved via the inclusion of higher order terms in the asymptotic series. These perturbative corrections have already been incorporated into the ring-polymer instanton (RPI) formalism (to give RPI+PC) for the calculation of molecular tunneling splittings.<cit.>
Future work will look to implement these corrections in the context of reaction rates to help describe systems where there is a significant change in anharmonicity along the reaction coordinate.

Such a combination of RPI+PC with the present theory would be a strong competitor to SCTST,<cit.> with the important advantage of providing a rigorous description of deep tunneling.    




















































































   
















































































































§ SUPPLEMENTARY MATERIAL

The supplementary material contains derivations of Eqs. (<ref>) and (<ref>). It also contains a comparison of the present approach to the recent proposal of Upadhyayula and Pollak that illustrates that their result is not a rigorous uniform asymptotic approximation to the rate constant in τ as ħ→0.



§ DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request. 




§ ACKNOWLEDGEMENTS

I would like to thank Jeremy Richardson and George Trenins for helpful discussions. 
This work was supported by an Independent Postdoctoral Fellowship at the Simons Center for Computational Physical Chemistry, under a grant from the Simons Foundation (839534, MT).










    





    
    







§ APPENDIX: DERIVATION OF MAIN RESULT IN ONE DIMENSION FROM KEMBLE'S UNIFORM SEMICLASSICAL TRANSMISSION PROBABILITY



In the main text we derive our uniform expression in the time domain, using methods for the asymptotic evaluation of integrals.
In one dimension an alternative but equivalent approach is to use WKB analysis. These two approaches are equivalent as the asymptotic parameter is formally the same. To arrive at a uniform expression for the thermal rate we can, therefore, begin with the uniform WKB expression for the transmission probability as a function of energy

    P_ SC(E) = 1/1+e^W(E)/ħ.

This expression was originally proposed by Kemble<cit.> in 1935 and was later rigorously derived by Fröman and Fröman.<cit.> 
It involves the reduced Euclidean action, W(E), which is related to S_ inst(τ) via a Legendre transform as

    S_ inst(τ) = W(E) + τ E

with

    τ = -W'(E)

and

    E = S_ inst'(τ).


When E>V^, then W(E)<0 and we can write

    P_ SC(E) = ∑_n=0^∞ (-1)^n  e^ n W(E)/ħ.

Similarly when E<V^ then W(E)>0 and we can write 

    P_ SC(E) = e^- W(E)/ħ∑_n=0^∞ (-1)^n   e^-n W(E)/ħ.


To obtain the correct uniform asymptotic expression for the thermal rate we therefore begin by writing

    k Z_r ∼1/2πħ∫_0^∞ P_ SC(E) e^-τ E/ħdE.

Then separating this into two parts

    kZ_r∼1/2πħ   ∫_0^V^ P_ SC(E) e^-τ E/ħdE 
    
        +1/2πħ   ∫_V^^∞ P_ SC(E) e^-τ E/ħdE

allows us to insert the expansions Eqs. (<ref>) and (<ref>) to give

    kZ_r∼1/2πħ   ∫_0^V^∑_n=1^∞ (-1)^n+1  e^-n W(E)/ħ-τ E/ħdE 
    
        +1/2πħ   ∫_V^^∞∑_n=0^∞ (-1)^n  e^ n W(E)/ħ-τ E/ħdE.


Consider first the integral up to the barrier height

    I_n(ħ)=1/2πħ∫_0^V^ e^-(n W(E)+τ E)/ħdE.

Integrating by steepest descent gives the stationary condition as

    n W'(E_n^⋆(τ)) = -τ.

Now for τ/n<τ_ c this stationary point will move outside of the integration range, E_n^⋆(τ<nτ_ c)>V^. Hence to obtain a uniform expression valid when E_n^⋆>V^, we again make use of Bleistein's method.<cit.> Hence, defining

    S_n,inst(τ) = n W(E_n^⋆)+τ E_n^⋆ = n  S_inst(τ/n)

and

    Δ S_n = τ V^ - S_n,inst(τ)

application of Bleistein's method results in the following expression as ħ→0

    I_n(ħ)   ∼ k_n,instZ_r1/2erfc(sgn(nτ_ c-τ)√(Δ S_n/ħ)) 
       +1/2πe^-τ V^/ħ1/nτ_ c - τ
       - 1/2πe^-τ V^/ħsgn(nτ_ c-τ) √(-S_n,inst”(τ)/2Δ S_n)

where

    k_n,instZ_r = e^-S_n,inst(τ)/ħ1/√(2πħ)(-d^2 S_n,inst/dτ^2)^1/2.



Next we consider the integral above the barrier height

    J_n(ħ)=1/2πħ∫_V^^∞ e^(n W(E)-τ E)/ħdE.

When evaluating this asymptotically, we note that W'(E) < 0. 
Hence, for τ>0 then there is no stationary point inside the integration region and the integrand is peaked about E=V^. The integral is, therefore,  approximated asymptotically by expanding the exponent linearly. Using,
W'(V^)=-τ_ c=-2π/ω we thus have

    J_n(ħ)    ∼1/2πħ∫_V^^∞ e^-(nτ_ c(E-V^)+τ E)/ħdE 
        = 1/2π e^-τ V^/ħ1/nτ_ c+τ

as ħ→0.

Combining these two results together we then obtain

    k Z_r   ∼∑_n=1^∞ (-1)^n+1 k_n,instZ_r1/2erfc(sgn(nτ_ c-τ)√(Δ S_n/ħ)) 
       -∑_n=1^∞ (-1)^n+11/2πe^-τ V^/ħsgn(nτ_ c-τ) √(-S_n,inst”(τ)/2Δ S_n)
       +∑_n=1^∞ (-1)^n+11/2π e^-τ V^/ħ1/nτ_ c-τ
       +∑_n=0^∞ (-1)^n1/2π e^-τ V^/ħ1/nτ_ c+τ.

This can then be simplified by noting that the final two sums can be combined to give

    ∑_n=1^∞ (-1)^n+11/2π n/ω-τ + ∑_n=0^∞ (-1)^n1/2π n/ω+τ=∑_n=-∞^∞(-1)^n/τ+2π n/ω

which is equivalent to 

    ∑_n=-∞^∞(-1)^n/τ+2π n/ω = ω/2sin(ωτ/2).

Hence, we arrive at our final expression for the uniform thermal rate in one dimension

    k Z_r   ∼∑_n=1^∞ (-1)^n+1 k_n,instZ_r1/2erfc(sgn(nτ_ c-τ)√(Δ S_n/ħ)) 
       -∑_n=1^∞ (-1)^n+11/2πe^-τ V^/ħsgn(nτ_ c-τ) √(-S_n,inst”(τ)/2Δ S_n)
       +ω/4πsin(ωτ/2) e^-τ V^/ħ.







