Multiple Superadiabatic Transitions and LandauZener Formulas
Abstract
We consider nonadiabatic systems in which the classical BornOppenheimer approximation breaks down. We present a general theory that accurately captures the full transmitted wavepacket after multiple transitions through either a single or distinct avoided crossings, including phase information and associated interference effects. Under suitable approximations, we recover both the celebrated LandauZener formula and standard surfacehopping algorithms. Our algorithm shows excellent agreement with the full quantum dynamics for a range of avoided crossing systems, and can also be applied to single full crossings with similar accuracy.
pacs:
I Introduction
The BornOppenheimer approximation (BOA) Born and Oppenheimer (1927) is one of the most widely used methods used to study the quantum dynamics of molecules. Intuitively, it is motivated by the fact that the electrons are much lighter, and therefore much faster, than the nuclei, and hence rapidly adjust their positions with respect to those of the nuclei. This scale separation allows, in many cases, for the electronic and nuclear dynamics to be decoupled. In particular, if the electrons start in a certain bound state, for a fixed set of nuclei positions, then they should remain in this bound state even though the nuclei are slowly moving. Hence the nuclear dynamics can be determined by considering their motion on only one (electronic) potential energy surface.
However, there are interesting situations in which the BOA breaks down Domcke, Yarkony, and Horst (2004); Domcke, Yarkony, and Köppel (2011); Nakamura (2012); Tully (2012). For example, in many photochemical processes the nuclear motion cannot be restricted to a single potential energy surface because, for some nuclear configurations, two such surfaces become close, or even cross. In the former case, known as an avoided crossing, the BOA is still valid to leading order (in the small parameter , which is the square root of the ratio of the electronic and nuclear masses), but the remaining corrections are of fundamental interest and, in fact, determine the associated chemistry. In the latter case, which generally takes the form of conical intersections, the BOA breaks down completely.
Here we are primarily interested in cases where the transmitted wavepacket is (exponentially) small Hagedorn and Joye (2001, 2005); Martinez and Sordoni (2002), for example when there is an avoided crossing, or when the wavepacket does not pass directly over the conical intersection. Such regimes are, in some sense, generic, as avoided crossings are generic in 1D Von Neumann and Wigner (1929), and in higher dimensions the probability of an arbitrary wavepacket exactly hitting a conical intersection is vanishingly small Tully (2012). In particular, we consider cases where the wavepacket passes through multiple avoided crossings, or repeatedly through the same crossing. In such cases the transmitted wavepackets can interfere, and thus it is necessary to understand their phases. This suggests that a full quantum mechanical treatment of the problem is required. However, in even moderate dimensions, such treatments are numerically intractable, especially for multiple, coupled electronic potential surfaces.
In order to overcome this, a range of coupled quantumclassical and semiclassical methods have been developed. These include the multiplespawning wavepacket method BenNun and Martínez (1998); BenNun, Quenneville, and Martínez (2000); Virshup, Chen, and Martínez (2012), the frozen Gaussian wavepacket method Heller (1991), Ehrenfest dynamics McLachlan (1964); Meyer and Miller (1979); Sawada, Nitzan, and Metiu (1985), and the semiclassical initial value representation Miller (1970); Kreek and Marcus (1974); Miller (2001). The main advantage of such schemes is the significantly reduced computational cost. The main disadvantage, at least with respect to the problem at hand, is the lack of phase information from almost all such schemes. Along with those mentioned above, one of the most widelyused quantumclassical approaches is surface hopping Tully and Preston (1971); Miller and George (1972); Stine and Muckerman (1976); Kuntz, Kendrick, and Whitton (1979); Tully (1990); HammesSchiffer and Tully (1994); Müller and Stock (1997); Fabiano, Groenhof, and Thiel (2008); Fermanian Kammerer and Lasser (2008); Lasser and Swart (2008); Belyaev, Lasser, and Trigila (2014); Belyaev et al. (2015), in which particles are evolved under classical dynamics on a single surface and can ‘hop’ to other surfaces with a specified probability. Perhaps the most common approach is to only allow hops at points in the trajectory where the gap between energy surfaces has a local minimum (i.e. at an avoided crossing), and the probability of the hop is given by a LandauZener (LZ) formula Zener (1932); Landau (1965). Such methods give good results for a single transition, especially when the transmitted wavepacket is reasonably large, but fail completely when multiple transitions are involved, due to the complete lack of phase information Fermanian Kammerer and Lasser (2017). We note here that there is at least one such scheme Chai et al. (2015) that does aim to retain the phase information, but this is limited to small gaps between the potential energy surfaces, which in turn leads to large transmitted wavepackets. The same restriction is true for other mathematical approaches that lead to explicit formulae for the transmitted wavepacket; see e.g. Ref. Hagedorn and Joye (2005). It has been shown that, if the gap scales with , then the transitions are of order one and dominated by the LandauZener factor Hagedorn (1994); Hagedorn and Joye (1998).
An alternative approach, inspired by the work of Berry on superadiabatic representations Berry (1990); Berry and Lim (1993), considers the full quantum mechanical wavepacket. These results, which are restricted to the semiclassical regime where the nuclei move classically, were later made rigorous Hagedorn and Joye (2004); Betz and Teufel (2005, 2005). It was later shown that, through the use of such superadiabatic representations (which are generalisations of the wellknown adiabatic representation), it is possible to derive a formula for the transmitted wavepacket, including phase, at an avoided crossing Betz, Goddard, and Teufel (2009); Betz and Goddard (2009, 2011); Betz, Goddard, and Manthe (2016). The associated algorithm requires only the quantum evolution of wavepackets on single energy surfaces. Whilst this is still computationally demanding if one wants to solve the full Schrödinger equation, there are approximate methods, such as Hagedorn wavepackets Hagedorn (1981, 1994); Lubich (2008) or standard quantum chemistry techniques such as MCTDH HuarteLarrañaga and Manthe (2009), which make small relative errors and are computationally much more tractable. The associated algorithm has so far been applied to single transitions through avoided crossings Betz, Goddard, and Teufel (2009); Betz and Goddard (2009, 2011), and to multiple transitions of a single crossing in the case of the photodissociation of NaI Betz, Goddard, and Manthe (2016). The main goals here are to extend the methodology to multiple transitions through different avoided crossings and to systematically study the effects of making various approximations that lead to a LZlike transition probability. We will also demonstrate that, although not designed to tackle such problems, the methodology can be successfully applied to single transitions of full crossings.
We present an algorithm that has a number of advantages. We have already mentioned: (i) Preservation of phase information, which allows the accurate study of interference effects; (ii) Only evolution on a single surface is required, which significantly reduces the computational cost when compared to a fullycoupled system, whilst also allowing the use of stateofthe art numerical schemes. The main other benefits are: (iii) Only the adiabatic surfaces (which are the most commonly obtained surfaces from quantum chemistry calculations) are required, in particular there is no need for a diabatization scheme, or the determination of the adiabatic coupling elements; (iv) Such surfaces are only required locally, and thus can be computed onthefly; (v) The transmitted wavepacket is created instantaneously, and hence there is no reliance on complicated numerical cancellations of highlyoscillatory wavepackets, which are generally present in the adiabatic representation; (vi) The methodology is easily extended to multiple adiabatic surfaces; (vii) The derived formula is accurate for a wide range of potential energy gaps and small parameters , and for any semiclassical wavepacket, i.e. one of typical width or order .
There are, of course, also some disadvantages when compared to the more widelyused schemes: (i) In order to capture the phase information, the onelevel dynamics must retain at least some of their quantum nature, and this is inherently more computationally demanding than the analogous classical dynamics; (ii) In the full formalism, it is necessary to be able to extend the potential surface into the complex plane, at least in the region of an avoided crossing. This is essential to be able to accurately compute the transition probabilities. However, in some regimes, for example when the LZ formula is accurate, we can bypass this requirement; (iii) The scheme is, in principle, restricted to wavepackets that are semiclassical near the avoided crossing. However, due to the linearity of the Schrödinger equation, and as demonstrated in Betz, Goddard, and Manthe (2016), it is possible to ‘slice’ the wavepacket at the crossing. However, this may be more problematic in higher dimensions; (iv) As it stands, the method is restricted to 1D. However, we have successfully extended it to higher dimensions through a slicing procedure [[REF 2D PAPER]].
To outline our approach, we will first review the standard model for nonadiabatic transitions (Section II) and avoided crossings (Section III). We will then, in Section IV, give a brief overview of existing surface hopping models and LZ formulas. We then outline the superadiabatic approach and give the resulting formula in Section V, before describing the associated algorithm in Section VI. In Section VII we systematically investigating its accuracy, and the effects of replacing the true transition probability by two LZlike approximations. Finally, in Section VIII, we summarize our results and discuss some open problems.
Ii The Model
The Schrödinger equation governing the quantum dynamics of a molecular system can be written as
where and are the nuclear and electronic positions, respectively, and the Hamiltonian is given by
Here the first two terms are the kinetic energies of the nuclei and electrons with masses and , respectively. Note that the masses of the nuclei may all be chosen to be the same by a rescaling of the nuclear coordinates. The potentials and denote the nuclear and electronic Coulomb repulsions, respectively, whilst is the attraction between the nuclei and electrons.
We now change to atomic units () and define and the electronic Hamiltonian for fixed nuclear positions ,
Suppose that are two eigenvalues of the electronic Hamiltonian (i.e. two adiabatic potential energy surfaces) of multiplicity one and wellseparated from the rest of the electronic spectrum. Then, from BornOppenheimer theory Hagedorn (1980); Spohn and Teufel (2001), the effective nuclear Schrödinger equation is
(1) 
where is a matrix with eigenvalues , i.e. a diabatic matrix. In general, is symmetric and has the form
For notational convenience, and to connect back to previous work Betz, Goddard, and Teufel (2009); Betz and Goddard (2009, 2011); Betz, Goddard, and Manthe (2016), we find it useful to define
and so
It is easy to see that the adiabatic surfaces are then given by and so is half the energy gap between the two surfaces.
Iii Avoided Crossings
In the adiabatic representation, an explicit unitary transformation is applied to the system such that the electronic Hamiltonian is diagonal at each choice of . Transitions between the adiabatic surfaces are then governed by the kinetic energy term, which introduces offdiagonal coupling elements, giving (for a onedimensional (1D) system) a Hamiltonian of the form
(2) 
Here is an explicit ‘kinetic coupling’ function and we have grouped the terms such that it is more obvious that the offdiagonal elements of the above matrix are of order . This can be seen from the fact that wavepackets typically oscillate with frequency (see Section III), so is of order one. Hence we see that, naïvely, the transitions are of order globally in time. However, as discussed previously, the transitions are exponentially small in away from the avoided crossings.
Typically, when the adiabatic potentials are wellseparated, the coupling elements are small and then two levels may be treated separately via the BornOppenheimer approximation. However, if the adiabatic surfaces become close, but do not cross, the coupling terms typically become large (but do not diverge). Such nuclear configurations are known as avoided crossings. As a result of the large coupling elements, a small, but not negligible, part of the nuclear wavepacket is transferred between the adiabatic surfaces.
Suppose, for clarity of exposition, that the wavepacket initially occupies the upper adiabatic level. The aim of this work is to determine the transmitted wavepacket (on the lower adiabatic level) well away from the crossing (in the scattering regime). Whilst one can, in principle, compute this by a standard numerical solution of the Schrödinger equatiion, there are a number of challenges that prevent this from being a realistic option for most systems of interest:

In order to compute the dynamics, one needs an accurate representation of the potential energy surfaces. Typically the adiabatic surfaces are calculated using quantum chemistry methods, such as Density Functional Theory, but it is computationally expensive to determine such surfaces, especially when the number of degrees of freedom (dimension of ) is large. In such cases, it is desirable to design methods that can utilise onthefly surfaces, determined only locally. Additionally, practical methods for determining surfaces for excited states are still in their infancy, and one also needs to determine the offdiagonal coupling elements. Finally, we note that diabatic representations are not unique, and those obtained in two and multiplelevel cases may differ significantly Belyaev et al. (2010).

The wavepackets we wish to compute are highly oscillatory, typically oscillating with frequency of order in space. This can be seen by comparing the kinetic and potential terms in (1). When using a standard numerical scheme, such as Strang splitting, correctly resolving such oscillations requires very fine grids in both position and momentum space. The curse of dimensionality (for points in dimensions, one requires points) results in such approaches being impractical for all but very small dimensional systems.

Away from avoided crossings, the transmitted wavepacket is typically exponentially small in both the gap size and . This can be seen from the formula (5) in Section V.2 or the standard LZ transition probabilities (3) and (4), where . In contrast, globally in time, the transitions in the adiabatic representation are of order , which we have already seen from the Hamiltonian 2. The necessary cancellations in the transmitted wavepacket occur through Stückelberg oscillations. See Figure 1 for an example. There are two challenges here. The first is to correctly resolve these cancellations, which can require very small time steps. The second is the more general challenge of computing an exponentially small quantity; any absolute errors in the numerical scheme must also be exponentially small or they will overwhelm the desired results.
Iv Existing Approaches and LandauZener
In this section we discuss some existing approaches to calculate the transition probability or the transmitted wavepacket.
iv.1 Surface Hopping Algorithms
Here we present a brief overview of surface hopping methods, which are one of the most successful approaches for simulating nonadiabatic dynamics. Surface hopping is a mixed quantumclassical approach, where particles are transported classically on the adiabatic surfaces and hop between them under certain conditions, which simulates the quantum effects. A general surface hopping algorithm consists of four steps:

Sampling of the initial condition.

Classical evolution via , .

Surface hopping.

Computation of observables.
There are many such schemes, both deterministic and probabilistic and we refer to Tully and Preston (1971); Miller and George (1972); Stine and Muckerman (1976); Kuntz, Kendrick, and Whitton (1979); Tully (1990); HammesSchiffer and Tully (1994); Müller and Stock (1997); Fabiano, Groenhof, and Thiel (2008); Fermanian Kammerer and Lasser (2008); Lasser and Swart (2008); Belyaev, Lasser, and Trigila (2014); Belyaev et al. (2015) for further details.
Of particular interest here is the surface hopping step. Typically this is performed when the gap between the two adiabatic surfaces is minimal along a classical trajectory. Whenever such a trajectory reaches a local minimum, a transition to the other surface is performed with a certain probability, usually derived from a simplified quantum mechanical model. The standard approach is to use a LZ formula, which we describe in the next section. The choice of this hopping probability is the main distinguishing feature of different surface hopping models.
The principal advantage of surface hopping algorithms is their simplicity. Due to their use of classical dynamics, which only require local properties of the potential energy surfaces, the methods can be applied in relatively high dimensions, using onthefly surfaces. As mentioned previously, such high dimensional systems are beyond the reach of full quantum mechanical methods.
The principal disadvantage is that they lose all phase information, and so cannot treat systems in which interference effects are important, or determine observables in which the relative phase of the wavepackets on the adiabatic surfaces is required Belyaev, Lasser, and Trigila (2014). Additionally, they are accurate only when the specified hopping probability is accurate; we will investigate this in Section VII.
iv.2 The LandauZener Formula
As mentioned previously, in order to compute the transition probability, it is common to use a LZ formula. Whilst the LZ model provides a simple formula for the transition probability, it is generally obtained from a onedimensional, twolevel model in the diabatic representation. However, practical applications occur in multiple dimensions and the potential energy surfaces are usually calculated in the adiabatic representation. There are a number of formulations of the LZ probability, including the extension to multiple dimensions in the diabatic formalism Fermanian Kammerer and Lasser (2008), and versions which only require knowledge of the adiabatic potentials Zhu, Teranishi, and Nakamura (2001); Belyaev and Lebedev (2011). Here we restrict ourselves to two such formalisms, the first is a diabatic representation. which requires knowledge of the diabatic matrix elements, whilst the second is an adiabatic representation, which only requires the gap between the adiabatic potentials. From now on, we consider only 1D systems; see Section VIII for some discussion of progress in higher dimensions.
Consider a classical particle with trajectory in phase space. Denote the position where attains a minimum by , and the momentum of the particle at the corresponding time by . Then the diabatic LZ transition probability is given by
(3) 
The corresponding adiabatic transition probability is given by
If one has knowledge of the diabatic matrix elements, and hence and , this can be rewritten as
(4) 
Note that in the corresponding multidimensional formula Belyaev, Lasser, and Trigila (2014) there is an additional term, which in 1D would be . However, since an avoided crossing is defined as a minimum of , and , this term is zero in 1D.
V Superadiabatic Representations and the Formula
In this section we will briefly review the ideas behind the use of superadiabatic representations to compute the transmitted wavepacket and refer the reader to the cited works for more details. We will then present a generalisation of a previouslyderived formula, which is applicable to 1D avoided crossings not centred at the origin.
v.1 Superadiabatic Representations
Superadiabatic representations were first introduced by Berry Berry (1990); Berry and Lim (1993), under the additional approximation that the nuclei move classically. More recently this has been extended to the full BOA Betz, Goddard, and Teufel (2009); Betz and Goddard (2009, 2011); Betz, Goddard, and Manthe (2016). As suggested by the name, superadiabatic representations are refinements of the adiabatic representation, which we described in Section refS:avoidedCrossings. In the adiabatic representation, transitions can be very complicated, as demonstrated by the population on the lower level during a typical transition, see Figure 1. This reliance on large cancellations to leave an exponentially small wavepacket suggests that the adiabatic representation may not be the ideal frame of reference in which to study transitions at avoided crossings.
Superadiabatic representations improve on the adiabatic one by simplifying the dynamics near an avoided crossing, at the expense of introducing computational complexities. The superadiabatic representations can be enumerated, and, initially, moving to successively higher superadiabatic representations reduces the spurious oscillations in the dynamics until the transmitted population builds up monotonically as the wavepacket travels through the avoided crossing. This is known as the optimal superadiabatic representation. However, moving to even higher representations results in the spurious oscillations returning. Previous results give a reliable method to determine the optimal superadiabatic representation Betz, Goddard, and Teufel (2009); Betz and Goddard (2011). However, computing the unitary operators for this representation is highly challenging, and performing the numerical computations in such a representation is similarly difficult.
The main benefit of superadiabatic representations for our purposes is that they allow the derivation of an explicit formula for the transmitted wavepacket in the optimal superadiabatic representation, without requiring the associated unitary matrix. By general theory Teufel (2003), all of the superadiabatic representations agree with the adiabatic one away from any avoided crossing. This leads to a simple algorithm to compute the transition through an avoided crossing in the adiabatic representation, as described in Section VI.
v.2 The Formula
Following Ref. Berry and Lim (1993), it is useful to introduce a nonlinear rescaling in which the adiabatic coupling elements obtain a universal form, known as the natural scale,
where is the position of the avoided crossing. We now extend and into the complex plane and, by the theory of Stokes lines Joye, Mileti, and Pfister (1991), the analytic continuation of has a pair of complex conjugate zeros, close to , at and . We define
Let be the incoming wavepacket on the corresponding adiabatic surface at time when the centre of mass coincides with an avoided crossing at . Then, for , the transmitted wavepacket on the other adiabatic surface can be approximated by
where are the BOA Hamiltonians for the two levels and is a wavepacket instantaneously created at time , which is more easily expressed in Fourier space via
(5) 
Here
is the Heaviside function and the Fourier transform needs to be performed under the correct scaling:
We note that the principle difference from previous presentations of the formula is the final exponential factor involving , the position of the avoided crossing. In previous work, this position has been taken to be zero, in which case the factor is simply 1. The new term arises from the approximation (which is a simple generalisation of the calculation for (Betz and Goddard, 2011, p. 2258)).
v.3 Analysis of the Formula
We now present a brief analysis of the formula in 5, which allows us to connect to surface hopping approaches, as well as to LZ formulas.
Firstly we note that the formula involves the same momentum adjustment that is phenomenologically introduced in surface hopping algorithms. We note that is precisely the classical incoming momentum required to give outgoing momentum when moving down/up, respectively, a potential energy gap of and requiring (classical) energy conservation. Relatedly, when passing from the upper to the lower level, the Heaviside function ensures that the transmitted wavepacket has (absolute) momentum at least , whereas when passing from the lower to upper level it is trivially 1, indicating no restriction on the transmitted momentum. The analogous restriction that a classical particle can only be transmitted to the upper level if it has sufficient kinetic energy is accounted for by the term.
We now discuss how, in appropriate limits, the formula essentially reduces to a LZ transition for each point in momentum space. We make a number of independent approximations:

.
For a single avoided crossing we may do this without loss of generality by shifting the space variable. 
.
This is the case, for example, when the potential is symmetric around the avoided crossing. 
is small.
This produces two simplifications to the formula using that :
The prefactor simplifies to ;

The factor in the exponential simplifies to .
Note that the small parameter in these expansions is actually , and so we expect these approximations to be more accurate for either small or large incoming momentum.


Second order expansion of .
It is well known Berry and Lim (1993) that a naïve second order expansion of is incorrect as the analytic continuation of must vanish like a square root at its complex zeros. We therefore approximate viawith a smooth function such that . Performing a second order expansion of then gives
In this case both and can be computed analytically to give
where . To connect purely to , we note that and hence
Finally, in order to connect to the LZ formulas, an explicit computation using that and gives
and so
Suppose now that we make all four approximations. Then the formula in (5) becomes
(6) 
It is now clear that the exponential factor corresponds precisely to the adiabatic LZ transition probability in (4), with the additional factor of accounting for the fact that we are determining the size of the transmitted wavepacket rather than the transition probability, which is proportional to the square of the wavepacket. The Heaviside function is also included indirectly in surface hopping models, which explicitly exclude classicallyforbidden transitions, see e.g. Belyaev, Lasser, and Trigila (2014).
Note that if we are interested solely in the transition probability then the first two approximations are irrelevant as they only affect the phase. However, when dealing with multiple transitions these terms are crucial in understanding interference effects. In Section VII we will investigate the effects of these approximations in some example systems.
After approximations (1)–(4) have been made, the resulting formula (6) can be thought of as a surface hopping algorithm that retains phase information. This can be seen by noting that the formula decouples in momentum space. Thus, if we replace the classical transport of individual particles, the ensemble of which represents the initial wavepacket, with quantum evolution of the initial wavepacket, and then replace particle hopping with hopping of momentum components of the wavepacket, then we have a clear analogue of the surface hopping methods. One promising avenue of further work is to investigate the use of formula 5 for the transmission probability (instead of the LZ one) in traditional surface hopping algorithms. Alternatively, we can recover the surface hopping methodology (but retaining phase information) by dividing the wavepacket into small pieces (the surface hopping particles), evolving them classically on the initial level (e.g. using Hagedorn’s wavepacket approach Hagedorn (1981, 1994); Lubich (2008)) until they reach an avoided crossing, and then applying the formula either with the full transition probability, or the LZ approximation, and reconstructing the wavepacket on the other level.
Vi The Algorithm
The general algorithm described below is similar to that presented in previous work, but here it is extended to multiple transitions and to different levels of approximation, which ultimately lead to an analogue of the LZ formula, but applied to wavepackets, rather than simply as a transition probability. The transmitted wavepacket is computed via the following algorithm. For clarity, we present the algorithm for two BOA surfaces, but its extension to multiple surfaces is straightforward due to the linearity of the Schrödinger equation.

Initial Condition: The initial wavepacket should be specified on either the upper or lower adiabatic level, well away from any of the avoided crossings. Note that, in such regions, the adiabatic, superadiabatic, and diabatic levels are very close, so one may instead specify the wavepacket on a single diabatic level. If the initial wavepacket is given close to an avoided crossing, for example as the result of a laser excitation, then it must be evolved away (into the scattering regime) on the corresponding adiabatic level under the BO approximation to obtain an appropriate initial wavepacket.

OneLevel Dynamics: The initial wavepacket is now evolved under the BOA until the final, specified time, or until another termination condition is satisfied (such as the wavepacket reaching a minimum distance from an avoided crossing). This can be done using any onelevel scheme that provides sufficient accuracy, such as Strang splitting, Hagedorn wavepackets Hagedorn (1981, 1994); Lubich (2008), or MCTDH HuarteLarrañaga and Manthe (2009). The wavepacket on the other BOA level is evolved simultaneously; the level is initially unoccupied.

Detection of Avoided Crossings: Here an avoided crossing is defined as a (local) minimum of the gap . Whenever the centre of mass of the wavepacket reaches such a minimum, apply the formula as described in the next step. Such local minima may be determined a priori, for example when the potentials are given analytically, or onthefly by monitoring .

Application of the Formula: Apply the formula (5) to the wavepacket at the avoided crossing and add the resulting wavepacket to the lower level. Note that the formula implicitly requires the potentials to be extended into the complex plane in order to compute . However, as described in the following Section, this requirement may be bypassed by using an analogue of the LandauZener formula, at a cost to accuracy which is investigated for some examples in Section VII.

Computation of Observables: At any time step the wavepackets on the two levels may be used to compute observables, such as mean position, momentum and the level populations, including those which require phase information such as interlevel observables. Note, however (as discussed in Section V.1), that these will only agree with the corresponding quantities computed for the adiabatic populations well away from any avoided crossings. An extreme example of this is that, before the wavepacket on the initial level reaches the avoided crossing, the other level us completely unoccupied; see Figure 8.
In the following section we will investigate the accuracy of this algorithm. One restriction for its application to multiple crossings is that the transmitted wavepacket must be small, or, more precisely, the wavepacket remaining on the original surface must not change significantly when compared to its evolution under the BOA. This is due to the perturbative nature of the derivation, which assumes that the original wavepacket is unchanged during a transition.
Vii Numerics
Note that the MATLAB code used to produce the results in this section is available from https://bitbucket.org/bdgoddard/qmd1dpublic/.
vii.1 JahnTeller
We consider first a simple example in order to demonstrate the effects of the approximations in Section V.3. We choose
where we have , , . There is a single avoided crossing at , with gap . It is clear that and a straightforward calculation shows that . Note, therefore, that assumptions (1), (2) and (4) of Section V.3 hold exactly. Furthermore, since , the diabatic and adiabatic LZ transition probabilities given in (3) and (4) are identical in this case. This simple model allows us to investigate the effects of approximation (3), i.e. the difference between the full formula (5) and the LZ approximation for a range of values of . From the arguments in Section V.3, we expect the two results to agree to high accuracy when is small, and hence the transition is large, but we expect the full formula result to be more accurate in the regime of interest (relatively large and small transitions).
We choose to specify the wavepacket at the avoided crossing, and determine the initial condition by evolving it backwards in time away from the crossing on a single adiabatic surface. This ensures that the wavepacket is semiclassical (i.e. of width order ) when it reaches the avoided crossing. As noted above, due to the linearity of the Schrödinger equation, if this were not the case then we could apply a slicing procedure to obtain similarly accurate results. In particular, we choose
(7) 
where, along with , the free parameters are and . For this example we fix , which is similar to the value chosen in surface hopping works e.g. Fermanian Kammerer and Lasser (2008); Belyaev, Lasser, and Trigila (2014); Fermanian Kammerer and Lasser (2017) (and approximately correct for realworld systems e.g. Betz, Goddard, and Manthe (2016)) and . We could, in principle, vary these parameters, and we will do so in later examples. We note that due to the nature of the potential, in order to start sufficiently far away from the avoided crossing (such that the adiabatic and superadiabatic representations agree) the initial potential energy must be reasonably large, leading to a correspondingly large minimum value of at the avoiding crossing. See Figure 3 for the adiabatic potentials.
Here we evolve backwards to a start time of with timestep and then forwards through the avoided crossing for time with the same timestep. We perform the numerics with a spatial grid with points and endpoints . We observe that halving the time step and doubling the number of grid points does not significantly affect the results.
In Figure 2 we show the relative error in the transmitted wavepacket and the transmitted mass. This clearly demonstrates that, for small (and large transmitted mass), both our formula (5) and the LZlike version (6) give very good results. However, as increases, the simplified version becomes increasingly inaccurate.
vii.2 Simple Avoided Crossing
We now consider a simple example, which will both allow us to systematically investigate the accuracy of our method for different parameter regimes, and provide a benchmark for the accuracy of a single transition; this is, at least heuristically, a lower bound for the accuracy for multiple transitions. We choose
where we have , , . See Figure 3 for the adiabatic potentials with .
As in the previous example, in order to control the (mean) momentum of the wavepacket when it reaches the crossing, we specify the wavepacket in momentum space at the avoided crossing and then evolve it backwards in time on a single adiabatic surface to obtain an initial wavepacket for the computations. In particular, we take a Gaussian wavepacket as given in (7) for a range of values of and . We compute the results for a single transition of the avoided crossing, both using the full formula (5) and the LZlike one (6). As can be seen from Figure 4, the relative error when using 5 is typically of the order of a few percent, with increasing accuracy as and/or increase. The deviation of the green curve, which corresponds to is a result of the asymptotic nature of the formula. The odd behaviour of the blue curve for , and seems to be a result of parts of the wavepacket becoming ‘trapped’ near the avoided crossing, which violates the assumption of a single transition.
Figures 5 and 6 demonstrate the effects of using the algorithm with the approximate formula (6). As can be seen from Figure 5, for moderate values of , the results become very poor. However, as expected, Figure 6 shows that, for small , the results are very similar to those using the full formula (5).
For completeness, we give the numerical details: The spatial grid uses points with limits . We use a time step of and obtain the initial wavepacket by evolving the wavepacket backwards from the crossing for time . The system is then evolved forwards for time . Again, we note that halving the time step and doubling the number of grid points does not significantly affect the results.
As a further test of the accuracy of the algorithm we perform the same calculation as for the Gaussian wavepackets in the previous example, but with a wavepacket on the upper level at the crossing given by
(8) 
where is a Gaussian as given by (7). We choose , , and . However, we note that the results are robust under these choices for a wide range of values. We show the resulting transmitted wavepacket in Figure 7 which for convenience of displaying the phase, we have evolved backwards to the avoided crossing on the lower level. Note that the relative error in this case is 0.0057. In particular, Figure 7 demonstrates that highermomentum wavepackets are more likely to make the transition.
We note here that the results for wavepackets starting on the lower level are very similar, and we will investigate such a situation in the following Section.
vii.2.1 Full Crossings
Here we consider if the algorithm is applicable to full crossings (with ). In such a case, the approximations made in Section V.3 lead to the conclusion that the transmitted wavepacket is approximately equal to the incoming wavepacket. Applying this in the case , and gives a relative of 0.0856 for both the formula (5) and LZ approximation (6). This is comparable to the relative error for small, but nonzero (see Figure 4). This indicates that the methodology can also be used for full crossings. This is important in higher dimensions, where part of the wavepacket may travel across a full crossing (conical intersection), whilst other parts experience an effective avoided crossing, in which case we need only one algorithm to accurately treat the whole wavepacket.
vii.2.2 Diabatic vs. Adiabatic LZ
In the previous examples, we have and hence the diabatic and adiabatic LZ transition probabilities, (3) and (4), respectively, are identical. However, here we briefly consider an example in which this is not the case. We note that such a situation was also investigated in Ref. Belyaev, Lasser, and Trigila (2014) for twodimensional crossings and the results when using the two formalisms were found to be very similar. We will now show that this is not always the case. We take a sightly perturbed version of the simple potential matrix above
i.e and . We choose , and where these parameters, are chosen such that we are in a regime where we expect both the formula (5) and the (adiabatic) LZ approximation to be reasonably accurate, whilst simultaneously the results are not dominated by being very small. We use the same numerical scheme as for the simple crossing above and, find that the relative errors when using the formula (5) and the adiabatic LZ approximation (4) (or (6)) are very similar at 0.0219 and 0.0217, respectively. In contrast, when using the diabatic approximation (3), the results are much worse, with a relative error of 0.1081. This, along with previous results, demonstrates a clear motivation to use the transition formula in surface hopping algorithms.
vii.3 Multiple Transitions of a Single Crossing
We now demonstrate the algorithm when the wavepacket makes multiple transitions of a single avoided crossing. Here we add a quadratic confining potential, which causes the wavepacket to oscillate backwards and forwards through the avoided crossing:
where , , and . We choose , which gives a relatively weak confining potential. We use the same grid and time step as for the simple case in Section VII.2 but here evolve back to and forwards to , which gives 3 complete transitions of the avoided crossing. Here we start with a wavepacket of the form (7) on the lower level with and . Again we choose . See Figure 3 for the adiabatic potentials.
As can be seen in Figure 8, the ‘exact’ dynamics require extreme numerical cancellations at each transition in order to produce the true wavepacket. Although not shown in the Figure, the maximum transmitted mass is 0.0028, which is around 200 times larger than the final mass. Note that the results of both formulas are of a similar accuracy to the results for a single crossing, with relative errors 0.0123 and 3.637 for (5) and (6), respectively. In particular, the agreement between the ‘exact’ and formula (5) results is excellent whilst, in this case, (6) significantly underestimates the size of the transmitted wavepackets.
Whist, in principle, we would expect the results of using (6) to improve when decreases (i.e. when the transmitted wavepacket is larger) this adds a complication to the algorithm: When the transmitted wavepacket is large, the transition significantly affects the wavepacket on the original level, which is used explicitly in the formula for the transmitted wavepacket at the next avoided crossing. Due to the perturbative nature of the derivation of the formula (5) (see e.g. Betz, Goddard, and Teufel (2009); Betz and Goddard (2009)), the wavepacket on the original level is not treated explicitly, and so we do not have access to this unless it can be assumed that it is largely unaffected by the transition. A necessary requirement for this, due to mass conservation, is that the transmitted wavepacket is small.
vii.4 Dual Avoided Crossings
As we have seen, for multiple transitions at a single avoided crossing, the algorithm described in Section VI works as expected, determining the correct phase of the wavepackets, and therefore also the correct interference effects. However, for transitions at separate avoided crossings there is an extra difficulty that arises from the definition of the diabatic and adiabatic potentials. As can be seen from Figure 9, in an example with two identical avoided crossings at , the diabatic eigenfunctions are even. Hence, treating the two crossings independently, the dynamics through the second crossing could be computed by flipping the surfaces in space (which gives the diabatic surfaces associated with the first crossing) and reversing the momentum of the wavepacket. From (5), we see that reversing the momentum introduces a sign change in the transmitted wavepacket, which must be taken into account when computing the total transmitted wavepacket. (This issue is related to the diabatic eigenfunctions only being defined up to their sign.) Note that this argument generalises to the case where there are multiple nonidentical crossings; the case here was chosen for clarity.
Here we choose potentials
where we have and . We choose , and an initial Gaussian condition using (7) with and . We use the same numerical scheme as in Section VII.2], first evolving backwards for and then forwards for . See Figure 3 for the adiabatic potentials.
As can be seen from Figure 10, the results using (5) are once again very good (with a relative error of 0.0295), whilst those using the approximate formula (6) are much poorer (relative error 2.193). Note that if we do not include the additional phase correction described above then the results using (5) are also very poor.
Viii Conclusions and Open Problems
We have presented a general scheme for the computation of wavepackets transmitted during multiple transitions through avoided crossings (at least when the transmitted wavepacket is small), which is also applicable to single transitions through full crossings. In fact, since, in the latter case, almost the entire wavepacket is transmitted, the scheme should also give accurate results for multiple transitions of full crossings, or combinations of a single full and multiple avoided crossings.
The principal advantage of our algorithm is that it produces the full quantum wavepacket, including its phase, in particular allowing the investigation interference effects during multiple transitions. This is in contrast to standard surfacehopping algorithms that lose all phase information, and cannot hope to treat systems with interference effects.
Open problems, which will be the subject of future works, are (i) Approximation of the wavepacket that remains on the original level when the transmitted wavepacket is not small, which would allow the study of multiple transitions of general crossings; (ii) Extension to higher dimensions. This can be done via a slicing algorithm; preliminary results for model systems are presented in Betz, Goddard, and Hurst (); (iii) Implementation of our more accurate transition rate in surface hopping models, which should extend their range of validity to systems when the transmitted wavepacket is significantly smaller than those which can be accurately captured by existing LZ schemes.
Acknowledgements.
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, HeriotWatt University and the University of Edinburgh.References
 Born and Oppenheimer (1927) M. Born and R. Oppenheimer, Ann. Phys. (Leipzig) 84, 457 (1927).
 Domcke, Yarkony, and Horst (2004) W. Domcke, D. Yarkony, and K. Horst, Conical intersections: electronic structure, dynamics and spectroscopy, Vol. 15 (World Scientific, 2004).
 Domcke, Yarkony, and Köppel (2011) W. Domcke, D. Yarkony, and A. Köppel, Conical intersections: theory, computation and experiment, Vol. 17 (World Scientific, 2011).
 Nakamura (2012) H. Nakamura, Nonadiabatic transition: concepts, basic theories and applications (World Scientific, 2012).
 Tully (2012) J. Tully, J. Chem. Phys. 137, 22A301 (2012).
 Hagedorn and Joye (2001) G. Hagedorn and A. Joye, Comm. Math. Phys. 223, 583 (2001).
 Hagedorn and Joye (2005) G. Hagedorn and A. Joye, Ann. Henri Poincaré 6, 937 (2005).
 Martinez and Sordoni (2002) A. Martinez and V. Sordoni, C. R. Math. Acad. Sci. Paris 334, 185 (2002).
 Von Neumann and Wigner (1929) J. Von Neumann and E. Wigner, Phys. Z. 30, 467 (1929).
 BenNun and Martínez (1998) M. BenNun and T. Martínez, J. Chem. Phys. 108, 7244 (1998).
 BenNun, Quenneville, and Martínez (2000) M. BenNun, J. Quenneville, and T. Martínez, J. Phys. Chem. A 104, 5161 (2000).
 Virshup, Chen, and Martínez (2012) A. Virshup, J. Chen, and T. Martínez, J. Chem. Phys. 137, 22A519 (2012).
 Heller (1991) E. Heller, J. Chem. Phys. 94, 2723 (1991).
 McLachlan (1964) A. McLachlan, Mol. Phys. 8, 39 (1964).
 Meyer and Miller (1979) H.D. Meyer and W. Miller, J. Chem. Phys. 70, 3214 (1979).
 Sawada, Nitzan, and Metiu (1985) S.I. Sawada, A. Nitzan, and H. Metiu, Phys. Rev. B 32, 851 (1985).
 Miller (1970) W. Miller, J. Chem. Phys. 53, 3578 (1970).
 Kreek and Marcus (1974) H. Kreek and R. Marcus, J. Chem. Phys. 61, 3308 (1974).
 Miller (2001) W. Miller, J. Phys. Chem. A 105, 2942 (2001).
 Tully and Preston (1971) J. Tully and R. Preston, J. Chem. Phys. 55, 562 (1971).
 Miller and George (1972) W. Miller and T. George, J. Chem. Phys. 56, 5637 (1972).
 Stine and Muckerman (1976) J. Stine and J. Muckerman, J. Chem. Phys. 65, 3975 (1976).
 Kuntz, Kendrick, and Whitton (1979) P. Kuntz, J. Kendrick, and W. Whitton, Chem. Phys. 38, 147 (1979).
 Tully (1990) J. Tully, J. Chem. Phys. 93, 1061 (1990).
 HammesSchiffer and Tully (1994) S. HammesSchiffer and J. Tully, J. Chem. Phys. 101, 4657 (1994).
 Müller and Stock (1997) U. Müller and G. Stock, J. Chem. Phys. 107, 6230 (1997).
 Fabiano, Groenhof, and Thiel (2008) E. Fabiano, G. Groenhof, and W. Thiel, Chem. Phys. 351, 111 (2008).
 Fermanian Kammerer and Lasser (2008) C. Fermanian Kammerer and C. Lasser, J. Chem. Phys. 128, 144102 (2008).
 Lasser and Swart (2008) C. Lasser and T. Swart, J. Chem. Phys. 129, 034302 (2008).
 Belyaev, Lasser, and Trigila (2014) A. Belyaev, C. Lasser, and G. Trigila, J. Chem. Phys. 140, 224108 (2014).
 Belyaev et al. (2015) A. Belyaev, W. Domcke, C. Lasser, and G. Trigila, J. Chem. Phys. 142, 104307 (2015).
 Zener (1932) D. Zener, Proc. Roy. Soc. London 137, 696 (1932).
 Landau (1965) L. Landau, Collected Papers of L.D. Landau (Pergamon Press, Oxford, 1965).
 Fermanian Kammerer and Lasser (2017) C. Fermanian Kammerer and C. Lasser, Commun. in Math. Phys. , 1 (2017).
 Chai et al. (2015) L. Chai, S. Jin, Q. Li, and O. Morandi, Multiscale Model. Sim. 13, 205 (2015).
 Hagedorn (1994) G. Hagedorn, Molecular propagation through electron energy level crossings, Vol. 536 (American Mathematical Soc., 1994).
 Hagedorn and Joye (1998) G. Hagedorn and A. Joye, Ann. Henri Poincaré 68, 85 (1998).
 Berry (1990) M. Berry, P. Roy. Soc. Lond. A Mat. 429, 61 (1990).
 Berry and Lim (1993) M. Berry and R. Lim, J. Phys. AMath. Gen. 26, 4737 (1993).
 Hagedorn and Joye (2004) G. Hagedorn and A. Joye, Comm. Math. Phys. 250, 393 (2004).
 Betz and Teufel (2005) V. Betz and S. Teufel, Ann. Henri Poincaré 6, 217 (2005).
 Betz, Goddard, and Teufel (2009) V. Betz, B. Goddard, and S. Teufel, Proc. Roy. Soc. A  Math. Phys. 465, 3553 (2009).
 Betz and Goddard (2009) V. Betz and B. Goddard, Phys. Rev. Lett. 103, 213001 (2009).
 Betz and Goddard (2011) V. Betz and B. Goddard, SIAM J. Sci. Comput. 33, 2247 (2011).
 Betz, Goddard, and Manthe (2016) V. Betz, B. Goddard, and U. Manthe, J. Chem. Phys. 144, 224109 (2016).
 Hagedorn (1981) G. Hagedorn, Ann. Phys. 135, 58 (1981).
 Lubich (2008) C. Lubich, From quantum to classical molecular dynamics: reduced models and numerical analysis (European Mathematical Society, 2008).
 HuarteLarrañaga and Manthe (2009) F. HuarteLarrañaga and U. Manthe, Multidimensional Quantum Dynamics: MCTDH Theory and Applications (Wiley, 2009).
 Hagedorn (1980) G. A. Hagedorn, Commun. Math. Phys. 77, 1 (1980).
 Spohn and Teufel (2001) H. Spohn and S. Teufel, Commun. Math. Phys. 224, 113 (2001).
 Belyaev et al. (2010) A. Belyaev, P. Barklem, A. Dickinson, and F. Gadéa, Phys. Rev. A 81, 032706 (2010).
 Zhu, Teranishi, and Nakamura (2001) C. Zhu, Y. Teranishi, and H. Nakamura, Adv. Chem. Phys. 117, 127 (2001).
 Belyaev and Lebedev (2011) A. Belyaev and O. Lebedev, Phys. Rev. A 84, 014701 (2011).
 Teufel (2003) S. Teufel, Adiabatic perturbation theory in quantum dynamics (Springer, 2003).
 Joye, Mileti, and Pfister (1991) A. Joye, G. Mileti, and C. Pfister, Phys. Rev. A 44, 4280 (1991).
 (56) V. Betz, B. D. Goddard, and T. Hurst, In preparation.