Tunneling splittings of vibrationally excited states using general instanton paths
aa r X i v : . [ phy s i c s . c h e m - ph ] A ug Tunneling splittings of vibrationally excited states using general instantonpaths
Mihael Eraković and Marko T. Cvitaš a) Department of Physical Chemistry, Ruđer Bošković Institute, Bijenička Cesta 54, 10000 Zagreb,Croatia (Dated: 27 August 2020)
A multidimensional semiclassical method for calculating tunneling splittings in vibrationally excited statesof molecules using Cartesian coordinates is developed. It is an extension of the theory by Mil’nikov andNakamura [
J. Chem. Phys. , 124311 (2005)] to asymmetric paths that are necessary for calculatingtunneling splitting patterns in multi-well systems, such as water clusters. Additionally, new terms areintroduced in the description of the semiclassical wavefunction that drastically improve the splittingestimates for certain systems. The method is based on the instanton theory and builds the semiclassicalwavefunction of the vibrationally excited states from the ground-state instanton wavefunction along theminimum action path and its harmonic neighborhood. The splittings of excited states are thus obtainedat a negligible added numerical effort. The cost is concentrated, as for the ground-state splittings, in theinstanton path optimization and the hessian evaluation along the path. The method can thus be appliedwithout modification to many mid-sized molecules in full dimensionality and in combination with on-the-flyevaluation of electronic potentials. The tests were performed on several model potentials and on the waterdimer.The following article has been submitted to Journal of Chemical Physics. After it is published, itwill be found at https://aip.scitation.org/journal/jcp
I. INTRODUCTION
Tunneling splittings of molecular energy levels arespectroscopic signatures of rearrangements that takeplace between degenerate symmetric wells via tunnelingmotion . These splittings can be detected in high-precision spectroscopic measurements and carry in-formation about the molecular structure and dynamicsalong the accessible tunneling paths . Dynamical theo-ries, in combination with potential energy surfaces (PES)or first principles electronic structure calculations, aim toreach an agreement with the measurements and providea physical interpretation .Computational studies of tunneling splittings concen-trated initially on the symmetric tunneling systems. Pro-ton transfer in malonaldehyde , collective migration ofhydrogen atoms in ammonia or concerted monomer mo-tion in the HF dimer are some examples of extensivelystudied systems. More recently, the splitting patterns inwater clusters have also come into focus, motivated bythe development of a universal water model that is ca-pable of predicting properties of liquid water from firstprinciples . Water clusters are multi-well systemsand exhibit multiple tunneling pathways. These tunnel-ing paths are often asymmetric, whereby tunneling atomstake on different roles in the minima they connect .The splittings vary over many orders of magnitudeeven in a single system. In water dimer, for instance, a) Author to whom correspondence should be addressed:[email protected] they vary over three orders of magnitude dependingon which of the five tunneling pathways is taken, all ofwhich reflect on the appearance of the splitting patternin the spectrum. Likewise, the experiments on watertrimer and pentamer show that the splittings of vi-brationally excited states differ by up to three orders ofmagnitude in comparison to the ground-state splittings,depending on which normal mode is excited. The in-terplay of different rearrangement pathways can lead toan increase in the width of a vibrational manifold and areduction in another , as contributions from differentpathways enter the splitting pattern with the same or op-posite signs, respectively. Qualitatively different tunnel-ing splitting patterns in water hexamer spectrum distin-guish the prism and cage structures of almost equalenergy. The contributions of different tunneling path-ways can be disentangled, by computation, to reveal theexperimental evidence of unexpected mechanisms, suchas the simultaneous double hydrogen-bond breaking inthe water hexamer prism. The investigations of tunnel-ing splitting patterns thus provide a sensitive test of boththe dynamical theories and the potentials at geometriesalong which the hydrogen bonds rearrange.Tunneling splittings can be determined by solv-ing the Schrödinger equation. Variational methodshave been used to determine the tunneling splittingsin, e.g., HF dimer , ammonia , vinyl radical ,malonaldehyde and water dimer , using time-independent methods, and, e.g., malonaldehyde ,using time-dependent methods. Both, ground- andexcited-state splittings are obtained in this way, how-ever, the cost of these methods scales prohibitively withthe basis set size and a different approach is neededfor larger systems. Diffusion Monte Carlo in combina-tion with the projection operator techniques has beenused to calculate tunneling splittings in water trimer and malonaldehyde . The recently-developed path-integral molecular dynamics method has been used toobtain the splittings in water trimer and hexamer infull dimensionality. However, the tunneling splittings ofvibrationally excited states, which are the the topic ofour investigations here, cannot be obtained using theseapproaches. The remaining options include resorting todynamical approximations , reduced-dimensionalityapproaches or semiclassical methods .The development in this paper belongs to theclass of semiclassical methods based on the in-stanton theory . In the standard instantonformulation , tunneling splitting is calculated from thezero-temperature limit of the quantum partition functionin the path-integral formalism. The dominant contribu-tion to the partition function comes from the minimumaction path (MAP) that connects the symmetry-relatedminima. The contribution from all other paths is esti-mated analytically using the parameters in a harmonicexpansion of the potential in the directions perpendic-ular to the MAP. Instanton theories of tunneling split-tings come in several variants. Some approaches useapproximate MAPs , determined from the stationarypoints on the PES, and approximate hamiltonians ,in which analytic expressions for vibrational couplingsare fitted to the PES. The present contribution belongsto the category that is based on the numerically exactMAPs. Mil’nikov and Nakamura use the exact MAPand Hessians along the MAP to obtain splittings via theintegration of Jacobi fields (henceforth reffered to as theJFI method). They employ internal coordinates in theirtreatment in order to separate the overall rotational mo-tion. Ring-polymer instanton (RPI) method likewiseuses the numerically exact MAP and evalulates the split-ting from the eigenvalues of the discretized functionaldeterminant of the action Hessian. This approach istherefore computationally more demanding than the JFImethod and recovering the rotational dependence ofthe splittings, when it is significant, becomes elaborate .Its advantage is that it can be applied without modifica-tion to any molecule of interest, as it works in Cartesiancoordinates, and it can readily be applied to systems thatexhibit asymmetric MAPs. The RPI method featuredprominently in the recent calculations of tunneling split-ting patterns in water clusters. It was used to obtainthe ground-state tunneling splitting pattern and revealmachanisms responsible for its formation in asymmetricsystems such as the water dimer, trimer , hexamer and octamer in full dimensionality.Standard instanton approaches for calculating tun-neling splittings suffer from the same drawback as theMonte-Carlo and path-integral based method mentionedabove in that they cannot provide the splittings of vibra-tionally excited states from the outset. It is well-knownthough that the instanton expression for the ground-state tunneling splitting can be obtained using a variant of theWKB theory and Herring formula . This link thusprovides a consistent route for calculating tunneling split-tings of vibrationally excited states , where this paperaims to contribute. In fact, the semiclassical methodsbased on the wavefunction along the classical trajectorythat connects the minima on the inverted potential en-ergy surface (PES), i.e., along the MAP, are regularlyreferred to as the instanton methods in literature .Tunneling splittings of vibrationally excited states havebeen obtained using the related methods in symmet-ric systems such as malonaldehyde , tropolone , 9-hydeoxyphenalenone , HO , formic acid dimer andthe vinyl radical .In our recent work , we generalized the JFI approachof Mil’nikov and Nakamura to obtain the ground-statetunneling splittings for asymmetric paths in Cartesian co-ordinates. We obtained an almost perfect agreement be-tween the JFI and RPI splittings for systems in whichrotations do not couple strongly to the internal degreesof freedom, like water trimer or malonaldehyde. The de-velopment enabled us to treat large asymmetric systemsthat exhibit slow motion of a heavy-atom skeleton, suchas the water pentamer , in full dimensionality. We wereable to calculate the 320-level ground-state splitting pat-tern of the pentamer, including the state symmetries, andto identify rearrangement motions responsible for its for-mation, in a treatment which would become extremelycumbersome in the RPI approach due to the large imag-inary time periods involved.Motivated by the effectiveness of our JFI approach,the present work aims to derive the tunneling splittingsof vibrationally excited states for general, symmetric andasymmetric paths, in a consistent approach. This is ac-complished by a WKB construction of wavefunction thatreproduces our JFI result in the ground state. In essence,our approach below follows the work of Mil’nikov andNakamura in which they extend their ground-state in-stanton theory of Ref. 52 to treat the low-lying vibra-tionally excited states. Distinctly, in our approach wecan readily treat asymmetric paths, that are regularly en-countered in the studies of clusters, and we again work inCartesian coordinates in order to make our approach gen-eral. Unlike Ref. 53, we treat the ‘longitudinal’ modes,that are parallel to the MAP at minima, and ‘transver-sal’ modes, that are perpendicular to the MAP at min-ima, on an equal footing. We achieve this by using adifferent form of the matching wavefunction near min-ima, which allows for a displacement of the wavefunctionnode away from the MAP. In particular, this means thatwe can treat the asymmetric paths in which the excitedmode is the longitudinal mode at one minimum and isa transversal mode near the other end of the MAP. Thestraightforward generalization of Ref. 53 to asymmetricpaths would give a zero splitting in that case. The the-ory thus includes newly added terms which for certaincases dramatically improve the splitting estimates even insymmetric systems. It is applicable to low vibrationallyexcited states.Instanton method evaluates the splittings with a mod-est number of potential evaluations (on the order of athousand) in comparison with the exact methods .This means that the computations can be performed onlarger systems or using more accurate electronic poten-tials. In certain circumstances, it can probably providethe best possible splittings in a compromise between theaccuracy of the dynamical theory and the level of elec-tronic structure theory that the dynamical treatment al-lows. Numerical effort is concentrated in the MAP opti-mization and the Hessian evaluation along the MAP .Since the calculations of splittings in vibrationally ex-cited states do not require any additional informationabout the molecular system, they too enjoy the same ad-vantages over the exact methods.The paper is organized as follows. In Section II, we usea semiclassical expansion to approximate the wavefunc-tion about the MAP. The wavefunctions that start fromthe ‘left’ and from the ‘right’ symmetry-connected min-ima along the MAP are constructed and used in Herringformula at the dividing surface to obtain the ground-state tunneling splitting, which is identical in form tothe JFI instanton expression from our previous work .The derivation follows Ref. 52, but does not assume themirror symmetry of the potential along the MAP. Weprove explicitly that the expression for the splitting doesnot depend on the position of the connection point be-tween the left- and right-localized wavefunctions alongthe MAP. Section II thus lays the groundwork for con-structing the wavefunctions of the excited states in Sec-tion III. Section III follows the work of Ref. 53, but ar-rives at a different expression for the tunneling splittingsof vibrationally excited states. As stated above, our for-mulation treats longitudinal and transversal excitationsin a unified approach. In certain cases, as the numer-ical exercises on symmetric and asymmetric model po-tentials in Section IV show, the contribution from thenewly added terms can dominate the splittings. Thedeuterated water dimer provides a real-life test systemthat exhibits asymmetric paths, including the path fea-turing the longitudinal-transversal excitation mode andthe vibrational modes that do not line up either in eitherparallel or perpendicular direction with respect to theMAP near minima. The importance of different termsin the semiclassical expansion is discussed in terms ofthe accuracy improvements that they bring to the split-tings and the stability with regards to the position of thedividing surface. Conclusions and outlook are given inSection V. Atomic units (¯ h = 1) are used throughoutunless indicated otherwise. II. GROUND-STATE TUNNELING SPLITTING
Tunneling splittings in molecular systems with mul-tiple symmetry-related minima can be expressed as theeigenvalues of a tunneling matrix in which rows and columns are numbered by the indices of the minima, us-ing group theoretic arguments. The tunneling matrixelement h connecting two minima, termed L and R forconvenience, is the transition amplitude between the de-generate states φ (L) and φ (R) , localized in their respec-tive wells, that neglect the presence of tunneling motion.The tunneling splitting of the isolated double-well systemconnecting minima L and R is thus ∆ = − h , the dif-ference between the tunneling matrix eigenvalues. Thetunneling matrix eigenvectors are comprised of the co-efficients of the energy eigenstates in the φ (L / R) basis.For a double-well system, they form the symmetric andantisymmetric linear combinations of φ (L) and φ (R) .In our previous work , we derived the tunneling ma-trix element h , or equivalently the tunneling splitting ∆ , using the JFI theory. The splitting is dominated bythe Euclidean action of the MAP, while the contribu-tions from all other paths in the harmonic neighborhoodof the MAP are collected into the fluctuation prefactor.The fluctuation prefactor is then evaluated via integra-tion of Jacobi fields . We now proceed along the linesof Refs. 53, 57, and 70 to derive an identical expressionusing the semiclassical WKB approach to construct thelocalized states φ (L / R) .Whenever the energy eigenstates are well approxi-mated by the symmetric and antisymmetric combina-tions of the localized state functions, φ (L / R) , the tunnel-ing splitting can be calculated using Herring formula , ∆ = R (cid:0) φ (R) ∂∂S φ (L) − φ (R) ∂∂S φ (L) (cid:1) δ ( f ( x ))d x R (cid:12)(cid:12) φ (L) (cid:12)(cid:12) d x , (1)where x is the molecular geometry in mass-scaled Carte-sian coordinates and f ( x ) = 0 is an implicit equation ofan arbitrary dividing plane, which separates the two min-ima. Variable S corresponds to the position on a localnormal to the dividing plane.We now construct the localized states φ (L / R) in thefamiliar WKB form as φ = e − h ( W + W ¯ h ) , (2)where we drop the labels (L/R) from this point onwardsas the equations are valid in both wells. In Eq. (2), W satisfies Hamilton-Jacobi equation ∂W ∂x i ∂W ∂x i = 2 V ( x ) , (3)where V ( x ) is the PES, and W satisfies the transportequation, ∂W ∂x i ∂W ∂x i − ∂ W ∂x i ∂x i + E = 0 . (4)We note here that E is approximated by the ground-stateenergy of the quantum harmonic oscillator and is of theorder ¯ h . The whole energy dependence is moved to thetransport equation, Eq. (4), following Ref. 52 and 57.Hamilton-Jacobi equation, Eq. (3), can be solved us-ing the method of characteristics that we briefly describein Appendix A. The characteristics of Hamilton-Jacobiequation are given by ¨ x ( τ ) = ∇ V ( x ( τ )) , (5)with τ as parameter. The form of Eq. (5) suggests thatthe characteristics represent classical trajectories on theinverted PES and that τ represents time. As shown inAppendix A, these trajectories must have zero energy inorder to satisfy Eq. (3). On a characteristic, W can beobtained by a simple integration, W (x( τ )) = W (x( τ )) + Z τ τ p ( τ )d τ, (6)where p = √ V corresponds to the mass-scaled momen-tum on the classical trajectory. It is convenient to chooseone point to correspond to the minimum of the PES anddefine W ( x min ) = 0 . The reason behind this choice isthat in the vicinity of the minimum, the wavefunctioncan then be matched to that of the harmonic oscillator,which will be used later on to determine its norm. Withthat choice, since the minimum on the PES is a max-imum on the inverted PES, all other points along thecharacteristic correspond to time τ > τ min and the inte-gral in Eq. (6) remains positive. However, by choosingthe first point at the minimum, the time to any otherpoint will be infinite, since it takes infinite time to moveaway from the minimum with zero energy. This presentsa problem in a numerical implementation, which can con-veniently be fixed by reparametrizing the characteristicsusing the arc length distance S from the minimum alongthe characteristic, d S d τ = r d x i d τ d x i d τ = p . (7)Using this transformation, Eq. (6) reduces to W ( x ) = Z S ( x )0 p ( S ′ )d S ′ . (8)We observe that W equals Jacobi action between theminimum and the point S on the characteristic. Thecharacteristic between the minimum and a point x , aswell as W , can both be determined by a Jacobi actionminimization. The gradient of W is therefore parallel tothe characteristic.In order to describe W in the vicinity of a given char-acteristic, we assume that the Hessian of the potential, H ( S ) , along the characteristic is known. The equationfor the Hessian of W , A ij = ∂ W ∂x i ∂x j , along a character-istic is then obtained, by differentiating Eq. (3) twice,as p ∂∂S A ( S ) = H ( S ) − A ( S ) . (9) Riccatti equation in Eq. (9) is identical to the equationthat emerges in the JFI method as the equation forthe log-derivative of a Jacobi field. The initial conditionfor Eq. (9) at the minimum, where p = 0 , is A = H (0) / . This identification later serves to match thesemiclassical wavefunction φ in Eq. (2) to that of theharmonic oscillator at the minimum.We can now expand W around the characteristic as W ( S, ∆ x ) = Z S p ( S ′ )d S ′ + 12 ∆ x ⊤ A ∆ x , (10)where { S, ∆ x i } is a set of local coordinates for an ar-bitrary point x . Coordinate S corresponds to the posi-tion of the point x on the characteristic which satisfies ( x i − x i ) p i = 0 . The coordinates ∆ x i define an orthog-onal shift from x to x , so that ∆ x i = x i − x i . Jacobianof the transformation is derived in Appendix A. The firstterm in the expansion is missing, since ∇ W is tangent tothe classical trajectory. Eq. (10) serves to describe W inthe neighborhood of the characteristic without the needto compute new characteristics.Transport equation in Eq. (4) can be solved on a char-acteristic by a simple integration W ( S ) = 12 Z S Tr ( A ( S ′ ) − A ) p d S ′ , (11)where we inserted the energy of harmonic oscillator E = Tr A into the expression. Using Eqs. (10) and (11),the localized wavefunctions in Eq. (2) take the followingforms in their respective wells, φ (L) ( S ) =e − R S p ( S ′ )d S ′ − R S A (L)( S ′ ) − A (L)0 ) p d S ′ e − ∆ x ⊤ A (L) ∆ x φ (R) ( ˜ S ) =e − R ˜ S p ( ˜ S ′ )d ˜ S ′ − R ˜ S A (R)( ˜ S ′ ) − A (R)0 ) p d ˜ S ′ e − ∆ x ⊤ A (R) ∆ x , (12)where S is the distance from the left minimum alongthe characteristic, while ˜ S denotes the corresponding dis-tance from the right minimum. In the harmonic regionsnear minima, these wavefunctions are matched to that ofthe quantum harmonic oscillator, as we describe in Ap-pendix B. From that identification, we obtain their normas Z | φ | d x = s π N det A . (13)Having obtained the localized wavefunctions, Eqs. (12)and (13), we are ready to compute the tunneling splittingvia Herring formula in Eq. (1). One could take an arbi-trary dividing surface and compute the surface integralin Eq. (1) numerically. However, this requires computingthe characteristics that connect the minima with everypoint at which the integrand is evaluated on the dividingsurface. An economical way to compute the integral isto choose one point on the dividing surface and use Tay-lor expansion of W around it to evaluate the integrandat other points. If the dividing surface is chosen to bea hyperplane and the gradient of W taken to be con-stant, the integral can be computed analytically. Sincethe integrand in Herring formula is proportional to theproduct φ (L) φ (R) , the integral will be best approximatedif the point on the dividing surface is chosen so that itmaximizes this product. This is equivalent to the mini-mization of Z S (L) p (L)0 ( S ′ )d S ′ + Z ˜ S (R) p (R)0 ( ˜ S ′ )d ˜ S ′ , (14)which is accomplished when the point lies on the clas-sical trajectory that connects the two minima. In thatcase, the characteristics that originate at two minima aresmoothly joined at the connection point S = S cp and ˜ S = S tot − S cp , where S tot is the total length of the MAPthat connects the two minima. The two joined charac-teristics coincide with the instanton trajectory . Thesum of W (L)0 and W (R)0 then becomes the Jacobi action ofthe instanton trajectory, W (L)0 + W (R)0 = R S tot p d S . Thedividing surface is taken to be orthogonal to the trajec-tory at the connection point and Herring formula givesthe ground-state tunneling splitting as ∆ = r det A π N e − R S tot0 p d S − W (L)1 − W (R)1 Z ∂W (L)0 ∂S − ∂W (R)0 ∂S ! e − ∆ x ⊤ A (L)+ A (R)2 ∆ x δ ( f ( x ))d x , (15)where ∂W (R)0 ∂S = − ∂W (R)0 ∂ ˜ S evaluates to p at the connectionpoint, and is kept constant in the surface integral.In order to solve the integral in Eq. (15), we note thatthe matrix ¯ A = A (L) + A (R) (16)possesses a zero eigenvalue, which corresponds to thetangent vector. This is easily proved by differenti-ating Hamilton-Jacobi equation, Eq. (3), which yields A (L / R) p (L / R)0 = ∇ V . Subtracting these two equationsand using the fact that p (L)0 = − p (R)0 , it follows that ( A (L) + A (R) ) p (L)0 = 0 . The eigenvectors of ¯ A whichcorrespond to its non-zero eigenvalues λ i then span thedividing surface. Transforming to the eigenvector basisreduces this integral to ∆ = 2 p r det A π N e − R S tot0 p d S − W (L)1 − W (R)1 Z e − λ i ξ i d ξ = 2 p r det A π det ′ ¯ A e − R S tot0 p d S − W (L)1 − W (R)1 , (17) where det ′ denotes the product of non-zero λ i ’s, and W (L / R)1 at S = S cp are calculated using Eq. (11). Theground-state tunneling splitting formula in Eq. (17) isidentical to the instanton formula, Eq. (33) in Ref. 66.The splitting in Eq. (17) does not depend on the posi-tion of the connection point on the instanton trajectory.This is evident from the derivation of Ref. 66, but thepresent treatment does not guarantee it and we prove itin Appendix D. III. EXCITED-STATE TUNNELING SPLITTING
The calculation of tunneling splittings in vibrationallyexcited states is approached in a consistent manner,following Ref. 53. We assume one quantum of vibra-tional excitation in the mode with frequency ω e and con-struct the WKB wavefunctions in Eq. (2) by solving theHamilton-Jacobi and transport equations, Eqs. (3) and(4), and finally insert them into Herring formula, Eq. (1),which remains valid for the excited states.Only the transport equation depends on the energyand is different for the excited state. We decompose W in form W = W (0)1 + w, (18)where W (0)1 is the ground-state function given byEq. (11), and insert Eq. (18) in Eq. (4). We then findthat w satisfies ∂W ∂x i ∂w∂x i + ω e = 0 . (19)In a crucial difference from Ref. 53, we seek the solutionof Eq. (19) along the characteristic in the following form w = − ln (cid:0) U ⊤ ∆ x + F (cid:1) . (20)The above form, when used in Eq. (2), allows the match-ing to a harmonic oscillator wavefunction in the neighbor-hood of minima for both, the longitudinally and transver-sally excited modes with respect to the MAP, in a uni-fied approach. We insert Eq. (20) into Eq. (19), multiplythrough with U ⊤ ∆ x + F and equate the terms of order ∆ x and ∆ x to obtain equations for F and U as p dd S F = ω e F, (21) p dd S U = ω e U − AU + 2 (cid:0) U ⊤ p − ω e F (cid:1) Ap p . (22)Eq. (22) can be simplified by noting that, by definition,components of U equal to U i = ∂∂x i e − w = ∂∂x i F, (23)where the second equality is due to the fact that thepartial derivative is taken on the characteristic. Thismeans that the projection of U onto the tangent is U ⊤ p = ∂F∂x i ∂W ∂x i = p dd S F = ω e F, (24)where Eq. (21) was used. Combining Eqs. (24) and (22),reduces the equation for U to p dd S U = ω e U − AU . (25)This is the same equation that Mil’nikov and Nakamura obtained in their treatment of transversal excitations.Here, however, we use it for both, longitudinal andtransversal excitations. As our test calculations belowdemonstrate, it is important to propagate both compo-nents of U simultaneously for best accuracy.Eqs. (21) and (25) have singularities at the minima ofPES. In order to avoid them, we need to start the propa-gation a small distance ε away from the minimum alongthe characteristic. If this distance is sufficiently small tofall into the harmonic region around the minimum, theinitial conditions at ε can be taken in form F ( ε ) = U ⊤ ( x ( ε ) − x (0)) , (26)as justified in Appendix B, and U ( ε ) = U , (27)where U is the excited normal mode at the minimum.Alternatively, we can solve Eq. (25) in the region [0 , ε ] using the same procedure that was used for solv-ing Eq. (9) in Refs. 52 and 66. We expand p , A and U around minimum as p = p (1)0 S, A = A + A S, U = X i C ( i ) S i . (28)We then insert Eq. (28) into Eq. (25) and equate theterms of the same order in S i to obtain the recurrencerelation for C ( i ) , A C (0) = ω e C (0) , (cid:16) A + ( ip (1)0 − ω e ) I (cid:17) C ( i ) = − A C i − . (29)Once U has been determined, F can be obtained fromEq. (24) as F ( S ) = Z S U ⊤ ( S ′ ) t ( S ′ )d S ′ = U ⊤ ( S ) p ( S ) ω e , (30)where t = p /p is the tangent vector at instanton tra-jectory. In this way, the anharmonicity of the PES nearminima is accounted for by A . Having obtained U ( ε ) , Eq. (25) is readily solved by a simple integrator, such asthe Runge-Kutta method .At the dividing plane, the wavefunction of the excitedstate in Eq. (2) takes the form φ (L) = (cid:16) U (L) ⊤ ∆ x + F (L) (cid:17) e − R S p ( S ′ )d S ′ e − R S A (L)( S ′ ) − A (L)0 ) p d S ′ − ∆ x ⊤ A (L) ∆ x ,φ (R) = (cid:16) U (R) ⊤ ∆ x + F (R) (cid:17) e − R ˜ S p ( ˜ S ′ )d ˜ S ′ e − R ˜ S A (R)( ˜ S ′ ) − A (R)0 ) p d ˜ S ′ − ∆ x ⊤ A (R) ∆ x . (31)By matching the above wavefunction to that of the har-monic oscillator at a minimum, one obtains the norm as Z | φ | d x = s π N det A ω e . (32)Wavefunctions in Eq. (31) are then inserted into Herringformula and the surface integral evaluated in a similarmanner to the ground-state case. This gives the tunnel-ing splitting of vibrationally excited states as ∆ = ∆ (2 ω e ) (cid:18) F (L) F (R) + 12 U (L) ¯ A − U (R) (cid:19) . (33)Since ¯ A possesses a zero eigenvalue, ¯ A − in Eq. (33) de-notes a pseudoinverse of ¯ A , defined by ¯ A ¯ A − = ¯ A − ¯ A = P , where P = I − tt ⊤ is a projector onto the orthogo-nal plane. The pseudoinverse has the same eigenvectorsas ¯ A , while its nonzero eigenvalues are reciprocals of theeigenvalues of ¯ A .It turns out, the tunneling splitting formula in Eq. (33)is dependent on the position of the connection point atwhich the dividing surface and the instanton trajectorycross. This undesirable behavior, which was not presentin the ground-state formula in Eq. (17), arises from the U (L) ¯ A − U (R) term, as shown in Appendix D. It can fur-ther be shown, by a similar analysis, that the terms whichcause this dependency cancel out if the next order termis included in the Taylor expansion of exp( − w ) , w = − ln (cid:18) F + U i ∆ x i + 12 Z ij ∆ x i ∆ x j (cid:19) . (34)However, the inclusion of Z in Eq. (34) brings new termsthat are again do depend on the connection point andto eliminate their dependence on S cp , it would be nec-essary to include higher order terms in the wavefunctionexpansion Eq. (2), such as the W term. The root of theproblem is that the expansion of exp( − w ) is inconsistentwith the expansion of W , as it gives rise to terms ofall orders in ∆ x in the expansion of w . Excluding thehigher order terms of w in Eq. (34), on the other hand,would degrade the quality of matching with the harmonicoscillator near minima.In fact, any improvement of the accuracy of the WKBwavefunction through the inclusion of extra terms in W and W necessarily requires the calculation of higher or-der derivatives of potential along the path. Calculationof the tensor of third derivatives of potential along thepath allows us to expand W in Eq. (10) up to the ∆ x term, W (0)1 in Eq. (11) up to ∆ x , and to include the ∆ x term in Eq. (34). The tensor of fourth derivativesof potential allows for the correction of the vibrationalenergy, the inclusion of the ∆ x term of W and thehigher order terms in W , W (0)1 and w . The calculationof higher order derivatives of the potential quickly be-comes computationally unfeasable for realistic potentialenergy functions and, in most cases, does not improvethe results significantly.In order to study the effect of anharmonicity that orig-inates from the inclusion of third derivatives of potentialon the tunneling splittings in numerical tests below, wederive the equation for Z along a characteristic in Ap-pendix C. It turns out that from all terms that can becomputed using the third derivatives of potential, this isthe only term that is meaningful to include in the tun-neling splitting formula, Eq. (C14), below. The inclu-sion of ∇ W (0)1 does not appreciably influence the results,whereas the inclusion of the ∆ x term in W in Eq. (10)does not result in convergent integrals on the dividingsurface.It can be shown, by using the Z contribution to thesplitting, derived in Appendix C, that the connectionpoint is best placed in the middle of the instanton pathfor symmetric systems, i.e., at the top of the barrier, be-cause, at this place, the Z contribution is the smallest.We found no such justification for the placement of theconnection point in asymmetric systems, so the safestplace to set it is at the barrier maximum as well.Alternatively, we can discard the terms that are re-sponsible for the connection point dependence of thesplittings in order to obtain an unambiguous formula-tion. For this purpose, we decompose the vector U intolongitudinal and transversal parts as U = U ⊥ + F ′ t , (35)where U ⊥ is the component of U that is perpendicularto the path. Since only the U ⊥ components contributeto the splitting in the U (L) ¯ A − U (R) term in Eq. (33),due to the fact that the tangent vector is an eigenvectorof ¯ A − with zero eigenvalue, it can be used instead ofthe vector U . We carry out the separation in Eq. (35) at S = ε and propagate U ⊥ and F independently towardsthe connection point from both minima. It can be shownthat U ⊥ satisfies the following equation p dd S U ⊥ = ω e U ⊥ − AU ⊥ − ω e F d t d S . (36)If we neglect the last term in Eq. (36), U ⊥ satisfies thesame equation as U . Vector U ⊥ remains perpendicu-lar to the instanton path , when it is propagated us-ing Eq. (25), and, as Appendix D shows, the splittingbecomes independent of the position of the connection point. Since the neglected term is proportional to thecurvature of the instanton path, it can safely be neglectedfor paths with small curvatures. For paths with a largecurvature, it turns out in Section IV, it is better to workwith the full vector U , as the deviations in the splittings,when the connection point is moved along the instantonpath, are smaller than the error introduced by the aboveapproximation. IV. NUMERICAL TESTS
We now perform tests of the above theory on a two-dimensional (2D) symmetric system, a 2D asymmetricsystem and the deuterated water dimer. Each calculationof the tunneling splitting in a vibrationally excited stateis preceded by a calculation of the ground-state tunnelingsplitting using the JFI method of Ref. 66. A JFI calcu-lation starts by an action minimization, using the stringor quadratic string method , followed by the evalua-tion of Hessians along the MAP, and, finally, it ends withthe computation of A by solving the Riccatti equation inEq. (9) along the MAP. Excited-state calculations addi-tionally require a propagation of U along the MAP usingEq. (25) for each vibrationally excited state of interest.In our tests below, we also evaluate Z along the MAP inorder to check the accuracy and convergence of the ob-tained results. To accomplish this, we first compute thetensor of third derivatives of potential along the MAP,we then use it to propagate Eq. (C2), and, finally, use B ,as well as A and U , to propagate Z along the MAP usingEq. (C11). The splittings are evaluated using Eqs. (17),(33) and (C14).In the tests, we discretized all instanton paths us-ing 600 equally spaced beads (or points) in mass-scaledCartesian coordinates and used the string method ofRef. 67 for the optimization of MAP. In the tests on waterdimer, the orientations of end beads were adjusted dur-ing optimization by minimizing the distance to the firstneighbor bead at every iteration . Convergence criterionwas taken to be the maximum value of the action gradi-ent orthogonal to the string as max (cid:8) S ⊥ i (cid:9) < − a . u . . Alarge number of beads and a tight convergence criterionwere used to ensure that the results do not depend onthe accuracy of the MAP. Hessians and third-derivativetensors were computed at all beads using fourth-order fi-nite difference method with the grid spacing of − a . u . .In water dimer calculations, we projected out the over-all translations and rotations, as described in Ref. 72.Molecular geometries, potential, Hessian matrix elementsand third derivative tensor elements were all interpolatedwith respect to the mass-scaled arc length distance S along the MAP using natural cubic splines. Eqs. (9),(C2) and (C11) were solved on the interval [0 , ε ] by lin-earization, as described previously in Ref. 52 and 66 andin Appendix C, while on the interval [ ε, S cp ] , they wereintegrated using Runge-Kutta method with the fixedstep length of − m / a . The parameter ε was takenas ε = 0 . m / a in all test systems.The normal modes were calculated at one minimumand obtained at the other minimum by utilizing the sym-metry operation that connects them in order to avoid signambiguity. Eq. (25) was then solved on the interval [0 , ε ] using the recurrence relation, Eq. (29). Taylor series of U in Eq. (28) was cut when the change in the norm of U ( ε ) fell below the threshold value of − . On theinterval [ ε, S cp ] , we used the exponential propagator tosolve Eq. (25), U ( S + h ) = e ( ω e I − A ) hp U ( S ) , (37)with fixed step length h = 10 − m / a . F values werecomputed from the tangent projection of the U vector,using Eq. (24). That procedure was found to be lesssensitive to the value of F ( ε ) than the direct integrationof Eq. (21), in Eq. (39). A. SYMMETRIC DOUBLE-WELL 2D POTENTIAL
We first test the theory on a model 2D double-well sys-tem. We call the system symmetric, since the potentialalong the MAP connecting two minima has a left-rightmirror symmetry with the maximum of the potential inthe middle of the path. The potential is given by thefollowing equations, V ( x ) = V V V + V ,V ( x ) = 12 (cid:16) x − x (1) (cid:17) ⊤ U (cid:18) α α (cid:19) U ⊤ (cid:16) x − x (1) (cid:17) ,V ( x ) = 12 (cid:16) x − x (2) (cid:17) ⊤ U (cid:18) α α (cid:19) U ⊤ (cid:16) x − x (2) (cid:17) , U = (cid:18) cos θ − sin θ sin θ cos θ (cid:19) , U = (cid:18) − cos θ sin θ sin θ cos θ (cid:19) , x (1 , = (0 , ± β ) ⊤ , (38)where x are not mass scaled. It has two minima, lo-cated at x (1 , , with normal modes given by matrices U , . The parameters were set to β = 2 , α = 1 . , α = 2 and m = 27 . Changing the angle θ changesthe angle between the normal modes of the two minima,as can be seen in Figure 1. With θ = 0 , the instantonpath is a straight line which connects the two minimaand, near minima, the path direction coincides with thelowest normal mode. As values of θ increase and nor-mal modes rotate, the instanton path does not rotate asquickly near minima. Instead, it picks up a non-zero dis-placement along the higher normal mode. It turns outthat this small displacement can significantly affect thesplitting. Combining Eqs. (21) and (26), we obtain F atthe dividing plane in the form F ( S cp ) = U ⊤ ( x ( ε ) − x (0))e ω e R S cp ε p d S ′ . (39) -3-2-10123-3 -2 -1 0 1 2 3-3-2-1012 -2 -1 0 1 2 3 00.511.522.533.5 x x x x V FIG. 1. Potential energy surfaces for model potential inEq. (38) ( α = 1 . , α = 2 , β = 2 ) for angles θ of, leftto right and top to bottom, 0, π/ , π/ and π/ . Super-posed on each potential energy surface are the correspondinginstanton pathways. The exponential growth of the F term in Eq. (39) is re-sponsible for this behavior. Even small displacementsalong the excited mode near minima can be magnifiedand result in an important contribution to the splitting.A useful parameter for quantifying the displacement nearminima is η = U ⊤ ( x ( ε ) − x (0)) /ε, (40)where the division with ε is made to cancel out the de-pendence on the step length ε , where it is observed. Thedependence of the displacement η on the angle θ is givenin Table I. It can be seen that the displacement is pre-dominantly along the lower mode for all angles θ in TableI. θ η (1 , η (0 , . . . . π/
12 0 . . . . π/ . . . . π/ . . . . TABLE I. Displacement η in Eq. (40) at ε = 0 . for the twonormal modes, (1,0) and (0,1), of the 2D symmetric potentialin Eq. (38). Fractional contribution of the F (L) F (R) term tothe tunneling splitting in Eq. (33), when the mode is excited,is given in parentheses. Table II shows the tunneling splittings in the groundstate and in the first two excited states, with the lower, (1 , , and the higher mode (0 , excited with one quan-tum of vibration. Convergence of the excited-state split-tings with the addition of F , U and Z terms in the exp( − w ) expansion is also shown. The exact quantum-mechanical results are obtained by the diagonalizationof Hamiltonian in the sine DVR basis with grid bound-aries at [ − . , . in both coordinates and 150 basis func-tions for each degree of freedom. They are given in TableII in parentheses for comparison. It can be seen that the Z term contribution is small for all the test cases. Thecontribution of F term is dominant for the longitudinalexcitation of the mode (1 , . On the other hand, whenthe higher mode (0 , is excited, the relative contributionof F and U terms changes with angle θ . Displacement η suggests that the excitation of (0 , is in the transversalmode. Indeed, at θ = 0 , F term does not contribute andthe U term determines the splitting, as in the theory ofRef. 53. But with an increase of θ , the F contributionquickly rises to account for more than of the split-ting at θ = π/ , while the displacement remains small at η = 0 . . This demonstrates that it is crucial to includethe F term in the expansion of exp( − w ) even when theexcited mode appears to be transversal. The contribu-tion from a small displacement can exponentially growand finally dominate the splitting. θ ∆ ∆ (1 ,
0) ∆ (0 , . −
8) 0 . . −
10) 1 . −
8) 5 . − . − . −
8) 5 . − . − . − π/
12 9 . −
9) 5 . − . −
10) 9 . −
9) 8 . − . − . −
9) 8 . − . − . − π/ . −
9) 4 . − . −
11) 1 . −
9) 4 . − . − . −
9) 4 . − . − . − π/ . −
11) 5 . − . −
12) 7 . −
11) 6 . − . − . −
11) 6 . − . − . − TABLE II. Tunneling splittings in the ground and first twovibrationally excited states for the potential in Eq. (38) atvarious angles θ obtained using instanton theory. The excited-state splittings are, top to bottom, obtained using the expan-sion of exp( − w ) to F , F + U i ∆ x i and F + U i ∆ x i + Z ij ∆ x i ∆ x j terms, respectively. The exact quantum-mechanical resultsare given in parentheses. The tunneling splittings are invariant with respect tothe position of the dividing plane when only F terms areconsidered, in accord with the analysis of Appendix D.The same is true for the splitting obtained with the inclu-sion of the U terms at θ = 0 . In this case, the instantonpath is a straight line and vector U remains perpendic-ular to the path. We can see that in Eq. (36), the lastterm disappears in that case, since the path curvature iszero. However, we observed in all other cases that thesplittings decrease as the position of the dividing planechanges from . S tot to . S tot . This decrease variesfrom . to . for the excitation in the lower, lon- gitudinal, mode and from to for the excitationin the higher, transversal, mode. This variation can beeliminated by using U ⊥ instead of U , in other words, byignoring the last term in Eq. (36). In this approach, the F term is still included, e.g., by using Eq. (39), while the U (L) ¯ A − U (R) contribution in Eq. (33) is computed with U ⊥ . This approach thus eliminates the dependence ofthe splitting on the position of the dividing plane, as dis-cussed in Appendix D. However, we noticed an increasein all computed splittings by as much as , which re-sulted in an overestimation of quantum-mechanical re-sults. Since the error introduced is larger than the varia-tion of splitting with the connection point position, usingthe full expression seems to be the preferable option.In Table III, we studied the dependence of splittingson the reduction of the mass of the system. Convergenceof the excited-state splittings with the addition of F , U and Z terms in the exp( − w ) expansion is again shown, aswell as the exact quantum-mechanical results in paren-theses. The reduction of mass causes an increase in theenergy of vibrational states, which provides an insightinto the limits of theory as the energy approaches thebarrier height. In the ground state, the effective barrierheight can be computed as V (0 , = V + 12 ( λ − ω − ω ) , (41)where V is the potential energy and λ is the nonega-tive eigenvalue of matrix A at the position of the barrier,whereas ω and ω are vibrational frequencies at the min-imum. If lower, longitudinal mode is excited, the effectivebarrier is lowered by ω and becomes V (1 , = V (0 , − ω , (42)while if the higher, transversal mode is excited, the effec-tive barrier changes as V (0 , = V (0 , − ω + λ . (43)As we reduce the effective barrier height, by varying themass in Table III, the instanton method starts to over-estimate the tunneling splittings. When V eff ≈ , theexcited-state splitting is overestimated by about a factorof 2, similarly to the earlier observations in the groundstate . This is mainly caused by the overestimation ofthe state energy in the harmonic approximation, whichis then used in the transport equation. Furthermore, asignificant effect comes from the underestimation of thenorm of the localized wavefunction in the harmonic ap-proximation, as it extends further on the other side ofthe barrier. Therefore, in the case of a ’shallow’ splittingor the ’over-the-barrier’ splitting, the estimates obtainedusing the instanton method should only serve as an upperlimit.0 m ∆ ∆ (1 , V (1 , ∆ (0 , V (0 , . . −
9) 5 . − . −
10) 9 . −
9) 8 . − . − . −
9) 1 .
273 8 . −
10) 1 . . − . − . . −
3) 2 . − . −
4) 4 . −
3) 4 . − . − . −
3) 0 .
731 4 . −
4) 1 . . − . − . .
214 1 . − . −
2) 0 .
215 3 . − . − .
219 0 .
051 3 . −
2) 0 . . . − . .
297 1 . − . −
2) 0 .
298 4 . − . − . − .
055 4 . −
2) 0 . . . − . .
740 0 . . −
2) 0 .
745 0 . . − . − .
445 0 .
141 0 . . . TABLE III. Tunneling splittings in the ground ( ∆ ) and first two vibrationally excited states ( ∆ ) for the potential in Eq. (38)at θ = π/ and various masses m obtained using instaton theory. The excited-state splittings are, top to bottom, obtainedusing the expansion of exp( − w ) to F , F + U i ∆ x i and F + U i ∆ x i + Z ij ∆ x i ∆ x j terms, respectively. The exact quantum-mechanical results are given in parentheses. For each excitation, the effective barrier heights V eff on the instanton path are alsogiven. B. ASYMMETRIC DOUBLE-WELL 2D POTENTIAL
We next perform tests on an asymmetric model 2D sys-tem. The potential profile along the MAP connecting anytwo minima does not have the left-right symmetry andthe maximum does not, in general, lie at the midpoint.The MAP can approach two minima along different nor-mal modes in an asymmetric system. The asymmetricpotential that we use in our tests is given by the follow-ing equations, V = 12 α ( x + β ) + 12 α ( x + β ) ,V = 12 α ( x − β ) + 12 α ( x + β ) ,V = 12 α ( x − β ) + 12 α ( x − β ) ,V = 12 α ( x + β ) + 12 α ( x − β ) ,V = V V V V V V V + V V V + V V V + V V V , (44)where x i are not mass scaled. The potential parametersin Eq. (44) are taken as β = 2 , α = 1 . α = 2 and m = 27 . The potential has four minima, and possesses a C symmetry axis, as shown in Figure 2. Instanton pathsconnect the neighboring minima as indicated in the fig-ure. The ’diagonal’ instanton paths have large actionsand are negligible. Energy levels split due to tunnelinginto a triplet, in which the middle level is doubly degen-erate. The tunneling splitting pattern consists of energy -4 -3 -2 -1 0 1 2 3 4-4-3-2-101234 012345 x x V FIG. 2. Potential energy surface for model potential inEq. (44) ( α = 1 . , α = 2 , β = 2 ) Superposed on potentialenergy surface are the instanton pathways that are responsi-ble for the formation of the tunneling splitting pattern. levels E = E − ∆ , E = E = E and E = E + ∆ ,where ∆ corresponds to the tunneling splitting betweenthe neighboring minima and E is the harmonic energy.We now label the minimum at ( − β, − β ) as ’left’ and theminimum at ( β, − β ) as ’right’. Each instanton path isalmost a straight line between two minima, however, be-cause of the anharmonicity, the path is slightly deflected1near minima. As a result of this deflection, it entersthe left minimum along the lower mode, instead of thehigher one, as explained in Appendix B. However, it alsopossesses a large displacement η in Eq. (40) along thehigher mode. The higher mode is therefore longitudinalat the left minimum, while the lower mode is longitudi-nal near the right minimum. As a result, when eitherof the modes is excited, it cannot be described as a lon-gitudinal or a transversal excitation with respect to theinstanton path. It represents the case of longitudinal-transversal excitation, where the excited mode is longi-tudinal at one minimum and trasversal to the path atthe other minimum. This case cannot be treated withthe method of Ref. 53. The localized wavefunction thatcorresponds to the longitudinal excitation is of the form p exp ( − / x ⊤ A ∆ x ) , which means that it is even inthe dividing plane. On the other hand, the wavefunc-tion that corresponds to the transversal excitation is ofthe form ( U ⊤ ∆ x ) exp ( − / x ⊤ A ∆ x ) , which is odd inthe dividing plane. As a result, the surface integral inHerring formula is odd and identically equal zero. It is -2-1.5-1-0.500.511.52-8 -6 -4 -2 0 2 4 6 8 Δ x Φ FIG. 3. Comparison of the wavefunctions in the dividingplane (line) obtained by using U ⊥ only (dotted line) and byusing F + U i ∆ x i (full line) in the preexponential factor of thelocalized wavefunction in Eq. (31). clear, however, from quantum-mechanical computationsthat the splitting is not zero, but can, in fact, even belarger than the splitting in the ground state, as can beseen in Table IV.In our treatment, the addition of F term breaks thesymmetry of the wavefunction in the dividing plane, andit moves the node away from the instanton trajectory,while the maximum of the Gaussian part in Eq. (31) stayson the trajectory, as shown in Figure 3. As a result, theintegral in Herring formula does not vanish. Results ob-tained using our approach are given in Table IV. Fromthe η values in the left minimum, it is clear that in itsvicinity, the instanton trajectory rapidly turns towards ∆ (1 ,
0) ∆ (0 , . −
11) 2 . − instanton . −
11) 3 . − . −
11) 3 . − QM . −
11) 6 . − η (L) . . η (R) . . TABLE IV. Tunneling splittings in first two vibrationally ex-cited states ( ∆ ) for the potential in Eq. (44) obtained usinginstanton theory. Displacements, η in Eq. (40), are given forthe left, ( − β, − β ) , and the right, ( β, − β ) , minimum. QM la-bels the exact quantum-mechanical results. The ground-statesplitting is ∆ = 9 . − , using the JFI method. Theexact result is ∆ = 8 . − . the direction of the second (higher) normal mode, whileit has to enter the minimum along the first (lower) mode.As a result of this sharp turn, F value for the left mini-mum is not zero and, in the end, gives rise to the non-zerotunneling splitting. Contribution of the Z term in bothexcited states is quite large compared to its contributionin the symmetric test case above. This is indicative of thepresence of non-negligible anharmonic effects in this sys-tem. The anharmonicity is also a probable reason for therelatively large discrepancies between the instanton andthe exact quantum-mechanical results (obtained on thesame grid as for the symmetric potential above), wherethe latter are and higher for the excitationof the first and second vibrational mode, respectively. Alarger discrepancy in the higher mode could be attributedto its larger energy, and the larger spread of its wavefunc-tion into the regions away from the instanton path whereanharmonicity is significant. C. WATER DIMER
The tunneling splitting pattern of water dimerhas been extensively studied both experimentaly andtheoretically , which makes it a good benchmarksystem to test our method. We chose the fully deuter-ated dimer over the non-deuterated one, because its vi-brational energies are lower. As a consequence, thereare more vibrational excitations which do not exceedthe barrier height, and can be treated with the instan-ton method. Analytical potential energy surface MB-pol was used in all calculations.Water dimer, shown labeled in Figure 4, has 8 equiva-lent symmetry-related and accessible minima, which cor-respond to the permutations of hydrogen and oxygenatoms that do not break the covalent H-O bonds. Per-mutations which do break the covalent bonds are con-sidered unfeasable. These minima are connected by fivedistinct tunneling rearrangement pathways . Accep-tor tunneling path (AT) corresponds to the permutation (34) . In the ground state, its effective barrier is relativelylow, V eff = 77 cm − ), so it gives rise to the largest tun-2 Mode AT GI AI BT DE1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . TABLE V. Left and right displacements η in Eq. (40) in deuterated water dimer for five instanton pathways and excitationsinto lowest five vibrational modes. Pathways are acceptor tunneling (AT), geared interchange (GI), antigeared interchange(AI), bifurcation tunneling (BT) and donor exchange (DE). AB4 3 1 2
FIG. 4. The minimum energy geometry of the water dimerlabeled to represent the reference version. neling matrix element. This matrix element is respon-sible for the splitting of energy levels into two groups,whose energy difference is called the acceptor splitting ∆(A) = 4 | h (AT) | . As seen in Table V, the displace-ments η for the AT path lie predominantly along thelowest mode at both minima. Next contribution to thesplitting pattern arises from the geared interchange (GI)and anti-geared interchange (AI) pathways, which cor-respond to the (AB)(1324) and (AB)(14)(23) permuta-tions. These pathways have larger effective barriers inthe ground state, V eff = 188 cm − and V eff = 227 cm − ,respectively. They cause the energy levels in both groups,formed by acceptor tunneling, to split into triplets, withthe energy width of the lower group called the lower inter-change ∆(LI) = 4 | h (GI) + h (AI) | , while the upper groupenergy width is called the upper interchange ∆(UI) =4 | h (GI) − h (AI) | . The AI path is mostly displaced alongthe lowest mode near minima as well, but has larger pro-jections onto the second and third mode. In contrast, theGI path is almost equally displaced along the first andsecond mode near minima, while it has to enter the min-ima along the lowest mode. Finally, the smallest contri-bution to the splitting pattern of water dimer arises fromthe bifurcation tunneling (BT) and donor exchange (DE)paths, which correspond to the (12)(34) and (12) permu- tations, respectively. These pathways possess the highesteffective barriers, V eff = 469 cm − and V eff = 581 cm − ,respectively. They cause the shifts in the energies of thetriplets by the amounts called the lower bifurcation andthe upper bifurcation, ∆(LB) = | h (BT) + 4 h (DE) | and ∆(UB) = | h (BT) − h (DE) | . Bifurcation tunneling pathis displaced mostly along the third mode near minima.Donor exchange path, on the other hand, is displacedmostly along the lowest mode near one minimum, whileit is displaced mostly along the second and third modenear the other minimum. Therefore, this path representsa realistic case of the asymmetric potential which featureslongitudinal-transversal excitations that we discussed inthe previous subsection on a 2D model potential.The lowest mode of vibration in the deuterated wa-ter dimer corresponds to donor torsion and has a fre-quency of ω = 84 cm − . In order to calculate the split-ting pattern with the excited donor torsion, we calculatethe matrix elements, h = − ∆ / , for all five rearrange-ment paths. The AT matrix element, obtained by theinstanton method, is 3 times larger than the experimen-tal value, as seen in Table VI. Since donor torsion is thelongitudinal mode of the AT path and its excitation fre-quency is larger than the effective barrier on the path,this represents a case of over-the-barrier tunneling. Theinstanton method is known to overestimate the splittingsby a factor of 2 − , as alsonoted in the previous subsection. The sign of the ac-ceptor splitting is found to be opposite to that of theground state, indicating that the groups of states associ-ated with the lower and upper interchange change places.This observation is in agreement with the experimen-tal measurements and the exact quantum-mechanicalcalculations .GI and AI matrix elements are found to be in goodagreement with the experimental results in their abso-lute values, but their relative sign appears to be wrong.This results in the wrong ordering of the LI and UI split-tings in magnitude, as seen in Table VII. We note thatthe contribution of the F term accounts for and of the matrix element in Eq. (33). A large contributionfor the AI path is expected, as donor torsion is its longi-3tudinal mode. However, for the GI path, which lies alonga combination of modes near minima, the contribution of F term is also important. We presume that the disagree-ment between the instanton and quantum-mechanical re-sults of Ref. 25 is caused by a large rotation-vibrationcoupling in the excited mode, which mixes the vibra-tional states of K a = 0 and K a = 1 and is not accountedfor in the instanton method. The values obtained forLI and UI ( . cm − and . cm − ) are, in fact, ina better agreement with the experimental values for K a = 1 , which are . cm − and . cm − , both inmagnitude and in ordering.Lower and upper bifurcations are underestimated forthe first excited vibrational mode, as can be seen in Ta-ble VII. For the DE path, this represents a longitudinal-transversal excitation, and it was shown for the modelpotential above that an underestimate is expected be-cause of the unaccounted anharmonicities. However, thedifference between the lower and upper bifurcation is notzero, as it would be in using the theory of Ref. 53, andeven though it is underestimated, a rough estimate ofits value is obtained. The exact quantum-mechanicalcalculations do not report it, probably due to the dif-ficulty in converging the values with sufficient accuracy.It is also worth mentioning that the UB and LB changesignificantly in the K a = 1 rotational state, to . − cm − for UB and . − cm − for LB. These valuesare again in better agreement with those that we com-puted, as in the case of the AT path, which provides fur-ther indication that the coupling of the first excited stateto rotations plays a significant role. Finally, the UB andLB are underestimated even in the ground vibrationalstate, which suggests the possibility that the BT and DEpathways are poorly described by the PES, either by toolarge potential energy barriers, or by slightly misplacedinstanton paths, both of which can have a drastic effecton the splittings.The second mode corresponds to the acceptor twist,with frequency ω = 100 cm − , while the third mode cor-responds to the acceptor wag, with frequency ω = 110 cm − . However, in quantum-mechanical calculations ,the order of these two motions changes, and the acceptorwag frequency drops to cm − , while the acceptor twistdrops to cm − . The large deviation of vibrational en-ergies from the harmonic frequencies is a strong indica-tion of large anharmonic effects in these two vibrationalmodes. Furthermore, since their energy difference is verysmall, it was noticed that these states interact througha Coriolis perturbation adding to the quantitative dis-agreement with the harmonic analysis. Nevertheless, thesplittings obtained from the second excited mode are ingood agreement with the experimental results. We notethat the F term on the AI path contributes with around to the matrix element, even though the displace-ments near minima along this mode are small. The over-estimation of the GI matrix element can be explainedby the fact that the path has a large projection ontothe second mode near minima, which means that the ef- fective barrier is significantly lowered. Discrepancy ofthe AI matrix element can be explained by the inaccu-racy of the PES, since quantum-mechanical results ona similar surface also overestimate this matrix element.Upper and lower bifurcations are again underestimated,probably for the same reasons as above, namely the in-adequate PES and the unaccounted anharmonic effectsin the longitudinal-transversal excitation.In the case of the third mode excitation, especially in-teresting is the AT path for which the contributions ofthe F term and the U term in the matrix element al-most cancel each other out, while the major contributionarises from the anharmonicity contained in the Z term.For this excitation, both GI and AI matrix elements areoverestimated. This can again be attributed to the rovi-brational coupling, since the quantum-mechanical resultsshow a significant increase in the lower and upper inter-change with the excitation to K a = 1 rotational state .Upper and lower bifurcations for this excitation show amuch better agreement with the experimental values than above.At larger excitation frequencies, the theory breaksdown. A probable cause of this breakdown is the factthat as the frequency increases, the contribution of the w term to the overall splitting rises significantly. This isdue to the fact that the F contribution depends exponen-tially on the frequency of excitation, while the η values donot compensate it. As a result, its contribution becomescomparable to that of W , while the WKB approach as-sumes ln F << W . A good test of the reliability of theobtained results is to redo the calculations with a differ-ent value of the initial ’jump’ parameter ε . As the valueof ε is reduced, the results should converge to the correctvalue. However, there is a limit to how much ε can be re-duced, as the propagation from the point too close to theminimum is not stable . If the results converge beforethis breakdown, they can be treated as reliable. Also, asthe value of ε is increased, values of the splittings shouldnot change by more than a few percent. This is the casefor the excitations in the first three lowest modes. For thefourth excited mode, if we change ε from 0.1 m / a to 1 m / a , the AT matrix element changes from .
98 cm − to − , which is an indication that the break-down of theory occured. Similar behaviour is presentfor the AI pathway, where the matrix element changesfrom . −
2) cm − to .
23 cm − . The change is not asdrastic as in the AT case, but it indicates that the errorbars on our results are very large, which also explainsthe discrepancies of results for the LI and UI splittings.Noticeable changes are also present for the DE pathway(from − . −
4) cm − to − . −
4) cm − ), while thevalues for other pathways do not change appreciably andcan be considered reliable.4 Mode AT GI AI BT DEGS .
766 9 . −
3) 4 . −
4) 1 . −
4) 3 . − − . − . −
2) 1 . −
2) 1 . −
9) 2 . − − . − . −
2) 1 . − − . −
5) 1 . − − . − . −
2) 1 . − − . −
5) 8 . − . . − . − − ) ( − )(3 .
92) (6 . − . − − ) ( − ) − . − .
256 4 . − − . − − . − − . − .
254 4 . − − . − − . − − . − .
261 5 . − − . − − . − . . . − − ) ( − )(0 . . . − − ) ( − ) .
15 0 .
147 2 . −
2) 5 . −
3) 3 . − . −
2) 0 .
141 2 . −
2) 5 . −
3) 2 . − .
469 0 .
143 2 . −
2) 5 . −
3) 5 . − . . − . − − ) ( − )(0 .
45) (2 . − . − − ) ( − ) . . −
2) 2 . −
2) 2 . − − . − .
98 3 . −
2) 2 . −
2) 2 . − − . − − . . − . −
2) 6 . −
2) 4 . − − ) ( − ) ( − ) ( − ) ( − )(1 .
23) (0 . . − − ) ( − ) TABLE VI. Tunneling matrix elements − h/ cm − for different tunneling pathways in deuterated water dimer (D O) obtainedusing instanton theory. Pathways described are acceptor tunneling (AT), geared interchange (GI), antigeared interchange(AI), bifurcation tunneling (BT) and donor exchange (DE). The excited-state splittings are, top to bottom, obtained using theexpansion of exp( − w ) to F , F + U i ∆ x i and F + U i ∆ x i + Z ij ∆ x i ∆ x j terms, respectively. The splittings given in parentheses areexperimental (top) and quantum-mechanical (bottom) results. (Ground-state (GS) experimental results are from Refs. 25and 78.) Mode A UI LI UB LBGS .
06 3 . −
2) 4 . −
2) 1 . −
4) 1 . − .
77) (3 . − . − . − . − .
78) (3 . − . − − ) ( − ) . .
257 0 .
109 1 . −
5) 1 . − . .
276 0 .
129 8 . −
5) 5 . − . .
290 0 .
134 7 . −
5) 6 . − . . . . − . − .
68) (0 .
20) (0 .
33) ( − ) ( − ) .
01 1 .
04 1 .
01 2 . −
4) 2 . − .
04 1 .
04 1 .
00 5 . −
5) 4 . − .
83 1 .
07 1 .
02 3 . −
5) 4 . − . . . . − . − .
03) (0 .
73) (0 .
39) ( − ) ( − ) .
60 0 .
503 0 .
673 5 . −
3) 5 . − .
109 0 .
479 0 .
649 5 . −
3) 5 . − .
88 0 .
484 0 .
660 5 . −
3) 5 . − . . . . − . − .
81) (0 .
11) (0 .
12) ( − ) ( − ) . .
012 0 .
22 3 . −
3) 1 . − . .
025 0 .
231 3 . −
3) 9 . − .
04 0 .
725 6 . −
2) 6 . − − ) ( − ) ( − ) ( − ) ( − )(4 .
9) (0 .
38) (1 .
0) ( − ) ( − ) TABLE VII. Acceptor, upper interchange and lower interchange splittings ( cm − ) in deuterated water dimer (D O) obtainedusing instanton method. The excited-state splittings are, top to bottom, obtained using the expansion of exp( − w ) to F , F + U i ∆ x i and F + U i ∆ x i + Z ij ∆ x i ∆ x j terms, respectively. The splittings given in parentheses are experimental (top) andquantum-mechanical (bottom) results. (Ground-state (GS) experimental results are from Refs. 25 and 78.) V. CONCLUSIONS
We developed a semiclassical theory for calculatingtunneling splittings of low-lying vibrationally excitedstates based on the instanton method. A WKB wave-function is constructed along the instanton path and itsharmonic neighborhood for each well, and inserted intoHerring formula to obtain the splitting that matches theJFI result in the ground state . The excited-state split-tings are then obtained constructing excited-state wave-functions analogously. The procedure closely follows thatof Ref. 53, but uses a more general boundary conditionnear minima and does not assume the left-right mirrorsymmetry of potential along the instanton path. In ourapproach, transversal and longitudinal excitations do notrequire separate treatments as in Ref. 53. This allows usto compute splittings in the systems where the excited vi-brational mode does not line up along the instanton pathnear minima, but has both longitudinal and transversalcomponents, or the systems in which the excited modeis longitudinal at one minimum and transversal at theother. Both components are propagated simultaneouslyalong the instanton path and cross interaction is kept inthe treatment.The tests on the symmetric double-well model poten-tial showed that a high accuracy can be expected forlow-lying states below the barrier. It was shown thatfor transversal modes, even a small longitudinal displace-ment near minima can dominate the tunneling splitting.We also observed that the longitudinal-transversal crossterms improve results. The tests on the asymmetricmodel potential showed that we can calculate splittingestimates for excited longitudinal-transversal modes, al-beit with somewhat reduced accuracy. Finally, we cal-culated the tunneling splitting pattern of the deuteratedwater dimer in vibrationally excited lowest three modesby computing contributions from five different rearrange-ment pathways. This is a particularly challenging sys-tem for treatment with partly harmonic theories. Addi-tionaly, the system exhibits significant rovibrational cou-plings, which are, at present, neglected in our treatment.We could nevertheless obtain reasonable agreement inmany cases in a system which showcases the situationsin which the present theory gives significantly differentresults from that of Ref. 53.Tunneling splittings in vibrationally excited states re-quire no additional information about the molecular sys-tem. All computational effort is concentrated, as for theground-state splittings, in determining the MAP by opti-mization and the evaluation of Hessians along the MAP.This allows us to compute and interpret splitting patternsin many mid-sized molecules using state-of-the-art po-tentials. The theory is applied in Cartesian coordinatesand requires no modification for treating different molec-ular systems. However, tunneling splittings in vibrationalstates with higher frequencies, such as the excitations oflibrational modes of water trimer and pentamer thatwere recently measured, cannot be treated with the the- ory in the present format. Also, many small tunnelingsystems exhibit large rotation-vibration coupling, whichis currently neglected and can affect the splittings. Acomputationally tractable theory for calculating split-tings in rotationally excited states would also be desir-able. These are some of the immediate challenges remain-ing in which the future efforts will certainly be directedin a quest to provide quantitative estimates for splittingpatterns for molecules and clusters that are out of reachto the exact quantum-mechanical treatments. ACKNOWLEDGMENTS
This work was supported by Croatian Science Foun-dation Grant No. IP-2016-06-1142, and in part by theQuantiXLie Centre of Excellence, a project cofinanced bythe Croatian Government and European Union throughthe European Regional Development Fund – the Compet-itiveness and Cohesion Operational Programme (GrantKK.01.1.1.01.0004).
DATA AVAILABILITY
The data that support the findings of this study areavailable from the corresponding author upon reasonablerequest.
Appendix A: Method of Characteristics and local coordinates
Method of characteristics is a technique for solving par-tial differential equations . It relies on locating curves,the characteristics, along which the gradient of the de-sired solution is tangential. As a consequence, the partialdifferential equation reduces to an ordinary differentialequation. For a non-linear partial differential equation ofthe form, F ( x , ..., x N , p , ..., p N , f ) = 0 , (A1)where p i = ∂f /∂x i , defining equations of the character-istics are d x i d τ = ∂F∂p i , d p i d τ = − ∂F∂x i − ∂F∂f p i , d f d τ = ∂F∂p i p i , (A2)where τ parametrizes the characteristic.Hamilton-Jacobi equation is a non-linear partial dif-ferential equation for which F = p i p i − V , where p i = ∂W /∂x i . Its characteristics are therefore d x i d τ = p i , d p i d τ = ∂V∂x i . (A3)6The characteristics describe classical trajectories on theinverted PES, while ∇ W is the momentum on the tra-jectory. The total energy of the classical motion is E tot = p i p i + ( − V ) = 0 . On characteristics, W isfound by solving d W d τ = p i p i = 2 V. (A4)The parameter τ represents time and, as the trajectoryapproaches minimum, its value τ → −∞ . This is nu-merically problematic, so we reparametrize characteris-tics with the arc length distance from the minimum, S ,using the transformation in Eq. (7).In order to expand W in Taylor series around the char-acteristic, it is convenient to define a set of local coordi-nates { S, ∆ x } . Since coordinate S parametrizes charac-teristic, it is only defined for the points lying on it. Inorder to assign a value S to the point that does not lieon the characteristic, a point x ( S ) which does lie on itis chosen so that ( x i − x i ( S )) p i = 0 , (A5)that is, x ( S ) is chosen so that the vector connectingit with the point x is orthogonal to the characteristic at x ( S ) . The value of S which corresponds to x ( S ) is thenassigned to x . The orthogonal coordinates ∆ x are thendefined as ∆ x = x − x ( S ) . Differentiation of Eq. (A5)gives ∂S∂x i = p i p − a ⊤ ∆ x p = p i p (cid:18) a ⊤ ∆ x p + ... (cid:19) , (A6)where a = d p d τ denotes the acceleration. From the differ-entiation of Hamilton-Jacobi equation, Eq. (3), we obtain a = Ap . And, finally, the differentiation of the definingequation of orthogonal coordinates in Eq. (A5) gives thetransformation ∂ ∆ x i ∂x j = δ ij − p i p j p (cid:18) a ⊤ ∆ x p + ... (cid:19) . (A7)Eqs. (A6) and (A7) are used throughout the paper totransform between Cartesian and local coordinates onthe characteristic as A , U , B and Z are all given indifferential form. Appendix B: Wavefunctions near minima
Near minima x min , the PES can be approximated bya harmonic oscillator potential V = 12 ω i q i , (B1)where q i = V ji ( x j − x min j ) are normal coordinates, and ω i corresponding harmonic frequencies. Since A = H / , we have V ⊤ A V = Ω , with ( Ω ) ij = ω i δ ij . In the harmonic region near minima, the equations of char-acteristics, Eq. (A3), become d q i d τ = ∂V∂q i , d q i d τ = ω i q i . (B2)The trajectory along the characteristic from the mini-mum to an arbitrary point q at τ = 0 inside the har-monic region is q i ( τ ) = q i e ω i τ . (B3)By considering the tangent vector of the characteristic, t i = p i p = ω i q i e ω i τ q ω j q j e ω j τ , (B4)we note that in the limit τ → −∞ , the tangent becomes t i = δ iM , where M denotes the lowest frequency normalmode for which q ,M = 0 . This means that all charac-teristics approach the minimum along the lowest normalmode with a non-zero projection upon entering the har-monic region.Function W in Eq. (A4) can be evaluated in the har-monic region at the characteristic as W ( τ ) = Z τ −∞ ω j q j e ω j τ ′ d τ ′ , (B5)or, making use of Eq. (B3), as W ( q ) = 12 ω j q j . (B6)Furthermore, since in the harmonic region A ≈ A ,the ground-state wavefunction corresponds to that of theharmonic oscillator, φ = e − ω j q j . (B7)Eq. (B7) is used to approximate the norm of the ground-state wavefunction in Herring formula Eq. (1).For vibrationally excitated states, the correct form ofthe wavefunction at the minimum is obtained by choosing ( U ) i = V i e , that is by equating the vector U with theexcited normal mode at the minimum. The wavefunctionthen has the form φ = q e e − ω j q j . (B8)For a point on the characteristic, which lies in the har-monic region, ∆ x i = 0 , so its form is φ = F ( ε )e ω j q j ( ε ) . (B9)Therefore, the initial condition for the F term at S = ε has to be F ( ε ) = q e = U ⊤ ( x ( ε ) − x (0)) , (B10)in order to yield the correct form of the wavefunction inEq. (B8).7 Appendix C: Anharmonicity about the instanton path
The anharmonicity of potential in the directions per-pendicular to the instanton path can be partially ac-counted for by including the higher derivatives of thePES along the instanton path, beyond Hessian, in thesemiclassical treatment of Section III. We assume belowthat the third derivative tensor of the PES, with elements c ijk = ∂ V∂x i ∂x j ∂x k along the instanton path has been deter-mined. This allows us to compute the third derivativesof function W , B ijk = ∂ W ∂x i ∂x j ∂x k , in Taylor expan-sion Eq. (10). The equation for propagation of tensor B is obtained by differentiating Hamilton-Jacobi equation,Eq. (3), three times as, ∂ W ∂x i ∂x j ∂x k ∂x l ∂W ∂x l + ∂ W ∂x i ∂x j ∂x l ∂ W ∂x l ∂x k + ∂ W ∂x i ∂x l ∂x k ∂ W ∂x l ∂x j + ∂ W ∂x l ∂x j ∂x k ∂ W ∂x l ∂x i = ∂ V∂x i ∂x j ∂x k . (C1)The first term in Eq. (C1) represents a directional deriva-tive of the tensor element B ijk along the instanton tra-jectory, while the other terms can be recognized as tensorelements of B and of Hessian A , which is determined bysolving Eq. (9). Eq. (C1) on the instanton reads p B ′ ijk + B ijl A lk + B ilk A lj + B ljk A li = c ijk . (C2)We proceed to determine the initial condition B ( ε ) inthe vicinity of the minimum. For that purpose we lin-earize Eq. (C2), following an analogous procedure to thatfor A in Refs. 52 and 66, as B = B (0) + B (1) S, c = c (0) + c (1) S, A = A (0) + A (1) S,p = p (1)0 S. (C3)Inserting the above expressions into Eq. (C2) and equat-ing terms of the same order in S yields equations for B (0) and B (1) as B (0) ijl A (0) lk + B (0) ilk A (0) lj + B (0) ljk A (0) li = c (0) ijk ,p (1)0 B (1) ijk + B (1) ijl A (0) lk + B (1) ilk A (0) lj + B (1) ljk A (0) li == c (1) ijk − B (0) ijl A (1) lk − B (0) ilk A (1) lj − B (0) ljk A (1) li . (C4)These are solved by transforming to the basis of normalmodes, the eigenvectors of A (0) , using the following rela-tions, ω i δ ij = V i ′ i V j ′ j A (0) i ′ j ′ , ˜ B (0) ijk = V i ′ i V j ′ j V k ′ k B (0) i ′ j ′ k ′ , ˜ c (0) ijk = V i ′ i V j ′ j V k ′ k c (0) i ′ j ′ k ′ . (C5) Inserting Eq. (C5) into Eq. (C4) yields equations ˜ B (0) ijk = ˜ c (0) ijk ω i + ω j + ω k , (C6) ˜ B (1) ijk = ˜ c (1) ijk − ˜ B (0) ijl ˜ A (1) lk − ˜ B (0) ilk ˜ A (1) lj − ˜ B (0) ljk ˜ A (1) li p (1)0 + ω i + ω j + ω k , (C7)that are needed to construct B ( ε ) . Eq. (C2) can now besolved in the interval [ ε, S ] using any differential equationsolver, such as the Runge Kutta method .Tensor B cannot be included in the wavefunction ofEq. (2) without the inclusion of fourth derivatives, asthe resulting wavefunction would not be integrable in thedividing plane. However, it is used below to compute the Z term in the expansion of exp( − w ) , Eq. (34), and thusindirectly account for a part of anharmonicity.We first note the following expressions are valid on theinstanton path, F = e − w ,U i = ∂∂x i e − w = − ∂w∂x i e − w ,Z ij = ∂ ∂x i ∂x j e − w = − ∂ w∂x i ∂x j e − w + ∂w∂x i ∂w∂x j e − w ,Z ij p j = (cid:18) ∂∂x j U i (cid:19) p j = p U ′ i . (C8)In the next step, we differentiate Eq. (19) twice to obtainuseful relations p k ∂ w∂x i ∂x k e − w = A ik U k ,p k ∂ w∂x i ∂x j ∂x k e − w = B ijk U k + A ik Z kj + A jk Z ki + ∂w∂x j A ik U k + ∂w∂x i A jk U k . (C9)Finally, we take the third derivative of exp( − w ) inEq. (34) to arrive at (cid:18) ∂∂x k Z ij (cid:19) p k = − p k ∂ w∂x i ∂x j ∂x k e − w + p k ∂w∂x k ∂ w∂x i ∂x j e − w + p k ∂w∂x i ∂ w∂x k ∂x j e − w + p k ∂w∂x j ∂ w∂x i ∂x k e − w − p k ∂w∂x k ∂w∂x i ∂w∂x j e − w , (C10)where we insert Eq. (C9) and recognize Eq. (C8) to ob-tain the equation for Z in the following form, p Z ′ ij + A ik Z kj + A jk Z ki + B ijk U k + ω e Z ij = 0 . (C11)This equation is again solved separately in the interval [0 , ε ] and [ ε, S ] , following the same procedure as for A and B . All objects are expanded up to linear terms in8 S and inserted into Eq. (C11). By equating terms of thesame order in S , we obtain equations for Z (0) and Z (1) , A (0) Z (0) + Z (0) A (0) + ω e Z (0) + B (0) U (0) = 0 ,p (1)0 Z (1) + A (0) Z (1) + Z (1) A (0) + ω e Z (1) + A (1) Z (0) + Z (0) A (1) + B (1) U (0) + B (0) U (1) = 0 , (C12)where matrices ( BU ) ij evaluate as B ijk U k . These equa-tions are again solved by transforming to the basis ofnormal modes, the eigenvectors of A (0) , as ˜ Z (0) ij = − ( ˜ B (0) ˜ U (0) ) ij ω e + ω i + ω j , ˜ Z (1) ij = − ( ˜ B (1) ˜ U (0) ) ij + ( ˜ B (0) ˜ U (1) ) ij + ( ˜ A (1) ˜ Z (0) ) ij + ( ˜ Z (0) ˜ A (1) ) ij p (1)0 + ω e + ω i + ω j . (C13)We are now in the position to compute Z ( ε ) , which servesas the initial condition for the propagation of Z in theinterval [ ε, S ] by solving Eq. (C11) using, e.g., Runge-Kutta method.When the Z term is included in the expansion of exp( − w ) , the tunneling splitting formula assumes the fol-lowing form ∆ =∆ (2 ω e ) (cid:18) F (L) F (R) + 12 U (L) ⊤ ¯ A − U (R) + 14 F (L) Tr (cid:16) Z (R) ¯ A − (cid:17) + 14 F (R) Tr (cid:16) Z (L) ¯ A − (cid:17) (cid:19) , (C14)where terms of the form Z ij Z kl ∆ x i ∆ x j ∆ x k ∆ x l in thesurface integral have been neglected, as their contributionwas found to be negligible. Appendix D: Invariance of tunneling splittings with respect tothe position of the dividing plane
Invariance of the ground-state tunneling splitting for-mula can be proved by differentiating Eq. (17) with re-spect to the position of the connection point S cp , wherethe dividing plane intersects the instanton path, ∂ ∆ ∂S cp = ∆ p (cid:18) ∂p ∂S cp − p det ′ ¯ A ∂∂S cp det ′ ¯ A − p ∂∂S cp W (L)1 − p ∂∂S cp W (R)1 (cid:19) . (D1)The derivative of determinant in Eq. (D1) is simplifiedusing Jacobi formula (Eq. (C4) in Ref. 66), while W L / R1 functions are differentiated in the upper/lower limit of the integral in Eq. (11), ∂ ∆ ∂S cp = ∆ p (cid:18) ∂p ∂S cp − Tr (cid:18) ¯ A − p ∂∂S cp ¯ A (cid:19) − Tr( A (L) − A ) + Tr( A (R) − A ) (cid:19) . (D2)Derivative of ¯ A can be shown to equal ∂∂S cp ¯ A = 12 ¯ A (cid:16) A (R) − A (L) (cid:17) + 12 (cid:16) A (R) − A (L) (cid:17) ¯ A , (D3)where use has been made of Eqs. (16) and (9). Further-more, since the tangent t is an eigenvector of ¯ A withzero eigenvalue and, by definition of the pseudoinverse, ¯ A − t = 0 , we have P ¯ AP = ¯ A and P ¯ A − P = ¯ A − ,where P = I − tt ⊤ is the operator that projects out thetangent of the instanton path. Using the above, one canshow that Tr (cid:18) ¯ A − p ∂∂S cp ¯ A (cid:19) = Tr (cid:16) P (cid:16) A (R) − A (L) (cid:17) P (cid:17) = Tr (cid:16) A (R) ⊥ − A (L) ⊥ (cid:17) . (D4)Thus, the derivative of the tunneling splitting becomes ∂ ∆ ∂S cp = ∆ p (cid:18) ∂p ∂S cp + Tr( A (L) ⊥ − A (R) ⊥ ) − Tr( A (L) − A (R) ) (cid:19) . (D5)Finally, since Tr A (L) = p ′ + Tr A (L) ⊥ and Tr A (R) = − p ′ +Tr A (R) ⊥ , as shown in Ref. 52, we have ∂ ∆ ∂S cp = ∆ p (cid:18) ∂p ∂S cp + Tr( A (L) ⊥ − A (R) ⊥ ) − Tr( A (L) ⊥ − A (R) ⊥ ) − ∂p ∂S cp (cid:19) = 0 , (D6)which proves that the ground-state tunneling splittingdoes not depend on S cp , the position of the dividingplane.Similarly, the invariance of the excited-state tunnelingsplitting on the position of the dividing plane is checkedby differentiating Eq. (33) with respect to S cp . If onlythe F terms are included in the expansion of exp( − w ) ,we have ∂ ∆ ∂S cp = ∆ ω e (cid:18) ∂F (L) ∂S cp F (R) + F (L) ∂F (R) ∂S cp (cid:19) , (D7)which together with Eq. (21) gives ∂ ∆ ∂S cp = ∆ ω e (cid:18) ω e p F (L) F (R) − ω e p F (L) F (R) (cid:19) = 0 . (D8)9If we include the U terms in the expansion of exp( − w ) ,the derivative of the splitting becomes ∂ ∆ ∂S cp =∆ ω e (cid:18) ∂∂S cp U (L) ⊤ ¯ A − U (R) + U (L) ⊤ ∂∂S cp ¯ A − U (R) + U (L) ⊤ ¯ A − ∂∂S cp U (R) (cid:19) . (D9)It can be shown that p ∂∂S cp ¯ A − = PA (L) ¯ A − − ¯ A − A (R) P , (D10)which can be used to rewrite Eq. (D9) as ∂ ∆ ∂S cp =∆ ω e p (cid:18) − ( U (L) − PU (L) ) ⊤ A (L) ¯ A − U (R) + U (L) ⊤ ¯ A − A (R) ( U (R) − PU (R) ) (cid:19) . (D11)In this form, it is evident that if U remains orthogonalto the instanton path, i.e., PU = U , the excited-statetunneling splittings become independent on the positionof the dividing plane. If that is not the case, however,Eq. (D9) can be further simplified to ∂ ∆ ∂S cp =∆ ω e p (cid:18) − ω e p F (L) t ⊤ A (L) ¯ A − U (R) − ω e p F (R) U (L) ⊤ ¯ A − A (R) t (cid:19) , (D12)which does not vanish and the tunneling splitting will, ingeneral, depend on the position of the dividing plane, asobserved in Section IV.If the same analysis is performed with the Z terms,there arise two factors which cancel out the U terms.However, a multitude of other factors also arise, whichagain cause the dependence on the position of the divid-ing plane. As mentioned above, the root of the problemis that the expansion of exp( − w ) is inconsistent with theexpansion of W , and it gives rise to terms of all ordersin ∆ x . However, in the case of a symmetric potential,all perpendicular components of the gradients, Hessiansand third-order tensors are the same for the left- andright-localized wavefunctions at the dividing plane in themiddle of the instanton path, while their tangent com-ponents differ in sign. Thus, it is possible to show thatthe derivative of the Z contribution with respect to S cp vanishes at the middle of the instanton path, and nu-merical tests show that its contribution is minimal there.Therefore, for symmetric systems, the middle of the pathrepresents the optimal position of the dividing plane. Forthe asymmetric paths, there is no such preferential pointon the instanton. However, good results are obtained bypositioning the dividing plane at the maximum of thebarrier, as at this point p is the largest, and the deriva-tives of the splitting are generally smallest, which meansthat, at this point, the splittings are relatively stable. R. P. Bell,
The Tunnel Effect in Chemistry (Chapman and Hall,London, 1980). L. H. Coudert and J. T. Hougen,J. Mol. Spectrosc. , 86 (1988). T. R. Walsh and D. J. Wales, J. Chem. Soc. Faraday Trans. ,2505 (1996). Y. Xu and W. Jäger, J. Chem. Phys. , 7968 (1997). F. N. Keutsch and R. J. Saykally,P. Natl. Acad. Sci. USA , 10533 (2001). K. Liu, M. G. Brown, J. D. Cruzan, and R. J. Saykally, Science , 62 (1996). M. T. Cvitaš and J. O. Richardson, in
Molecular Spectroscopyand Quantum Dynamics , edited by R. Marquardt and M. Quack(Elsevier, 2020) Chap. 10. T. Hammer, M. D. Coutinho-Neto, A. Viel, and U. Manthe,J. Chem. Phys. , 224109 (2009). C. Fábri, R. Marquardt, A. G. Császár, and M. Quack,J. Chem. Phys. , 014102 (2019). P. M. Felker and Z. Bačić, J. Chem. Phys. , 024305 (2019). M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie,A. Michaelides, M. A. Morales, and T. E. Markland,Chem. Rev. , 7529 (2016). Y. Wang, X. Huang, B. C. Shepler, B. J. Braams, and J. M.Bowman, J. Chem. Phys. , 094509 (2011). S. K. Reddy, S. C. Straight, P. Bajaj, C. Huy Pham, M. Riera,D. R. Moberg, M. A. Morales, C. Knight, A. W. Götz, andF. Paesani, J. Chem. Phys. , 194504 (2016). J. O. Richardson, S. C. Althorpe, and D. J. Wales,J. Chem. Phys. , 124109 (2011). F. N. Keutsch, R. S. Fellers, M. R. Viant, and R. J. Saykally,J. Chem. Phys. , 4005 (2001). W. T. S. Cole, R. S. Fellers, M. R. Viant, and R. J. Saykally,J. Chem. Phys. , 014306 (2017). M. T. Cvitaš and J. O. Richardson,Phys. Chem. Chem. Phys. , 1035 (2019). B. J. Mhin, J. Kim, S. Lee, J. Y. Lee, and K. S. Kim, J. Chem.Phys. , 4484 (1994). C. Pérez, M. T. Muckle, D. P. Zaleski, N. A. Seifert, B. Temelso,G. C. Shields, Z. Kisiel, and B. H. Pate, Science , 897 (2012). J. O. Richardson, C. Pérez, S. Lobsiger, A. A. Reid, B. Temelso,G. C. Shields, Z. Kisiel, D. J. Wales, B. H. Pate, and S. C.Althorpe, Science , 1310 (2016). C. Léonard, N. C. Handy, S. Carter, and J. M. Bowman,Spectrochimica Acta Part A , 825 (2002). M. Neff and G. Rauhut, Spectrochimica Acta Part A , 100 (2014). J. Šmydke, C. Fábri, J. Sarka, and A. G. Császár,Phys. Chem. Chem. Phys. , 3453 (2019). F. Wu, Y. Ren, and W. Bian,J. Chem. Phys. , 074309 (2016). C. Leforestier, K. Szalewicz, and A. van der Avoird,J. Chem. Phys. , 014305 (2012). X.-G. Wang and T. Carrington,J. Chem. Phys. , 074108 (2018). M. Schröder, F. Gatti, and H.-D. Meyer,J. Chem. Phys. , 234307 (2011). M. Schröder and H.-D. Meyer,J. Chem. Phys. , 034116 (2014). T. Hammer and U. Manthe, J. Chem. Phys. , 224305 (2011). D. Blume and K. B. Whaley, J. Chem. Phys. , 2218 (2000). A. Viel, M. D. Coutinho-Neto, and U. Manthe,J. Chem. Phys. , 024308 (2007). Y. Wang, B. J. Braams, J. M. Bowman, S. Carter, and D. P.Tew, J. Chem. Phys. , 224314 (2008). C. L. Vaillant, D. J. Wales, and S. C. Althorpe,J. Phys. Chem. Lett. , 7300 (2019). D. J. Nesbitt and F. Dong,Phys. Chem. Chem. Phys. , 2113 (2008). C. Qu and J. M. Bowman,Phys. Chem. Chem. Phys. , 24835 (2016). S. C. Althorpe and D. C. Clary, J. Chem. Phys. , 4390 (1995). I. Matanović, N. Došlić, and B. R. Johnson,J. Chem. Phys. , 084103 (2008). E. Kamarchik, Y. Wang, and J. Bowman, J. Phys. Chem. A , 7556 (2009). T. D. Sewell, Y. Guo, and D. L. Thompson,J. Chem. Phys. , 8557 (1995). C. S. Tautermann, A. F. Voegele, T. Loerting, and K. R. Liedl,J. Chem. Phys. , 1967 (2002). M. Ceotto, Mol. Phys. , 547 (2012). T. A. H. Burd and D. C. Clary,J. Chem. Theory Comput. , 3486 (2020). N. Makri and W. H. Miller, J. Chem. Phys. , 4026 (1989). S. Coleman, in
Proc. Int. School of Subnuclear Physics (Erice,1977) also in S. Coleman,
Aspects of Symmetry , chapter 7, pp.265–350 (Cambridge University Press, 1985). A. I. Vainshtein, V. I. Zakharov, V. A. Novikov, and M. A.Shifman, Sov. Phys. Uspekhi , 195 (1982), also in Instantonsin Gauge Theories , edited by M. Shifman, pp. 468 (Singapore:World Scientific, 1994). W. H. Miller, J. Chem. Phys. , 1899 (1975). V. A. Benderskii, D. E. Makarov, and C. A. Wight,
Chemical Dynamics at Low Temperatures , Adv. Chem. Phys.,Vol. 88 (Wiley, New York, 1994). W. Siebrand, Z. Smedarchina, M. Z. Zgierski, and A. Fernández-Ramos, Int. Rev. Phys. Chem. , 5 (1999). Z. Smedarchina, W. Caminati, and F. Zerbetto, Chemicalphysics letters , 279 (1995). Z. Smedarchina, W. Siebrand, and A. Fernández-Ramos,J. Chem. Phys. , 224105 (2012). V. A. Benderskii, E. V. Vetoshkin, L. Von Laue, and H. P.Trommsdorff, Chem. Phys. , 143 (1997). G. V. Mil’nikov and H. Nakamura,J. Chem. Phys. , 6881 (2001). G. V. Mil’nikov and H. Nakamura, J. Chem. Phys. , 124311(2005). J. O. Richardson and S. C. Althorpe,J. Chem. Phys. , 054109 (2011). C. Vaillant and M. T. Cvitaš,Phys. Chem. Chem. Phys. , 26809 (2018). J. O. Richardson, D. J. Wales, S. C. Althorpe, R. P.McLaughlin, M. R. Viant, O. Shih, and R. J. Saykally,J. Phys. Chem. A , 6960 (2013). A. Garg, Am. J. Phys. , 430 (2000). C. Herring, Rev. Mod. Phys. , 631 (1962). L. D. Landau and E. M. Lifshitz,
Quantum Mechanics: Non-Relativistic Theory , 2nd ed. (Pergamon Press, Oxford, 1965). V. Benderskii, E. Vetoshkin, L. von Laue, and H. Trommsdorff,Chem. Phys. , 143 (1997). V. Benderskii, E. Vetoshkin, I. Irgibaeva, and H. Trommsdorff,Chem. Phys. , 393 (2000). Z. Smedarchina, W. Siebrand, and M. Z. Zgierski,J. Chem. Phys. , 1203 (1996). A. Fernández-Ramos, Z. Smedarchina, M. Z. Zgierski, andW. Siebrand, J. Chem. Phys. , 1004 (1998). G. Mil’nikov, O. Kühn, and H. Nakamura, J. Chem. Phys. ,074308 (2005). G. V. Mil’nikov, T. Ishida, and H. Nakamura,J. Phys. Chem. A , 5430 (2006). M. Eraković, C. L. Vaillant, and M. T. Cvitaš,J. Chem. Phys. , 084111 (2020). M. T. Cvitaš and S. C. Althorpe,J. Chem. Theory Comput. , 787 (2016). M. T. Cvitaš, J. Chem. Theory Comput. , 1487 (2018). H. Kleinert,
Path Integrals in Quantum Mechanics, Statistics,Polymer Physics and Financial Markets , 5th ed. (World Scien-tific, Singapore, 2009). V. A. Benderskii, S. Y. Grebenshchikov, and G. V. Mil’nikov,Chem. Phys. , 1 (1995). W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flan-nery,
Numerical Recipes: The Art of Scientific Computing , 3rded. (Cambridge University Press, Cambridge, 2007). T. Kawatsu and S. Miura, J. Chem. Phys. , 024101 (2014). J. C. Light, I. P. Hamilton, and J. V. Lill, J. Chem. Phys. ,1400 (1985). V. Babin, C. Leforestier, and F. Paesani,J. Chem. Theory Comput. , 5395 (2013). V. Babin, G. R. Medders, and F. Paesani,J. Chem. Theory Comput. , 1599 (2014). Y. Watanabe, T. Taketsugu, and D. J. Wales, J. Chem. Phys. , 5993 (2004). L. B. Braly, J. D. Cruzan, K. Liu, R. S. Fellers, and R. J.Saykally, J. Chem. Phys. , 10293 (2000). E. N. Karyakin, G. T. Fraser, and R. D. Suenram, Mol. Phys. , 1179 (1993). T. T. Nguyen, E. Székely, G. Imbalzano, J. Behler,G. Csányi, M. Ceriotti, A. W. Götz, and F. Paesani,J. Chem. Phys. , 241725 (2018). R. Courant and D. Hilbert,