Non-adiabatic transitions in multiple dimensionsSubmitted to the editors DATE.      Funding: T. Hurst was supported by The Maxwell Institute Graduate School in Analysis and its Applications, a Centre for Doctoral Training funded by the UK Engineering and Physical Sciences Research Council (grant EP/L016508/01), the Scottish Funding Council, Heriot-Watt University and the University of Edinburgh.

# Non-adiabatic transitions in multiple dimensions††thanks: Submitted to the editors DATE.      Funding: T. Hurst was supported by The Maxwell Institute Graduate School in Analysis and its Applications, a Centre for Doctoral Training funded by the UK Engineering and Physical Sciences Research Council (grant EP/L016508/01), the Scottish Funding Council, Heriot-Watt University and the University of Edinburgh.

V. Betz Fachbereich Mathematik, Technische Universität Darmstadt, 64289 Darmstadt, Germany (betz@mathematik.tu-darmstadt.de, http://www.mathematik.tu-darmstadt.de/b̃etz/).    B. D. Goddard School of Mathematics and the Maxwell Institute for Mathematical Sciences, University of Edinburgh, Edinburgh, UK, EH9 3FD (b.goddard@ed.ac.uk, http://www.maths.ed.ac.uk/b̃goddard/).    T. Hurst Maxwell Institute for Mathematical Sciences, School of Mathematics, University of Edinburgh, Edinburgh, UK, EH9 3FD (t.hurst@sms.ed.ac.uk, http://www.maxwell.ac.uk/migsaa/people/tim-hurst).
###### Abstract

We consider non-adiabatic transitions in multiple dimensions, which occur when the Born-Oppenheimer approximation breaks down. We present a general, multi-dimensional algorithm which can be used to accurately and efficiently compute the transmitted wavepacket at an avoided crossing. The algorithm requires only one-level Born-Oppenheimer dynamics and local knowledge of the potential surfaces. Crucially, in contrast to standard methods in the literature, we compute the whole wavepacket, including its phase, rather than simply the transition probability. We demonstrate the excellent agreement with full quantum dynamics for a range of examples in two dimensions. We also demonstrate surprisingly good agreement for a system with a full conical intersection.

AMS subject classifications. 35Q40, 81V55

## 1 Introduction

Many computations in quantum molecular dynamics rely on the Born-Oppenheimer Approximation (BOA) [13], which utilises the small ratio of electronic and reduced nuclear masses to replace the electronic degrees of freedom with Born-Oppenheimer potential surfaces. When these surfaces are well separated, the BOA further reduces computational complexity by decoupling the dynamics to individual surfaces.

However, there are many physical examples [15],[16],[27] and [31] where the Born-Oppenheimer surfaces are not well separated or even have a full intersection. In these regions the BOA breaks down, and the coupled dynamics must be considered; when a wavepacket travels over a region where the surfaces are separated by a small but none-vanishing amount, a chemically crucial portion of the wavepacket can move to a different energy level via a non-adiabatic transition. The existence of the small parameter introduced several challenges when attempting to numerically approximate the dynamics. First, and independent of the existence of an avoided or full crossing, the wavepacket oscillates with frequency and hence a very fine computational grid is required. Furthermore, in the region of an avoided crossing, the dynamics produce rapid oscillations; the transmitted wavepacket very close to the crossing is , but in the scattering regime the transmission is exponentially small. It is therefore necessary to travel far from the avoided crossing with a small time-step to accurately calculate the phase, size and shape of the transmitted wavepacket. In order to calculate the exponentially small wavepacket, one must ensure that the absolute errors in a given numerical scheme are also exponentially small, or they will swamp the true result. Finally, the number of gridpoints in the domain increases exponentially as the dimension of the system increases. Thus standard numerical algorithms quickly become computationally intractable.

Many efforts have been made to avoid computational expense by approximating the transmitted wavepacket while avoiding the coupled dynamics. Surface hopping algorithms discussed in [32, 25, 28, 22, 30, 20, 26, 17, 18, 24, 4, 3] approximate the transition using classical dynamics, where the Landau-Zener transition rate [33], [23] is used to determine the size of the transmitted wavepacket. This method has enjoyed some success, and has recently been applied to higher dimensional systems [24, 4]. However, the full transmitted quantum wavepacket is not calculated; phase information is lost. Such information is crucial when considering systems with interference effects, e.g. ones in which the initial wavepacket makes multiple transitions through an avoided crossing. Recently, there have been efforts to include phase information in surface hopping algorithms [14]. In contrast, in [10] and [7], a formula is derived to accurately approximate the full transmitted wavepacket, in one dimension, using only decoupled dynamics. The formula has been applied to a variety of examples with accurate results, including the transmitted wavepacket due to photo-dissociation of sodium iodide [9].

In this paper we construct a method to apply the formula derived in [10] and [7] to higher dimensional problems. We begin in LABEL:Section:1DFormula and LABEL:Section:Superadiabatic by outlining the derivation of the formula in one dimension [10], which involves deriving and approximating algebraic differential recursive equations for the quantum symbol of the coupling operator in superadiabatic representations. We extend these derivations to dimensions in LABEL:Section:NDCouplingOperators. In LABEL:Section:SlicingAlgorithm we create an -dimensional formula for systems which are slowly varying in all but one dimension, then extend this result via a simple algorithm to obtain a general -dimensional formula. We provide some examples and results in LABEL:Section:Numerics and note conclusions and future work in LABEL:Section:Conclusions.

Consider the evolution of a semiclassical wavepacket in dimensions with at time , , governed by the following equation:

 iε∂tψ(x,t)=Hψ(x,t), (2.1)

where is the ratio between an electron and the reduced nuclear mass of the molecule, i.e. . This system is derived after a rescaling of a two level Schrödinger equation [19]. Note that as we are considering semiclassical wavepackets, the derivatives of which are of order . The Hamiltonian of a two level system is given by [7]

 H =−ε22∇2xI+V(x)+d(x)I, (2.2)

where

 V(x)=(Z(x)X(x)X(x)−Z(x)) (2.3)

and is the part of the potential operator with non-zero trace. In general can be given by a Hermitian matrix, but as noted in [5], any Hermitian can be transformed into real symmetric form. It is useful to define , so that

 cos(θ(x))=Z(x)√X(x)2+Z(x)2,sin(θ(x))=X(x)√X(x)2+Z(x)2.

Then, defining , gives

 V(x) =ρ(x)(cos(θ(x))sin(θ(x))sin(θ(x))−cos(θ(x))). (2.4)

This is known as the diabatic representation of the system. We define and as the two diabatic potentials, with the diabatic coupling element as the off-diagonal element . Consider the unitary matrix which diagonalises the potential operator :

 U0(x)=⎛⎜⎝cos(θ(x)2)sin(θ(x)2)sin(θ(x)2)−cos(θ(x)2)⎞⎟⎠. (2.5)

If we define , then we arrive at the adiabatic Schrödinger equation

 iε∂tψ0(x,t)=H0ψ0(x,t). (2.6)

Here is given by

 H0 (2.7)

The adiabatic potential surfaces are given by the diagonal entries of the adiabatic potential matrix to leading order,

 VU=ρ(x)+d(x),VL=−ρ(x)+d(x), (2.8)

where is the upper adiabatic potential surface, and is the lower adiabatic potential surface. Assuming the initial wavepacket is purely on the upper level, the adiabatic representation approximates the transmitted wavepacket to leading order by the perturbative solution [29]

 ψ−0(t)=−iε∫t−∞e−iε(t−s)H−κ−1(x)⋅(ε∂x)e−iεsH+ϕds, (2.9)

where

 H±=−ε22∇2x±ρ(x)+d(x),κ±1(x)=±∂xθ(x)2. (2.10)

The perturbative solution in the adiabatic representation does not offer much explanation to the properties of the transmitted wavepacket. For instance, the constructed wavepacket at first looks to be . However due to the adiabatic coupling operator , fast oscillations and cancellations between upper and lower transmissions occur near the avoided crossing, so that far from the crossing the transmitted wavepacket is much smaller than the transition at the crossing point (LABEL:fig:MassDown).

For this reason, the transmitted wavepacket is better approximated using the perturbative solution from the superadiabatic representation [10], for some optimal choice of . The superadiabatic representation is produced by creating and applying unitary pseudodifferential operators , such that the off-diagonal elements of the potential operator have prefactor , and the diagonal elements are the same as in the adiabatic representation. Existence of such operators is discussed in [10]. In the superadiabatic representation the perturbative solution gives

 ψ−n(t)=−iεn∫t−∞e−iε(t−s)H−nK−n+1(x)e−iεsH+nϕds, (2.11)

where is the Hamiltonian in the superadiabatic representation, given by

 (2.12)

for some pseudodifferential coupling operators , and

 H±n =−ε22∇2x±ρ(x)+d(x). (2.13)

Unfortunately, the need to compute to compute the pseudodifferential operators and prevent this from directly producing a practical numerical scheme. However, as we now demonstrate, we may make use of the superadiabatic representations to obtain a simple and accurate algorithm.

## 3 Approximating the transition in one dimension

The derivation in [7] requires and to be analytic in a strip containing the real axis. The formula is derived in one dimension using the superadiabatic perturbative solution by

1. Finding algebraic recursive differential equations to calculate the quantum symbol , where is the Weyl quantisation of ,

 (W(κ±n+1)ψ)(x) =12π∫R2ndξdyK±n+1(ξ,12(x+y))eiξ⋅(x−y)ψ(y). (3.1)
2. Rescaling by

 τ(q)=2∫q0ρ(r)dr, (3.2)

then approximating in an analogous way to the time-adiabatic case in [11].

3. Assuming the potential surfaces are approximately linear near the avoided crossing, i.e. , to utilise the Avron-Herbst formula [1].

4. Applying a stationary phase argument to evaluate the remaining integral.

The result is presented in scaled momentum space: in dimensions the wavepacket in scaled momentum space is given using the -scaled Fourier transform

 ˆfε(p)=1(2πε)d/2∫Rdf(x)exp(−iεp⋅x)dx. (3.3)

Following this derivation leads to an approximation of the transmitted wavepacket, far from the avoided crossing:

 ˆψ−ε(k,t) =e−iεtˆH−(k)ν(k)+k2|ν(k)|e−τc2δε|k−ν(k)|e−iτr2δε(k−ν(k))χk2>4δˆϕ+0ε(ν(k)), (3.4)

where

• There is no dependence on the superadiabatic representation used in the formula derivation.

• is the wavepacket on the upper level evolved to the avoided crossing using uncoupled dynamics.

• , half the distance between the two adiabatic potential surfaces at the avoided crossing.

• , the initial momentum a classical particle would need to have momentum after falling down a potential energy difference of , i.e. the distance between the potential surfaces at the avoided crossing, which shifts the wavepacket in momentum space.

• , where is the smallest complex zero of , when extended to the complex plane. The prefactor determines the size of the transmitted wavepacket. In [19], we show that under appropriate approximations of the momentum and potential surfaces, this prefactor is comparable to the Landau-Zener transition prefactor used in surface hopping algorithms such as in [4]. An additional change in phase occurs due to , which is present when the potential is not symmetric about the avoided crossing.

The constructed formula LABEL:eq:1DFormula allows us to approximate the size and shape of the transmitted wave packet due to an avoided crossing, and avoid computing expensive coupled dynamics. The method for applying the algorithm is as follows:

1. Begin with an initial wave packet on the upper adiabatic energy surface, far from the crossing, with momentum such that the wave packet will cross the minimum of (LABEL:fig:1DToyExampleFormulaA).

2. Evolve according to the BOA on the upper adiabatic level until the centre of mass is at the avoided crossing, at time (LABEL:fig:1DToyExampleFormulaB), say ,

3. Apply the one dimensional formula to the -Fourier transform of the wave packet at the crossing (LABEL:fig:1DToyExampleFormulaC):

 ˆψ−ε(x,tcz)=ν(k)+k|ν(k)|e−τc2δε|k−ν(k)|e−iτr2δε(k−ν(k))χk2>4δˆϕ+0ε(ν(k),tcz), (3.5)
4. Evolve the transmitted wave packet far away enough from the crossing, say to time , using the BOA (LABEL:fig:1DToyExampleFormulaD): . At this time the transmitted wave packet in momentum space should be accurately approximated by the formula result, as long as the wave packet is evolved away from the area in which oscillations occur.

Note that we have assumed that the avoided crossing is centred at . When the avoided crossing is not centred at 0, an additional shift term , where is the position of the avoided crossing, is obtained when changing variables in the Fourier transform of the coupling symbols [19].

Applications of the one dimensional formula have been widely successful on a variety of examples. In [9], the formula is used to accurately approximate the transmitted wavepacket for sodium iodide. Tilted avoided crossings have also been examined, and a formula developed which is dependent on the value , so the optimal superadiabatic representation must be calculated. The formula has also been sucessfully applied to model interference effects in multiple transitions [8].

Finally, following the above derivation for reverse transitions (from lower to upper surface), the formula

 ˆψ+ε=−~ν+k2|~ν|e−τc|k−~ν|/(2δε)e−iτr(k−~ν)/(2δε)ˆψ+0ε(~ν(k)), (3.6)

where , can be used to approximate the wavepacket transmitted to the upper surface, far from the avoided crossing.

## 4 Coupling operators in higher dimensions

The first step in deriving LABEL:eq:1DFormula in [10] was to approximate the superadiabatic coupling operators . We now consider these operators in higher dimensions. We restrict the calculations here to two dimensions for clarity, but they can easily be adapted to dimensions.

###### Lemma 4.1.

In two dimensions, is given by

 κ±n+1(p,q)=−2ρ(q)(xn+1(p,q)±yn+1(p,q)). (4.1)

where are given by the following algebraic recursive differential equations:

 x1=z1=w1=0,y1=−i4ρ(p1∂q1θ+p2∂q2θ). (4.2)

and

 yn=0, n even,xn=zn=wn=0, n odd. (4.3)

For odd, we have

 xn+1=−12ρ[1i(p⋅∇qyn)−2n∑j=11(2i)jj!∑|α|=j∂αp(bαzn+1−j−aαxn+1−j)], (4.4)

and for even

 yn+1=−12ρ[1i((p⋅∇qxn)−zn(p⋅∇qθ))−2n∑j=11(2i)jj!∑|α|=j∂αp(−aαyn+1−j+bαwn+1−j)], (4.5)
 1i((p⋅∇qzn)−xn(p⋅∇qθ))=n∑j=11(2i)jj!∑|α|=j∂αp(bαyn+1−j+aαwn+1−j), (4.6) 1i(p⋅∇qwn)=2n∑j=11(2i)jj!∑|α|=j∂αp(aαzn+1−j+bαxn+1−j), (4.7)

where and , and and are given by the recursions

 a0=ρ(q1,q2),b0=0, a(α1+1,α2)=∂q1a(α1,α2)+(∂q1θ)b(α1,α2),b(α1+1,α2)=∂q1b(α1,α2)−(∂q1θ)a(α1,α2), a(α1,α2+1)=∂q2a(α1,α2)+(∂q2θ)b(α1,α2),b(α1,α2+1)=∂q2b(α1,α2)−(∂q2θ)a(α1,α2).

###### Proof.

The method is a straightforward extension of the theory in [10].

As in [10], by general theory, the coefficients are polynomials in of order . We therefore write

 xn(p,q)=n∑m=0m∑k=0pk1pm−k2xk,m−kn(q1,q2), (4.8)

for some , and similarly for . For a given , we write for each . Note that:

 ∂αjppk1pm−k2 ={k!(m−k)!(m−k)!(m−k−j+α)!pk−α1pm−k−j+α2,k≥α and m≥j,0,otherwise (4.9)

Then

 ∂αjpxn+1−j =n+1−j∑m=0m∑k=0(∂αp1pk1)(∂j−αp2pm−k2)xk,m−kn+1−j(q1,q2), =n+1−j∑m=jm−α+j∑k=αk!(k−α)!(m−k)!(m−k−j+α)!pk−α1pm−k−j+α2xk,m−kn+1−j(q1,q2),

so that

 A\coloneqqn∑j=11(2i)jj!j∑α=0∂αp1∂j−αp2aαjxn+1−j, (4.10)

can be rewritten as

 n∑j=11(2i)jj!j∑α=0aαjn+1−j∑m=jm+α−j∑k=αk!(k−α)!(m−k)!(m−k−j+α)!pk−α1pm−k−j+α2xk,m−kn+1−j.

We now want to extract and from the final two summations, so that we can compare coefficients on either side of the results of LABEL:Lemma:2DRecursions to construct recursive equations for for . Consider terms where . By the limits of the third summand, we find that , and that , a contradiction. Therefore we restrict the limits of first summand. Defining , and , we find

 A =⌊n+12⌋∑j=1j∑α=0n+1−2j∑c=0c∑b=0aαj(2i)jj!(b+α)!b!((c+j)−(b+α))!(c−b)!pb1pc−b2xb+α,(c+j)−(b+α)n+1−j.

We now want to switch the order of summation. We note that, for an arbitrary ,

 ⌊n+12⌋∑j=1n+1−2j∑c=0Bc,j=n+1∑c=0⌊c2⌋∑j=1Bn+1−c,j,

which can be shown directly (note that the terms where , are zero). Using this, we finally have that

 A=n+1∑c=0n+1−c∑b=0pb1pn+1−c−b2×⌊c2⌋∑j=1j∑α=0aαj(2i)jj!(b+α)!b!(n+1−c+j−b−α)!(n+1−c−b)!xb+α,(n+1−c+j)−(b+α)n+1−j. (4.11)

Importantly, and have been extracted from two of the summations. Note that taking and , or and in LABEL:eq:ASum, we return the 1D result in [10] for and respectively. We obtain the following result.

###### Proposition 4.2.

The coefficients to are determined by the following algebraic-differential recursive equations. We have

 xA,B1=zA,B1=wA,B1=0,A+B∈{0,1}, (4.12) y0,01=y1,11=0,y1,01=−i4ρ∂q1θ,y0,11=−i4ρ∂q2θ. (4.13)

Further,

 xA,Bn+1=−12ρ[1i(∂q1yA−1,Bn+∂q2yA,B−1n)−2⌊n+1−(A+B)2⌋∑j=1j∑α=01(2i)jj!×(A+α)!A!(B+j−α)!B!(bαjzA+α,B+j−αn+1−j−aαjxA+α,B+j−αn+1−j)], (4.14)

When is odd. When is even, we have

 yA,Bn+1=−12ρ[1i((∂q1xA−1,Bn+∂q2xA,B−1n)−(zA−1,Bn∂q1θ+zA,B−1n∂q2θ))−2⌊n+1−(A+B)2⌋∑j=1j∑α=01(2i)jj!(A+α)!A!(B+j−α)!B!×(−aαjyA+α,B+j−αn+1−j+bαjwA+α,B+j−αn+1−j)], (4.15)
 0=1i((∂q1zA−1,Bn+∂q2zA,B−1n)+(xA−1,Bn∂q1θ+xA,B−1n∂q2θ))−2⌊n+1−(A+B)2⌋∑j=1j∑α=01(2i)jj!(A+α)!A!(B+j−α)!B!×(bαjyA+α,B+j−αn+1−j+aαjwA+α,B+j−αn+1−j), (4.16)
 0=1i((∂q1wA−1,Bn+∂q2wA,B−1n)−2⌊n+1−(A+B)2⌋∑j=1j∑α=01(2i)jj!(A+α)!A!(B+j−α)!B!×(aαjzA+α,B+j−αn+1−j+bαjxA+α,B+j−αn+1−j). (4.17)

###### Proof.

We substitute LABEL:eq:xyzwExpansion into the results of LABEL:Lemma:2DRecursions and compare coefficients in powers of on either side, using LABEL:eq:ASum.

As with the coefficients and in LABEL:eq:KappaWithxy, has polynomial form:

 κ±n+1(p,q)=n∑m=0m∑j=0pj1pm−j2κ(j,m−j)±n+1(q1,q2). (4.18)

The coupling operators only act on wavepackets near the avoided crossing, so when considering the effect of the coupling operator on a wavepacket we assume that the path of the wavepacket is approximately linear. Furthermore, by a change of coordinate system, we may assume that the direction of travel of the wavepacket is independent of (i.e. we rotate the frame of reference so that the wavepacket is moving in the direction). We also assume that the leading order term in dominates : . This is similar to the assumption made in the one dimensional case, where it can be shown to be accurate for sufficiently large , but in practice holds for much smaller values. Then the 2D algebraic differential recursive equations then reduce to the one dimensional case in [10]:

 xn+1,0n+1≈i2ρ(∂q1yn,0n), yn+1,0n+1≈i2ρ((∂q1xn,0n)′−(∂q1θ)zn,0n),0≈∂q1zn,0n+(∂q1θ)xn,0n. (4.19)

To ease notation, redefine , and similar for . It is unclear what the analogue of LABEL:eq:1DTau, introduced initially in [6] for the time-adiabatic case, would be for multidimensional systems. We introduce the natural scaling in the first dimension

 τ(q1,q2)=2∫q10ρ(r,q2)dr. (4.20)

Defining the recursive relations LABEL:eq:xyzwLeading then become

 ~x0n+1=i~y0n+1,~y0n+1=i((~x0n)′+~θ′~z0n),0=(~z0n)′+~θ′~x0n, (4.21)

where . These recursive equations also occur in [11], where they are solved in one dimension, under the assumption that

 ddτ~θ(τ) =iγτ−¯τcz−iγτ−τcz+~θ′r(τ), (4.22)

where is a first order complex singularity of , and has no singularities closer to the real axis than . Assuming that the avoided crossing occurs at , we can write , for some analytic function such that , and is quadratic in the neighbourhood of . Therefore a Stokes line (i.e. a curve with Im) crosses the real axis perpendicularly [21], and following this line leads to a pair of complex conjugate points which are complex zeros of . Defining , it is shown in [6] that first order complex singularities of the adiabatic coupling function arise at these complex zeros. This derivation is still valid in our case, for each . The recursive algebraic differential equations solved in [11] then give us to leading order:

 κ−n(q)≈inπρ(q)(n−1)!(i(τ(q)−¯τcz(q2))n−i(τ(q)−τcz(q2))n). (4.23)

It is clear that the results of this section can be extended to higher dimensions, by assuming the direction of travel of the wavepacket is in the first dimension. We will now use this observation to design an algorithm for multi-dimensional transitions using only the 1D transition formula.

## 5 The multi-dimensional algorithm

The derivation of a multidimensional formula, under the assumptions above, follow similarly to the one dimensional case. First we define the multidimensional Weyl quantization [2].

###### Definition 5.1.

For a symbol , given a test function , we define the Weyl quantization of by

 (WεHψ)(x)=1(2πε)d∫R2ddξdyH(ε,ξ,12(x+y))eiε(ξ⋅(x−y))ψ(y). (5.1)

We want to approximate the pseudodifferential operator , which is given by the Weyl quantisation of . The particular form of allows us to simplify the Weyl quantisation as follows.

###### Proposition 5.2.

Let , for . Then

 ˆ(Wεκ)(ψ)ε(k)=1(2πε)d/2∫Rdˆgε(k−η)d∏i=1(ki+ηi2)Aiˆψε(η)dη. (5.2)

###### Proof.

Firstly, using that ,

 Wεκψ(x)= 1(2πε)d∫R2ddξdy(d∏i=1ξAii)g(x+y2)eiε(ξ⋅(x−y))ψ(y), = 1(2πε)3d/2∫R3ddξdydη(d∏i=1ξAii)g(x+y2)eiε(ξ⋅(x−y)+η⋅y)ˆψε(η).

Now define , . Then

 Wεκψ(x) =2d(2πε)3d/2∫R3ddξd~ydη(d∏i=1ξAii)g(~y)eiε(ξ⋅x+(2~y−x)⋅(η−ξ))ψ(η), =2d(2πε)d∫R2ddξdη(d∏i=1ξAii)eiε(x⋅(2ξ−η))ψ(η)^gε(2(ξ−η)).

We perform a second change of variables and find

 Wεκψ(x) =1(2πε)d