## Abstract

We study the dynamics of a molecule’s nuclear wave function near an avoided crossing of two electronic energy levels for one nuclear degree of freedom. We derive the general form of the Schrödinger equation in the *n*th superadiabatic representation for all . Using these results, we obtain closed formulas for the time development of the component of the wave function in an initially unoccupied energy subspace when a wave packet travels through the transition region. In the optimal superadiabatic representation, which we define, this component builds up monotonically. Finally, we give an explicit formula for the transition wave function away from the avoided crossing, which is in excellent agreement with high-precision numerical calculations.

## 1. Introduction

The time-dependent Born–Oppenheimer (BO) theory is the single most important tool for studying the quantum dynamics of molecules, just as the time-independent BO approximation (Born & Oppenheimer 1927) is for the theory of molecular bound states. The basic physical idea is that the electrons, being much lighter and thus faster than the nuclei, quickly adjust their position with respect to the latter. When starting in the *n*th bound state (for some fixed positions of nuclei), they should remain in the *n*th bound state even though the nuclei are slowly moving; this will then be the *n*th bound state with respect to the updated nuclear position. In turn, the nuclear quantum dynamical motion is determined by an effective potential given by the energy level of the *n*th bound state of the electrons, as a function of nuclear position. The theory is expected to be valid up to errors of order *ε*, the square root of the mass ratio between electrons and nuclei. It was first proposed by London (1928), and a version of it was first proved by Hagedorn (1980). Since then, there have been many extensions and improvements. Instead of trying to give a full account of all of them, we refer to the excellent review article (Hagedorn & Joye 2007).

Despite the success of the time-dependent BO theory, there are important situations where it fails. This happens when electronic energy levels for different bound states are not well separated. Then, either the energy levels actually cross (Hagedorn 1994; Cederbaum *et al.* 2005; Lasser & Teufel 2005), in which case the BO approximation breaks down completely, or they approach each other very closely without actually crossing, in which case the BO approximation is still valid to leading order (and, as it will turn out, even beyond all orders in the scattering regime); but the remaining deviations are still of great interest. In fact, it is an avoided crossing situation that leads to the photo-dissociation of NaI, one of the paradigmatic chemical reactions in photochemistry (Zewail 1994). The present article studies the part of the wave function that performs a transition between electronic energy levels at an avoided crossing.

We treat a two-level system with one nuclear degree of freedom. Its Hamiltonian is
1.1** I** is the 2×2 unit matrix, and

*H*acts on wave functions , i.e. square integrable functions with values in . is the nuclear coordinate,

*ε*>0 is the square root of the mass ratio and

*V*is the (adiabatic) potential energy matrix, written in polar form. For the moment, we assume that

*ρ*(

*x*) and

*θ*(

*x*) are analytic in a strip containing the real axis. More detailed conditions will follow later. We also assume that

*(*

*ρ**x*)≥

*δ*>0 for all . This corresponds to an avoided crossing situation with gap at least 2

*δ*. The time evolution is given by the time-dependent Schrödinger equation 1.2 in the time scale where nuclei move a distance of order 1 in a time of order 1.

It is natural to study equation (1.2) in a representation where the matrix *V* becomes approximately diagonal. The most fundamental of these is the adiabatic representation. It is implemented by the unitary transform *U* acting on , where *U**f*(*x*)=*U*_{0}(*x*)*f*(*x*), with
1.3
Putting *ψ*_{a}(*x*,*t*)=*U*_{0}(*x*)*ψ*(*x*,*t*), we obtain the Schrödinger equation in the adiabatic representation,
1.4
with
1.5
To leading order in *ε*, the dynamics in the two components of *ψ*_{a} decouple and are Schrödinger evolutions with potentials given by the energy bands ±* ρ*(

*x*): this is the BO approximation. Its proof is not trivial as there is also an

*ε*on the left-hand side of equation (1.4), see Teufel (2003) for details. Couplings between the components are given, to leading order, by the first terms in the off-diagonal of equation (1.5). Note that because acts on semiclassical wave functions that oscillate with frequency 1/

*ε*, the operator

*ε*∂

_{x}is actually of order 1. These couplings are called derivative couplings in quantum chemistry, for obvious reasons.

While transitions between the components of *ψ*_{a} are of order *ε* globally, they are usually exponentially small in the scattering regime. More precisely, it is known that, under suitable assumptions, for every there exist unitaries *U*_{n} on such that the components of solutions *ψ*_{n} of the corresponding transformed equation decouple up to errors of order *ε*^{n}. By optimizing over *n*, decoupling up to exponentially small (in *ε*) errors can be shown. The representation induced by *U*_{n} is sometimes called the *n*th superadiabatic representation. If *V* approaches a constant near sufficiently quickly, *U*_{n} agrees with *U*_{0} up to errors involving the derivative of *V* . We refer to Hagedorn & Joye (2001) and Martinez & Sordoni (2002, 2009) for details.

The weakness of the above results is that they only give upper bounds. Thus, they cannot answer the following natural and practically important questions: assume that the initial condition of equation (1.4) lies in the first adiabatic subspace, i.e. *ψ*_{a,0}(*x*)=(*ψ*_{+,0}(*x*),0)^{t}, and that *ψ*_{+,0}(*x*) is concentrated away from regions where *V* varies significantly. After one transition through a region where *V* varies, how large is the transition probability into the second adiabatic subspace? i.e. how large is the *L*^{2} norm of the second component *ψ*_{−}(*x*,*t*) of *ψ*_{a}(*x*,*t*)? What is the precise form of *ψ*_{−}(*x*,*t*)? How does the second component of the wave function in the optimal superadiabatic representation evolve in time? In this article, we give answers to all of these questions.

Our main rigorous result is theorem 3.4, where we derive the leading-order symbol *H*_{n}(*ε*,*p*,*q*) of the Hamiltonian in the *n*th superadiabatic representation. We do not currently control the full asymptotic behaviour of *H*_{n}(*ε*,*p*,*q*) for . However, the asymptotics of a non-rigorous approximation of that symbol are known and lead to the result of ours that we consider the most relevant in practice: an explicit formula for the exponentially small wave packet that makes the non-adiabatic transition in terms of data that are local in space and time. We specialize to the case * ρ*(

*x*)=

*δ*in equation (1.1), i.e. constant eigenvalues, and assume that the adiabatic coupling function

*θ*′(

*x*) appearing in equation (1.5) has the generic form (4.9). Let

*x*

_{t}be where

*θ*′(

*x*) has a global maximum in the absolute value. We will refer to

*x*

_{t}as the transition point and change variables to achieve

*x*

_{t}=0. Let be the Fourier transform of the wave packet moving according to the BO approximation in the upper electronic surface, at the time when its maximum reaches the transition point. We will call this time the transition time. The evolution in the lower electronic surface can then be approximated as follows: at the transition time, we start a wave packet with Fourier transform given by 1.6 in the lower electronic surface and evolve it according to the (decoupled) BO approximation for that surface. Above, is the momentum before the transition and 2

*δ*the energy gap, and both γ and

*q*

_{c}can be calculated from

*θ*′. We refer to §4 and especially equation (4.12) for more details.

In the *n*th superadiabatic representation, the approximation (1.6) can be used shortly after the transition time. In the adiabatic region, it is valid to the extent that the adiabatic and superadiabatic representations agree, which will usually be away from the avoided crossing region. Thus, we can extract from it exact asymptotic adiabatic transition probabilities, as well as the asymptotic form of the transition wave packet. Consequently, equation (1.6) provides a practical way to correctly include non-adiabatic transitions at avoided crossings into the BO approximation, for the case of constant eigenvalues. An analogous formula for Landau–Zener-like avoided crossings will be the content of a future publication.

As indicated above, the derivation of equation (1.6) is not rigorous. In fact, equation (1.6) is most likely not even asymptotically correct as *ε*→0, for reasons that we will discuss. For moderate values of *ε* and the incoming momentum, it is, on the other hand, an excellent approximation to the true dynamics. We show this by comparing the transmitted wave function as per equation (1.6) with high-precision *ab initio* numerical calculations. Not only does equation (1.6) correctly describe the transition probability, but it yields the shape of the wave function after the transition to such a high accuracy that, e.g. for the situation plotted in figure 2, the relative error for from formula (1.6) and from a highly accurate numerical solution of the Schrödinger equation is of the order 10^{−5} uniformly in those *k* on which is essentially supported. Put differently, figure 2 shows both the true solution and the approximation by our formula (1.6).

Equation (1.6) addresses the questions of transition probability and the shape of the transmitted wave function, but it does not give a hint about the time development of the wave packet during the transition. This question only really makes sense in superadiabatic representations, as, for example, in the adiabatic basis, transitions are of order *ε* and highly oscillatory near the transition point, and only decay later. In the time-adiabatic simplification of the problem (see below), it is known (Berry 1990; Betz & Teufel 2005*b*) that there exists an optimal superadiabatic representation in which transitions are (uniformly in time) exponentially small, and have the universal shape of an error function. For the present problem, it is not obvious how to identify optimality of superadiabatic representations. We will propose a novel, natural criterion for optimality, cf. definition 4.2; when our criterion is met, transition wave packets in the Fourier representation have the universal error-function shape as functions in time near their maximum. We also show that in a simple but important case our criterion can be met.

Owing to the intricacy of the problem of transitions at avoided crossings, not many mathematical approaches exist to date. In Hagedorn & Joye (1998) and Rousse (2004), a situation is considered where the gap between the energy levels shrinks with the parameter *ε*; this can be seen as an intermediate case between a crossing and a true avoided crossing. The most relevant work with respect to ours is by Hagedorn & Joye (2005), where another formula is given (and proved) for the asymptotic shape of a non-adiabatic scattering wave function at a true avoided crossing. Their formula looks quite different from ours. We do not state it here, as this would introduce too much additional notation, and instead refer to theorem 5.1 of Hagedorn & Joye (2005). But we want to comment on the connections between their result and ours. The first question should certainly be whether the two formulae agree. However, equation (1.6) is for constant * ρ*, while Hagedorn & Joye (2005) cover the Landau–Zener situation, and so a direct comparison is not possible at present. But even for a Landau–Zener version of equation (1.6) the answer would likely be no because Hagedorn and Joye’s formula is asymptotically correct while ours is not. But we know that equation (1.6) works extremely well for physically relevant parameters, and the great advantage over the results of Hagedorn & Joye (2005) is that equation (1.6) is much easier to implement in practice. The formula in Hagedorn & Joye (2005) contains complex contour integrals of the analytic continuation of some function of

*V*, and assumes a certain form (semiclassical wave packet) for the incoming wave function. When their assumptions are satisfied, the transmitted wave packet is always approximately Gaussian. By contrast, equation (1.6) is a simple multiplication in Fourier space, and works well for general incoming wave functions, including cases where the transmitted wave function is clearly non-Gaussian. Thus, while the results of Hagedorn & Joye (2005) are better than ours from the point of view of asymptotic analysis, ours seem better suited for practical exploitation.

Let us finally put our results into perspective with respect to a theory that has been extensively studied over many decades. It is the time-adiabatic version of equation (1.2), which is usually just referred to as the model of adiabatic quantum transitions. The simplification consists in prescribing a (classical) path *x*(*t*) for the nuclear motion of the two-level system, instead of considering its quantum evolution. The equation then reduces to
1.7
which is to be read as an equation for the occupation probabilities of the two relevant electronic energy bands given a certain motion of the nuclei. Zener (1932) investigated an explicitly solvable instance of equation (1.7), namely *X*(*t*)=*t*/2 and *Z*(*t*)=*δ*/2. He observed that by solving equation (1.7) with an initial condition parallel to the eigenvector corresponding to +* ρ*(

*t*) at , the solution at would have a component of magnitude e

^{−πδ2/(4ε)}in the eigenspace corresponding to −

*(*

*ρ**t*). In other words, in the scattering regime, the transition amplitude is exponentially small. Shortly after, Landau (1965) argued that the same exponentially small expression should describe the scattering regime also for general analytic

*X*(

*t*) and

*Z*(

*t*) such that, at the minimum

*t*

_{0}of

*ρ*^{2}(

*t*)=

*X*

^{2}(

*t*)+

*Z*

^{2}(

*t*),

*V*is approximated to second order by the one considered by Zener. Models of this type are now called Landau–Zener models, and the formula that predicts the magnitude of the exponentially small transition is known as the Landau–Zener formula. It has attracted much research, including Demkov

*et al.*(1978) where more general situations leading to different pre-factors were considered in generality, and Joye (1994), where a rigorous proof of generalized Landau–Zener transitions was given. We refer to Hagedorn & Joye (2009) for further information on the subject.

The exponentially small scattering amplitude cannot be explained easily by any method involving just the adiabatic subspaces, i.e. the instantaneous eigenspaces of *V* (*t*). Indeed, when starting the evolution in one of the eigenspaces and monitoring the component ϕ_{2}(*t*) parallel to the other one, one will typically observe a build-up of |ϕ_{2}(*t*)| up to order *ε*, and later an eventual decay to the exponentially small final value; see Lim & Berry (1991) and Betz & Teufel (2006) for a numerical illustration of the phenomenon. One way to understand better what is actually going on is the complex WKB method: one solves equation (1.7) not on the real line, but on a curve in the complex plane that approaches the real line at and passes through the zeroes of the complex continuation of * ρ*(

*t*). Then, the solution is then exponentially small all the way, and by this method Joye (1994) proves the Landau–Zener formula. For the true BO time evolution (1.2), Hagedorn and Joye extend these methods and obtain the transition formula given in Hagedorn & Joye (2005).

A second approach to the understanding of the exponentially small transition amplitudes was found by Berry (1990). In an influential paper (Berry 1990), he expands the solution of equation (1.7) into a formal power series, and by truncating the resulting asymptotic series after *n* terms, he obtains a time-dependent basis of , called the *n*th superadiabatic basis. For all *n*, these bases agree with the adiabatic basis in the scattering regime, i.e. when the eigenspaces are approximately constant. Choosing *n* so that the remainder term in the asymptotic expansion is minimal, Berry shows that not only are transitions between the corresponding subspaces exponentially small, but they are also universal in the sense that to leading order they are described by an error function for a wide class of matrices *V* (*t*). These results have later been made rigorous in Hagedorn & Joye (2004) and Betz & Teufel (2005*a*) for special, non-generic cases and in Betz & Teufel (2005*b*) for a large class of Hamiltonians including the generic case. In the latter two papers, a slight variation of Berry’s ideas is used: instead of expanding the solution of equation (1.7), the equation itself is transformed using the adiabatic perturbation theory as described in Teufel (2003), leading to superadiabatic representations. The present paper follows that same methodology, but the complications arising from trading a singularly perturbed ODE for a singularly perturbed PDE are considerable.

One might therefore ask why to bother and treat equation (1.2) when a nice effective system like equation (1.7) is available and, by now, pretty much well understood. One answer is that as equation (1.2) is the fundamental evolution equation of the quantum system, we had better understand what it does. More practically important is that the predictions made by the Landau–Zener formula are simply wrong even when applied to Gaussian wave packets. As is discussed after equation (4.12), and as was also observed by Hagedorn & Joye (2005), in true BO transitions, the outgoing momentum is larger than can be expected by energy conservation, and the transition amplitude is higher than predicted by the Landau–Zener formula. Therefore, a treatment that respects the quantum nature of nuclear movement, as carried out in the present work, is unavoidable in order to get quantitatively correct results.

## 2. Superadiabatic representations

We seek superadiabatic representations for equation (1.2). Intuitively, a superadiabatic representation of order *n* is implemented by an operator *U*_{n} in such that *U*_{n} is unitary and is an operator-valued diagonal matrix, each up to errors of order *ε*^{n+1}. Then, with **ϕ**_{n}=*U*_{n}**ϕ**, the equation decouples into two scalar equations, up to errors of order *ε*^{n+1}. The best known instance of this procedure is the adiabatic transformation (1.4). In that case, the *U*_{n}=*U*_{0} is just a matrix multiplication. In contrast, for *n*>0, *U*_{n} will be a pseudo-differential operator. Moreover, some care has to be exercised when speaking of errors of order *ε*^{n+1} above; as mentioned in the introduction, e.g. *ε*∂_{x} is not a small quantity in our context. For these reasons, it is necessary to work with semiclassical symbols. A full account of symbolic calculus can be found in Dimassi & Sjöstrand (1999) and Martinez (2002). Here, we only outline the main formulae that we are going to use, and do not touch the topic of symbol classes. Working in the symbolic representation means that we first replace *x* by and i*ε*∂_{x} by an independent variable in the definition (1.2) of *H*. The factor *ε* takes into account the semiclassical scaling. We then obtain
2.1
where *V* is as in equation (1.1). The requirement that translates into
2.2
for the symbol . Here # denotes the *Moyal product*: for two symbols *A*(*p*,*q*) and *B*(*p*,*q*), their Moyal product is defined through with
2.3
The Moyal product is the natural product for semiclassical symbols and accounts for the non-commutativity of the operators they represent. It can be extended to symbols depending on *ε*, like *U*_{n}, by representing *U*_{n} as a formal power series and collecting powers of *ε*, but we will only ever need formula (2.3) in explicit calculations. Once *H*_{n}(*ε*,*p*,*q*) is constructed, it can be recast as a pseudo-differential operator using the *Weyl quantization*: the operator corresponding to *H*_{n}(*ε*,*p*,*q*) acts on a test function ϕ by
2.4

For actually constructing the superadiabatic unitaries *U*_{n}, we use the method of *Superadiabatic projections*, that is, we seek symbols
2.5
such that
2.6
and
2.7
Above, [*A*,*B*]_{#}=*A*#*B*−*B*#*A* denotes the Moyal commutator.

The general theory (Teufel 2003) guarantees the existence of π^{(n)} for each *n*, and the same general theory states the following lemma.

## Lemma 2.1

There exists a semiclassical symbol *U*_{n}(*ε*,*p*,*q*) such that
2.82.9*and*2.10*where π*_{r}*is the projection onto the first component of*.

It is easy to see that equation (2.2) follows from equations (2.8)–(2.10) above, and thus *U*_{n} implements a superadiabatic representation. But we need to know more than that for understanding transitions between the electronic levels: we need to actually determine the leading order of the off-diagonal elements. This looks hopeless at first because the *U*_{n}, although given via a recursive scheme (Teufel 2003), are very tricky to calculate for large *n*. Fortunately, there is no need to do so owing to the following result.

## Proposition 2.2

*Define H*_{n}(*ε*,*p*,*q*) *as in equation* (2.2), *and assume the U*_{n}*given in that formula fulfils equations* (2.8)–(2.10). *Let F*_{n}*be given through equation* (2.7). *Denote by**the upper off-diagonal element of U*_{0}*F*_{n}*U*_{0}, *and let**be the lower off-diagonal element. Then*,

## Proof

The diagonal terms are immediate from first-order adiabatic perturbation theory. For the off-diagonal terms, we multiply equation (2.10) with *U*_{n} from the right and use equation (2.9) in order to find
2.11
for some symbol *R*, and similarly for *U*_{n}π_{r}. Iterating the above reasoning, we find
2.12
We now use equation (2.11) in order to get
The next to the last line above is *O*(*ε*^{2n+2}) by equation (2.12) and the fact (1−π_{r})π_{r}=0, and multiplying with π_{r} from the right, we find
This is the result for the upper-right off-diagonal element. For the lower left one, we interchange the roles of π_{r} and 1−π_{r} and obtain an additional minus sign from equation (2.7). ▪

## 3. Superadiabatic projections

Proposition 2.2 means that we can focus our attention entirely on superadiabatic projections. Here, we present increasingly explicit recursive schemes for calculating them. We start by deriving a matrix recursion for superadiabatic projections. It is easy to check that for π_{0}(*p*,*q*)=π_{0}(*q*)=(1/2)(1+*V* (*q*)/* ρ*(

*q*)), we have π

_{0}

*V*=

*π*

*ρ*_{0}. Hence, π

_{0}is the adiabatic (zeroth superadiabatic) projection, i.e. the projection on the ‘upper adiabatic’ subspace of

*V*corresponding to +

*. The same construction works when starting with the projection onto the lower adiabatic subspace.*

*ρ*In order to obtain the higher superadiabatic projections, we note that
Here, is the coeffcient of corresponding to *ε*^{k}. Using equation (2.3), we find
3.1
Similarly, equation (2.10) gives
3.2
Finally, we can calculate π_{n+1} through
3.3
The proof of equation (3.3) is given in Betz & Teufel (2005*a*, proposition 1), for the special case * ρ*=1/2 and

*F*

_{n+1}=−iπ

_{n}′(

*q*). As that proof applies verbatim to the general case, we do not repeat it here. Similarly, the important relations 3.4 and 3.5 follow as in Betz & Teufel (2005

*a*).

We need a more explicit recursive scheme in order to feasibly calculate the quantities *F*_{n}, *G*_{n} and π_{n} given in equations (2.5)–(2.10). We now introduce such a scheme, following Betz & Teufel (2005*a*), but using a different notation than there, for reasons to be discussed below. Let
3.6
be the Pauli matrices. We will need their representations in the adiabatic basis, given by
3.7
Above, note that . The usual algebraic relations of the Pauli matrices give
3.8
where **1** is the unit matrix. Moreover, the special relations
3.9
can be easily checked. Here and henceforth, primes denote derivatives with respect to *q*.

## Remark

We will use the Pauli matrices as a basis to represent the π_{n}. The basic idea is the same as in Betz & Teufel (2005*a*), but there the basis matrices *X*,*Y*,*Z* were chosen in an ad hoc manner. It happens that the matrices *X*,*Y*,*Z* from Betz & Teufel (2005*a*) are linked to **σ**_{x},**σ**_{y},**σ**_{z}, but unfortunately not in the most convenient way: we have *X*=i**σ**_{y}, *Y* =−**σ**_{z} and *Z*=−**σ**_{x}. This will lead to a serious clash of notation between the present paper and that of Betz & Teufel (2005*a*), but we feel that the formulation of the problem in the widely used Pauli matrices justifies this.

Note that **σ**_{z}(*q*)=*V* (*q*)/* ρ*(

*q*), and thus the adiabatic projection is given by 3.10 Indeed, owing to , and 3.11 The final equality above follows from the fact

**σ**

_{z}′(

*q*)=−

*θ*′(

*q*)

**σ**

_{x}(

*q*). In particular,

*F*

_{1}(

*p*,

*q*)=(i/2)

*p*

*θ*′(

*q*)

**σ**

_{x}(

*q*), and proposition 2.2 then gives which is just the adiabatic representation (1.5) in symbolic language.

We define the coefficients *x*_{n}, *y*_{n}, *z*_{n} and *w*_{n} through
3.12
and emphasize that these coefficients have swapped names with those given in Betz & Teufel (2005*a*). The pre-factor i in front of *y*_{n} will make for some more elegant formulas later, and slightly reduce the clash of notation with Betz & Teufel (2005*a*). For our derivation of the recursions for the coefficients *x*_{n} to *w*_{n}, we first need to treat the derivatives of *V* appearing in equation (3.1).

## Lemma 3.1

*We have**where a*_{n}(*q*) *and b*_{n}(*q*) *are given by the recursions*3.13

## Proof

Use the fact *V* =*ρ***σ**_{z} together with equation (3.9). ▪

The first step towards the coefficient recursion is as follows.

## Proposition 3.2

*Define x*_{n},*y*_{n},*z*_{n}*and w*_{n}*through equation (3.12). Then*,

## Proof

This is just a calculation. The former parts of each bracket stem from the first term of equation (3.1), which computes to
The latter terms of each bracket come from the various terms of the second part of equation (3.1); for even *k*, those are given by , and a direct calculation gives
for odd *k*, they are given by , with [*A*,*B*]_{+}=*A**B*+*B**A*. Noting that
and collecting coefficients give the result. ▪

It is remarkable that from this knowledge of *F*_{n+1} alone, we can derive a set of recursive differential equations that, together with zero boundary conditions at infinity, determine the coefficients *x*_{n} to *w*_{n}.

## Proposition 3.3

*The coefficients x*_{n}*to w*_{n}*defined in equation (3.12) are determined by the following recursive algebraic-differential equations: we have*3.14*Moreover*,
3.15*For n odd, we have*3.16*while for n even, we have*3.173.18*and*3.19

## Proof

For equation (3.14), note that
We now use equation (3.3) and proposition 3.2 in order to prove the recursive formulae: from equation (3.5), we can deduce that [*G*_{n+1},π_{n}]=0, and thus equation (3.9) implies that *G*_{n+1} is proportional to **σ**_{z} and **1**. Consequently, the parts of π_{n+1} proportional to **σ**_{x} and **σ**_{y} arise from the last term of equation (3.3) alone, and comparison with equation (3.9) and proposition 3.2 shows
and
Now, from equation (3.4), we can deduce that *F*_{n} is proportional to **σ**_{x} and **σ**_{y} only, and so proposition 3.2 immediately gives
and
Now equation (3.15) follows inductively, and after that equations (3.16)–(3.19) are immediate from the above equations, noting that all we did is to add some terms that are zero in order to get a more closed expression. ▪

We now cast proposition 2.2 into a more specific form, yielding our first main result.

## Theorem 3.4

*Define H*_{n}(*ε*,*p*,*q*) *as in proposition 2.2, x*_{n}(*p*,*q*) *and y*_{n}(*p*,*q*) *as in proposition 3.3 and put*3.20*Then*,

## Proof

From proposition 2.2, we see that we need to calculate *U*_{0}*F*_{n+1}*U*_{0}. By definition, *U*_{0}**σ**_{x}*U*_{0}=σ_{x}, and similarly for **σ**_{y}. As *F*_{n+1} is proportional to **σ**_{x} and **σ**_{y}, by comparing propositions 3.2 and 3.3, we obtain
Comparison with equation (3.6) and proposition 2.2 gives the result. ▪

In order to make use of theorem 3.4, we need to control the coefficients *x*_{n} and *z*_{n}. As they are polynomials (of order *n*) in *p*, it makes sense to consider coefficients. We put
3.21
with similar expressions for the other coefficients. Here, the index *m* in is an upper index rather than a power, and the choice of using for the highest power is a matter of convenience. Differentiating now gives
and thus the latter terms of expressions like equation (3.16) are of the form
3.22

## Proposition 3.5

*The coefficients**to**defined in equation (3.21) are determined by the following recursive algebraic-differential equations: we have*3.23*Moreover*,
*for n odd, while for n even, we have*

## Remark

From the above equations, it is obvious that for odd *m*. But more is true. By induction and using equation (3.15), we find that
3.24
for .

It appears to be a very hard problem to determine the asymptotic behaviour of the quantities etc. as , let alone prove it. Only in the case *m*=0, i.e. for the term of the highest order (in *p*), the asymptotics of , and are known, even rigorously. In that case, the sums on the right-hand sides in proposition 3.5 are empty, and we retain
After changing to the natural scale
3.25
these are (apart from a change of notation) just the recursions appearing in the time-adiabatic case, which have been solved in Betz & Teufel (2005*a*,*b*).

The relevant results in that case read as follows: we introduce the natural scale (3.25) and define for a given function *f*. Furthermore, we assume that
3.26
for some and τ_{c}>0, where has no singularities in , and only singularities of order smaller than one at ±iτ_{c}. As has first been observed in Berry & Lim (1993), equation (3.26) is the correct form of for a large class of models, including the generic ones. For further discussion, we refer to Berry & Lim (1993) and Betz & Teufel (2005*b*).

From equation (3.23), we now conclude
Let us write for the coefficient of belonging to *p*^{n}. Then by equation (3.20), we have
By the results of Betz & Teufel (2005*a*,*b*), there exists β>0 such that
3.27 is the universal pre-factor for the time-adiabatic transitions.

Unfortunately, the terms and appear not to constitute the leading-order contribution to *x*_{n}(*p*,*q*) and *y*_{n}(*p*,*q*) as . We expect this to be true in general, but it can be verified in the Landau–Zener case (Betz & Goddard in preparation), where we have for finite *m* and large *n*. We do not know the behaviour of when *m* is of the order of *n*; however, because numerical calculations clearly show that the latter case is where is maximal for fixed *n*, this is the regime that needs to be understood in order to do exponential asymptotics. On the other hand, for all models that we studied, the constant *c* above is very small. Therefore, if *n* is not too large, the *m*=0 term dominates. This is what will guide us in the next section, where we replace *x*_{n}(*p*,*q*) with without further justification in order to obtain our transition formulae. But let us remark that at least in the high-momentum regime (more precisely, for *p*=*ε*^{−β} with β>1/3), an adaption of the methods introduced in Betz & Teufel (2005*a*) leads to a proof that is indeed asymptotically dominant. The proof is not short and would distract from our main results, so it will be given elsewhere. The full asymptotic behaviour of the double sequence for finite *p* is, however, an intriguing and probably difficult open problem.

## 4. Transitions

We consider the Schrödinger equation in the *n*th superadiabatic representation,
4.1
where is the Weyl quantization of the symbol
Here, we split *H*_{n}(*p*,*q*), and in the same way also , into a diagonal part of order 1 and an off-diagonal ‘coupling’ part of order *ε*^{n+1}, and is as given in theorem 3.4. We assume that the potential *V* (*x*) approaches a constant matrix at spatial infinity quickly enough. This guarantees that all of the superadiabatic subspaces approach the adiabatic subspaces for large *x*. More precisely, we assume for the moment that the limits exist, that and that the limits are approached sufficiently fast. We write *ψ*_{+}(*t*,*x*) and *ψ*_{−}(*t*,*x*) for the first and second components of ** ψ**, respectively, and study solutions that are in the upper adiabatic subspace for . The first-order time-dependent perturbation theory determines

*ψ*

_{−}(

*t*,

*x*) up to errors of order

*ε*

^{n+1}. The perturbative solution of equation (4.1) with initial datum

*ψ*

_{+}(0,

*x*)=

*ψ*

_{+,0}(

*x*) and boundary condition is and 4.2 Here is as given in theorem 3.4, and

In view of formula (4.2), we need the Weyl quantization of . Given that is a polynomial in *p*, we will be interested in a general formula for the Weyl quantization of a symbol of the form *p*^{m}*g*(*q*). Moreover, we will later need the Fourier representation of the operator . To this end, we define the *scaled Fourier transform*4.3

## Lemma 4.1

*Let κ*(*p*,*q*)=*p*^{m}*g*(*q*) *be a semiclassical symbol for some*. *Then*,
4.4

## Proof

We use equation (2.4) with in order to get In the second line, we changed variables , and in the fourth . We now apply the scaled Fourier transform to both sides of the above equations and use 4.5 in order to obtain the result. ▪

## Remark

In position space, the Weyl quantization of κ(*p*,*q*)=*p*^{m}*g*(*q*) is given by
4.6
which can be seen from equation (2.4) using integration by parts with respect to *y* and equation (4.5).

We now give the momentum space version of equation (4.2) in a fairly explicit form. We write for the unitary propagator in the Fourier picture, and write
according to our results from §3. Combining equation (4.2) and lemma 4.1 now immediately shows that to leading order in *ε*,
4.7
Here, the operator *J*_{n+1} is given by , with
4.8

Although general, equation (4.7) is not very helpful when trying to calculate actual superadiabatic transition wave functions. To make further progress, we will make two simplifying assumptions. Firstly, we will drop all terms of with *m*>0. As discussed at the end of the previous section, this is necessary because we do not know the asymptotics of the for large *n* and *m*. We re-emphasize that while is not, in general, asymptotically dominant, it is expected to accurately describe the high-momentum regime, and our numerical results confirm that it is an excellent approximation for a wide range of parameters. The second simplification is to assume constant energy levels. This means that in equation (1.1), we will take
4.9
with *θ*_{r} having no singular points of order greater or equal to one in the strip . Constant * ρ* has the effect of trivializing at the same time the propagator in Fourier space,
and the transformation (3.25) to the natural time scale, with τ(

*q*)=2

*δ*

*q*. Then, equation (3.26) becomes with τ

_{c}=2

*δ*

*q*

_{c}, and equation (3.27) yields Using the residual theorem, the scaled Fourier transformation in

*q*of is given by By equation (4.8), we find We now insert this into equation (4.7). With the above conventions, and our calculations so far, equation (4.7) reads 4.10 For the moment, we are interested in for

*t*≫0. The integral in equation (4.10) converges as , and in the limit can be calculated explicitly. Consider the general integral where

*I*

_{±}is the integrand integrated over {±η>0}. Setting gives and , and so where χ

_{J}is the characteristic function on .

As , where is the inverse Fourier transform of *g*, we obtain
Applying these calculations to equation (4.10) reveals that, for large positive *t*, we have
4.11
The two terms in the square bracket in equation (4.11) are clearly connected to positive and negative incoming momenta, respectively. The second line will be negligible if either *k*<0 or is concentrated on the negative half-axis, while the third line will be negligible if *k*>0 or if is concentrated on the positive half-axis. This shows the intuitively obvious fact that the transmitted wave packet will travel in the same direction as the incoming one. It also shows that we can replace with without changing the leading-order result. We further streamline equation (4.11) by replacing α(*n*+1) by its asymptotic value , and by introducing
As asymptotically all the superadiabatic subspaces agree, equation (4.11) actually gives the asymptotic adiabatic transition.

As a result, after the transition, the wave function in the initially unoccupied adiabatic subspace is given by 4.12 A few comments are in order.

(i) The occurrence of the indicator function χ_{{k2>4δ}} can be interpreted in terms of energy conservation. Any part of the wave packet that makes the transition obtains, in addition to its kinetic energy, the potential energy difference between the electronic energy levels. Thus, there cannot be any kinetic energy *k*^{2}/2 smaller than 2*δ*. Recall also that these transitions are radiation-less: instead of being radiated away from the molecule in the form of a photon, the energy is transferred into kinetic energy of the nuclei.

(ii) We can read off a momentum shift from equation (4.11). We assume that is a semiclassical wave function, and write . Let us assume for convenience that the absolute minimum *k*_{*} of *M*(*k*) is on the positive real line. We find
4.13
Purely by energy conservation, one would expect the transition wave packet to be maximal when *v*=*k*_{*}. However, as is decreasing, the minimum of is shifted to the right. One can quantify this shift when *M* is given explicitly.

(iii) The last point also shows that the term χ_{{k2>4δ}} in equation (4.11) is of little consequence in practice. As *k*^{2}>4*δ* is equivalent to *v*>0, and since only a small region around its maximum matters for the semiclassical wave function , we can safely leave out the factor χ_{{k2>4δ}} in equation (4.11) without changing the leading-order result.

(iv) Equation (4.11) depends on *n* only through the convergent pre-factor α(*n*). In particular, we do not need to know the value of the optimal *n* in order to obtain the correct leading-order transitions. This suggests that, as in the time-adiabatic problem, one could obtain equation (4.11) by naive perturbation theory in the adiabatic basis; the renormalization of the pre-factor is then the one known from the time-adiabatic theory. However, this situation is special to the case where the potential energy matrix *V* has no trace, and fails in more general cases (Betz & Goddard in preparation).

(v) Let us compute the transition rate from equation (4.12), in the limit of large momentum and small momentum uncertainty. We choose When *C* is large, the minimum in equation (4.13) is taken very close to *v*=*p*_{0}, implying that . The value of the exponent at the maximum is then given by which is the transition rate for momentum *p*_{0}. For large *p*_{0}, we have , so that the transition probability (up to a pre-factor) in this regime is given by The latter is precisely the Landau–Zener transition probability for the parameters chosen in equation (4.9), cf. e.g. Betz & Teufel (2006), where one has to replace *ε* by *ε**p*_{0} throughout.

Equation (4.12) only uses local information: the parameters *q*_{c}, *δ* and γ are determined by the derivatives of the potential at the coupling point, and the incoming wave function is only needed at the crossing time. This immediately suggests an algorithm for computing non-adiabatic transitions. As in surface-hopping models, the wave packet evolves on the initial adiabatic surface under the uncoupled BO dynamics, until a local extremum in the coupling functions *θ*′(*q*_{0}(*t*)) is detected, where *q*_{0}(*t*) is the centre of the wave packet at time *t*. One then extracts the parameters *q*_{c}, *δ* and γ from the local shape of the potential energy matrix, starts a wave packed according to equation (4.8) on the initially unoccupied energy surface and evolves it under the uncoupled BO approximation of that surface. The same algorithm is expected to work in the general case of non-constant energy surfaces, but equation (4.12) needs to be modified to account for the changed dynamics in the transition region. This is not trivial and is currently under investigation.

We now show that the results of the above algorithm are in excellent agreement with highly precise numerical solutions of the full two-band Schrödinger equation. For solving the latter, we use standard methods, including a symmetric Strang splitting. We write ϕ_{−}(*q*,*t*) for the projection of the numerical solution onto the lower adiabatic eigenspace, and *ψ*_{−} the result of applying equation (4.12) and proceeding as indicated above. We compare our results in the Fourier representation, calculating the inverse *ε*-Fourier transform of ϕ_{−} with a standard fast Fourier transform (FFT). The final time *t* is chosen so that ∥ϕ_{−}(.,*t*)∥_{2} is constant under further time evolution in the exact calculation.

For our potential, we choose
4.14
in equation (1.1). This gives , with singularities closest to the real line at ±i*q*_{c}=±i(π/2α), the residue at which γ=−*c*/2α. We take *c*=−π/3, α=π/2, giving *q*_{c}=1 and γ=1/3. The choice of *θ* over the case (which would make *θ*_{r}=0 in equation (3.26)) is made to increase the rate at which the potential becomes flat. This reduces the necessary computation time for the numerical solution. If anything, we would expect the asymptotic results for to be better.

Our first choice for the wave function in the initially occupied band is a Gaussian wave packet. At time *t*=0, it is given by
4.15
As the crossing point is at *x*=0, by our choice of potential, *ψ*_{+,0} is sitting right at the middle of the crossing region. As the eigenvalues are constant, it may be evolved backwards to *t*_{0}<0 exactly on the upper level
4.16
We use equation (4.16) along with as initial conditions for our numerical solution, and take *t*_{0} sufficiently negative in order for the wave packet to be well away from the crossing region.

We find that equation (4.12) is in excellent agreement with the numerical solution for a wide range of *ε*, ranging from as large as 1/5 to 1/50, at which point a further reduction in *ε* makes the numerically exact calculations very time consuming. Figure 1 shows the relative error in the *L*^{2} norms between the numerical calculation and equation (4.12), i.e. ∥*ψ*_{−}−ϕ_{−}∥_{2}/∥ϕ_{−}∥_{2}. In each case, the step size in the numerical calculation was reduced until the difference between two subsequent numerical solutions was at least one order of magnitude smaller than the error to the solution obtained from equation (4.12).

Figure 1*a* also gives evidence that, as discussed before, equation (4.12) is not asymptotically correct for fixed *p*: after an initial increase in accuracy owing to the decrease in *ε*, the relative error becomes larger again as *ε* decreases further. This is unlikely to affect the practical usefulness of equation (4.12), which becomes clear when we consider orders of magnitude: for *p*_{0}=2 and *ε*=1/50, we have ∥*ψ*_{−}∥_{2}≈6×10^{−10}. While it is not clear whether quantities of such size are relevant in practice, the relative error is still excellent at about 4×10^{−3}, albeit deteriorating. On the other hand, for *p*_{0}=2 and *ε*=1/10, we have ∥*ψ*_{−}∥_{2}≈0.014, which is certainly of a physically measurable size, and the relative error is still of the order 10^{−3}. Finally, for *ε*=1/5 and *p*_{0}=2, we have ∥*ψ*_{−}∥_{2}≈0.11, and the relative error is still below 0.03. As is often the case in asymptotic formulae, the actual error is much better than could be expected from the *a priori* error estimates. This is particularly obvious in figure 1*b*, where the relative error initially decreases like e^{−c/ε} before saturating, while theory only predicts a decrease. Orders of magnitude in this case range from ∥*ψ*_{−}∥_{2}≈0.138 for *ε*=1/10, with relative error below 0.025, to ∥*ψ*_{−}∥_{2}≈6×10^{−5} for *ε*=1/50, with relative error below 2×10^{−5}.

Our second numerical example is an incoming wave packet that is strongly non-Gaussian. We choose where *Z* normalizes the *L*^{2} norm of . The potential *V* is as before, and we choose *p*_{0}=5 and *ε*=1/50. Figure 2 shows the absolute value of the transmitted wave packet in the scattering region in Fourier representation. The relative error in the *L*^{2}-norm in this case is around 5×10^{−5}, similar to the one in the Gaussian case. In particular, the pointwise relative error is extremely small in the region where is concentrated; thus, while figure 2 shows the result of applying equation (4.12), the plot showing the numerical calculation would be indistinguishable from the one given. The momentum shift relative to the energy conservation value is clearly visible, and the transmitted wave packet is strongly non-Gaussian, with a surplus of high-momentum components having made the transition.

Let us finally turn to the time development of the transition. From the time-adiabatic theory, we know that in the optimal superadiabatic representation, the transmitted part of the wave function builds up monotonically and has the shape of an error function as a function of time. We are now going to show that a similar property holds for the BO transitions after suitable modifications.

We start with theorem 4.7, which gives a formula for transitions in the *n*th superadiabatic representation. Our first task is to find a precise meaning to optimality of a superadiabatic representation. One natural idea would be to require that in the optimal superadiabatic representation the map increases monotonously during the transition. But given the complexity of formulas (4.7) or (4.10), this condition is difficult to check. Instead, we choose another condition that, as we will argue, should be equivalent. The basic idea is that through equation (4.7), is given by an integral where the integrand is expected to be both highly oscillatory and sharply peaked in *s* and *k*. The main contribution then occurs at points where the integrand has either a stationary phase or a maximal absolute value. The locations of both of these depend on *n*, and if we want to have any chance of seeing a ‘nice’ transition history, we have to choose (if possible) *n* such that these locations coincide. Thus, we define as follows.

## Definition 4.2

*Let*, *as given in equation* (4.7). *We write*. *Let* (*s*_{*}(*n*),*k*_{*}(*n*)) *be the location of the minimum of X on the real line for given n. We say that n is an optimal superadiabatic representation if* ∂_{s}*Y* (*s*_{*}(*n*),*k*_{*}(*n*))≈0.

The ≈ sign in the above definition means that, in principle, *n* has to be an integer, and thus equality might not hold anywhere. On the one hand, this problem should be less severe when *ε* is extremely small, and on the other hand, usually *J*_{n+1} is given by an explicit formula, in which case we can interpolate and take *n* to be real. In that case, we will usually be able to fulfil ∂_{s}*Y* (*s*_{*}(*n*),*k*_{*}(*n*))=0 exactly.

A simple calculation shows that the above definition indeed gives ‘optimal’ transition histories. We write etc., and , . Expanding the exponent of equation (4.7) around *k*_{*} and *s*_{*} gives
At the maximum of the transmitted wave function, i.e. for , we have
We see that the transmitted wave function, when adjusted for the propagation in the lower band, has the shape of an error function at its maximum *k*_{*}. The only unusual feature is that this error function is actually evaluated along a ‘diagonal’ in the complex plane, rather than on the real line. Nevertheless, the resulting shape will be close to monotone unless *Y*_{ss} is much larger than *X*_{ss}, in which case there are some oscillations around *t*=0. We do not know whether this case is likely to happen in practice. For , the error-function behaviour deteriorates, but as is a semiclassical wave function, we are only interested in , in which case the behaviour is still similar to the one at the maximum.

We now show that optimal superadiabatic representations exist in particular cases. We pick again the situation of constant eigenvalues, and rewrite the integral in equation (4.10) as
4.17
Here, e^{−M/ε} is the combined modulus of the transition kernel originating from and the wave function, and ϕ is the phase of the wave function, hence depending only on η. We can treat more general forms of the coupling function than the one leading to equation (4.10), including the full coupling function with all powers of *p* included. The only requirement is that, to leading order, κ should be symmetric, so that its Fourier transform is real. This is generically true: high derivatives of the pair of first-order poles in the complex plane determines the shape of κ, and these are either symmetric or antisymmetric, giving either purely real or purely imaginary Fourier transforms.

Given equation (4.17), we now let η_{*}=η_{*}(*n*) and *k*_{*}=*k*_{*}(*n*) be the place where *M*(*k*,η) is minimal. We expand *M* to second order around (*k*_{*},η_{*}). In order for the phase to be also stationary, we need ∂_{η}ϕ(η_{*})=*s*η_{*}, which is the equation determining the transition time *s*_{*}. We should keep in mind that any further explicit calculations will only make sense for *s* close to *s*_{*}. Now, we expand the phase to second order around η_{*} as well, and compute the resulting Gaussian integral in η. The result is
4.18
Again, we use the notation , etc. As before, we concentrate on the case . Let us assume that we can find *n* such that . It is then easy to check that for the remaining integrand both real and imaginary parts of the exponent are stationary at *s*=*s*_{*} and that *s*_{*} is a maximum of the real part. Thus, for constant eigenvalues, we can find an optimal superadiabatic representation if we can solve the equations
4.19
simultaneously, and if the resulting pair (*k*,η) is a minimum of *M*. As the above equations also depend on *n*, we are solving a system of three equations with three free variables, which means that we can hope for a solution. We specialize further to an example where we can show equation (4.19) to be solvable. We choose the upper-band wave function as a Gaussian wave packet with momentum *p*_{0},
4.20
Choosing real-valued amounts to choosing the packet to be at *x*=0 at time *t*=0; putting the avoided crossing at *x*=0 in addition ensures that the transition occurs at time *s*_{*}=0. Now, using the integrand in equation (4.10), we get, to leading order,
By the third equation in (4.19), we know *k*>η, which removes the absolute value above. Taking derivatives, we get 0=2*n**ε*η/(*k*^{2}−η^{2})−*q*_{c}+(η−*p*_{0})/σ^{2} and 0=−2*n**ε**k*/(*k*^{2}−η^{2})+*q*_{c}. Together with the equation *k*^{2}−η^{2}=4*δ*, these lead to
4.21
The first two equations are independent of *n* and *ε*, and determine η_{*},*k*_{*}. This is a special feature of the constant eigenvalue situation, and should not be expected in general. The third equation determines the optimal superadiabatic representation. Note that *n* is connected to the optimal superadiabatic *n* for the time-adiabatic situation: there, we have *n*_{ta}=2*δ**q*_{c}/*ε*, with the ‘momentum’ (i.e. speed on the time axis) normalized to one. But the tricky point that cannot be easily guessed is which value of *k* to pick in the formula for *n*_{*}; the naive guess of using the incoming momentum *p*_{0} would be totally wrong. The more sophisticated guess of using the mean momentum at the crossing point, (η_{*}+*k*_{*})/2, would be closer for small *δ* because the true value fulfils ; but it would still be far off for finite *δ*, which are of main interest here.

We close by comparing the effective formula (4.18) with numerical results for the transition histories in various superadiabatic representations. In order to obtain the optimal superadiabatic representation at both low and close to integer values of *n*, we have to choose the parameters somewhat carefully. Our choices are the potential (4.14) with *c*=−π/3, α=2π/5 and *δ*=3/32. This gives *q*_{c}=5/4 and γ=5/12. We take *ε*=0.02923, and in the incoming wave function (4.20), we take σ^{2}=2 and *p*_{0}=2.5. Solving equation (4.21) then yields η_{*}≈2.57, *k*_{*}≈2.64, and *n*_{*}≈3.04. Thus, we expect to obtain the optimal superadiabatic representation when *n*=3, and this is clearly confirmed by numerical simulations. Figure 3 shows the *L*^{2} norm of the transmitted wave function (calculated by a numerical solution of the Schrödinger equation), in the adiabatic and all superadiabatic representations up to *n*=5, as a function of time *t*. The dashed line is the prediction of formula (4.18) with η_{*} and *k*_{*} inserted from above. We see that indeed the optimal superadiabatic representation is at *n*=3, while below *n*=1 and above *n*=4 oscillations grow. The reader should note the similarity with the plots shown in Lim & Berry (1991). This is another confirmation that the time-adiabatic approximation is appropriate to qualitatively understand the mechanism of non-adiabatic transitions, as far as the population of the lower level is concerned. On the other hand, in order to obtain quantitatively correct results, or more detailed information about the transmitted wave packet (e.g. momentum spread of phase shift), the full quantum mechanical treatment, as given in the present work, is indispensable.

## Acknowledgements

V.B. is supported by the EPSRC fellowship EP/D07181X/1.

## Footnotes

- Received June 29, 2009.
- Accepted August 4, 2009.

- © 2009 The Royal Society