Nonadiabatic decay of metastable states on coupled linear potentials

Avoided crossings of level pairs with opposite slopes can form potential energy curves for the external degree of freedom of quantum particles. We investigate nonadiabatic decay of metastable states on such avoided crossings (MSACs) using diabatic and adiabatic representations. The system is described by a single scaled adiabaticity parameter, $V$. The time-independent two-component Schr\"odinger equation is solved in both representations, and the nonadiabatic lifetimes of MSACs are determined from a wave-function flux calculation and from the Breit-Wigner formula, leading to four lifetime values for each MSAC. We also solve the time-dependent Schr\"odinger equation in both pictures and derive the MSAC lifetimes from wave-function decay. The sets of six non-perturbative values for the MSAC lifetimes agree well, validating the approaches. As the adiabaticity parameter $V$ is increased by about a factor of ten, the MSAC character transitions from marginally to highly stable, with the lifetimes increasing by about ten orders of magnitude. The $\nu$-dependence of the lifetimes in several regimes is discussed. Time-dependent perturbation theory is found to yield approximate lifetimes that deviate by $\lesssim 30\%$ from the non-perturbative results, while predictions based on the semi-classical Landau-Zener tunneling equation are found to be up to a factor of twenty off, over the ranges of $V$ and $\nu$ studied. The results are relevant to numerous atomic and molecular systems with quantum states on intersecting, coupled potential energy curves.

In this paper, we present a non-perturbative, quantummechanical analysis of nonadiabatic decay of low-lying metastable states at avoided crossings (MSACs). Similar descriptions have previously been employed to model wave-packet dynamics on intersections [33,36] and in Rydberg-ground molecules [37,38]. Here, we concentrate on the nonadiabatic lifetimes of quasi-stationary MSACs, which are important in the aforementioned applications. After explaining our model and the utilized techniques in Sec. II, in Sec. III we obtain solutions of the time-dependent and time-independent Schrödinger equations in both diabatic and adiabatic representations. We extract nonadiabatic MSAC lifetimes from six nonperturbative methods, and compare and interpret the results. The analysis is performed for a range of coupling strengths between the diabatic potentials, for MSACs with vibrational quantum numbers ranging up to about 15. In Sec. IV, we compare the non-perturbative MSAC lifetime results with estimates based on time-dependent perturbation theory, and with semi-classical estimates based on the LZ formula. The paper is concluded in Sec. V.

A. System under study
In the system of interest, the physical Hamiltonian in the diabatic representation, acts on a two-component wave function (ψ 1 (x), ψ 2 (x)) with a position-independent, internal state space denoted as {|1 , |2 }, in that order. The constants α p and V p are chosen positive and real. The effective particle mass is denoted M , and the external degree of freedom has a spatial coordinate x p . The diabatic energies of the internal states {|1 , |2 } as a function of x p are V 1 = −α p x p /2 and V 2 = α p x p /2, respectively, with differential slope α p , and the constant coupling between these states is V p . For convenient description of different physical systems, we use the following units for length, energy, time and frequency, length : l 0 = Expressing length and energy in these units, the Hamiltonian in Eq. 1 transforms into the scaled Hamiltonian in diabatic representation, with scaled position x = x p /l 0 and scaled coupling strength The characteristic half width of the crossing region in scaled units is x w = 2V ; in physical length units it is x wp = V l 0 = 2V p /α p . The scaled coupling V serves as an adiabaticity parameter: the larger V , the more adiabatic a system will behave, and the less affected the MSACs will be by nonadiabatic decay. In the following, we will use the scaled units defined in Eq. 2. The x-dependent adiabatic-state basis {|u , |d } for the internal degree of freedom, and the adiabatic potentials V u and V d , are defined byĤ D |u = V u (x)|u andĤ D |d = V d (x)|d , with u and d standing for "up" and "down" in energy, V u positive, and V d = −V u . With the notation |u (x) = i=1,2 χ u,i (x)|i and |d (x) = i=1,2 χ d,i (x)|i , the first-and second-order nonadiabatic couplings are There, the index i denotes diabatic and the Greek letters adiabatic basis states. The χ α,i can be chosen real. The 2x2 matrix A α,β then is anti-symmetric at any value of x, with α and β being u or d. The diagonal elements of B α,β (x) are compounded with the adiabatic potentials to yield the potential energy curves (PECs) Hamiltonian, which is a special case of the Born-Huang representation [39,40] for the case studied in our paper, then writeŝ This Hamiltonian acts on the adiabatic wave functions, (ψ u (x), ψ d (x)). As visualized in Fig. 1, the nonadiabatic A-and B-couplings vanish for x ≫ V , with V from Eqs. 3 and 4.

Diabatic representation
A straightforward method to arrive at a nonperturbative solution is to solve the time-independent Schrödinger equation (TIDSE) for the Hamiltonian from Eq. 3. Here, we are interested in the energy range W > V , where MSAC resonances exist. The MSAC resonance energies and corresponding two-component wave functions are obtained numerically. As a spatial integration method, we have chosen a 4-th order Runge-Kutta (RK) method, which allows for first-derivative terms (needed in the adiabatic representation discussed in Sec. II B 2). In the following, relevant details are explained.
It is well-known from textbooks that for a spin-less particle on a linear potential the wave-function solutions are given by Airy functions (see, e.g., [41]). In the case of two coupled linear potentials, as in Eq. 3, the matching of the boundary conditions in the classically forbidden regions turns out to be numerically delicate due to the coupling V between the classicallyallowed, Airy-function-like solutions to the co-located classically forbidden ones. In the asymptotic regions, the allowed solutions are, locally, approximately given by a(x) cos(kx + φ), with a slowly-varying local amplitude a(x), wave number k(x) = 2(|x|/2 + W ), the quantum state's scaled energy W , and a phase φ. For large |x|, the classically-forbidden solutions are then approximately given by − a(x)V |x|/2+W cos(kx + φ). The amplitudes of the forbidden solutions drop off quite slowly in |x|, because V is fixed and never "turns off". In the numerical implementation, this exacerbates the tendency of the classically-forbidden solutions to exponentially diverge at large |x|. The issue is addressed by choosing sufficiently small values for the spatial step size, ∆x, and for the slope iteration parameter, s, explained in the next paragraph. The issue is less pronounced in the adiabatic approach, because the nonadiabatic A-and B-couplings both do "turn off" at large |x| (see Sec. II B 2). The energy spectrum of the Hamiltonian in Eq. 3 is continuous and ranges from −∞ to ∞. The numerical treatment is simplified by the symmetry of the realvalued solutions. For each energy W there exists an even and an odd solution. Even solutions, which are associated with even-parity MSACs, are of the form ψ 1 (x) = 1 − sx and ψ 2 (x) = 1 + sx for |x| → 0, with a slope parameter s. The odd solutions are of the form ψ 1 (x) = 1 + sx and ψ 2 (x) = −1 + sx for |x| → 0. Further, for any x it is ψ 2 (−x) = ψ 1 (x) for the even and ψ 2 (−x) = −ψ 1 (x) for the odd solutions. For any energy W , this leaves only one parameter -the slope s -to be iterated. In both even and odd cases, the slope parameter s is iterated to minimize the classically forbidden wavefunction components at the chosen spatial-range limit, |x| = x max . We vary x max depending on V and W , so as to allow for maximum outward propagation before the wave functions diverge due to numerical inaccuracies. For each energy W , this procedure yields exactly one even and one odd solution.

Adiabatic representation
TIDSE in the adiabatic picture has the Hamiltonian from Eq. 6. The even adiabatic solutions are of the form ψ u (x) = 1 and ψ d (x) = sx for |x| → 0, and the odd ones are of the form ψ u (x) = sx and ψ d (x) = 1 for |x| → 0. As in Sec. II B 1, the slope parameter s is iterated to minimize the classically forbidden wave-function components ψ u (x) at the spatial-range limits, |x| = x max . For each energy value W , there exist exactly one even and one odd solution. In the numerical treatment, the tendency of the classically forbidden solutions on the respective PECs to exponentially diverge at large |x| is less pronounced in the adiabatic representation than it is in the diabatic representation (Sec. II B 1), because the nonadiabatic Aand B-couplings "turn off" at large |x|. In contrast, in the diabatic representation the constant coupling V does not "turn off" at large |x|.

MSAC resonances
The energies W ν of the MSAC resonances, labeled by an integer vibrational quantum number ν, can be determined iteratively by locating the energy values at which the amplitudes of the sinusoidal wave-function tails in the respective classically-forbidden regions become minimal near the edges of the spatial integration range, |x| = x max . We label the resonances starting with ν = 0 for the MSAC ground state. The coupling parameter V is varied between 0.3 (least adiabatic) and 2.8 (most adiabatic), in scaled units as defined in Eq. 2. We find all MSAC resonances within an energy range of about V < W V + 3.8. For V ranging between 0.3 and 2.8, the number of MSACs with V < W ν V + 3.8 ranges from ν + 1 = 9 to 15, respectively. The integration limit, x max , is shifted outward with increasing V and ν in order to locate the MSAC energies as accurately as possible over the entire V -and ν-range studied. Here, x max is varied between x max = 13 at the lowest V and ν, and x max = 19 at the largest V and ν.
For illustration, in Fig. 1 we show plots of the adiabatic potentialsṼ d andṼ u and the A-and B-potentials for V = 0.306 and V = 1.5275, as well as the obtained lowest MSACs. In addition to Fig. 1 illustrates the rapid drop in amplitude of both the Aand B-potentials with increasing V . The effect of the diagonal B-potentials, B dd (x) and B uu (x), only becomes apparent in the V = 0.306-case in the form of small humps in the range |x| 1. The A-couplings generally appear to be more important than the B-couplings, as confirmed directly in Sec. IV A. A feature that becomes important in the interpretation of the ν-dependence of the MSAC lifetimes in Sec. III B is that at low V the approximate reach of the A-and B-potentials, given by the crossing half width, x w = 2V , is smaller than the typical wave function extents, whereas at large V the Aand B-potentials are spread out over the entire typical wave-function extent.

MSAC lifetimes from flux calculation
The main interest in the present work is in the nonadiabatic lifetime of the lowest MSAC resonances. To that end, we compute even and odd solutions on a dense grid of the continuous energy W , and determine the resonance centers, W ν , as described in Sec. II B 3. In either representation, the MSAC resonances correspond with wave-function solutions that minimize the amplitudes of the sinusoidal wave-function tails in the respective classically forbidden regions at the edges of the spatial integration range, ±x max . As seen in Fig. 2, in the asymptotic regions both classically allowed and forbidden wave-function tails are locally of the form ψ(x) = a(x) cos[k(x)x + φ(W )], with a slowly-varying amplitude a, an energy-dependent phase φ, and a slowly-varying wave number k. In the diabatic representation, the oscillatory behavior of the classically-forbidden tails results from the fixed coupling V , which induces π-out-of-phase classically-forbidden tails. Denoting the amplitude of the classically allowed tail a a , and that of the co-located classically forbidden tail a f , at the spatial integration boundary x max the amplitude a f ∼ aaV xmax/2+W can be on the order of 10% of a a [see, for instance, the inset of Fig. 2 (a)]. In the adiabatic representation, the oscillatory behavior of the classically-forbidden tails primarily results from the diminishing nonadiabatic A-coupling, which induces π/2-out-of-phase forbidden tails with much smaller amplitudes than in the diabatic representation [note that in the inset in Fig. 2 (b) the classically forbidden tail of ψ u is magnified by a factor of 100].
The steady-state solutions can be viewed as superpositions of in-going "pump" waves with out-going backscattered waves, forming perfect standing waves on either side of the potential. For a scalar wave function with oscillatory tails of amplitude a in the positiveand negative-x domains, ψ( , summed over the positive-and negative-x domains, is 2k|a| 2 .
If the wave function contains a metastable resonance in a central potential well, the decay rate of the resonance upon turning off the pump waves is Γ = 2k|a| 2 /P 0 , where P 0 is the wave-function norm. Since P 0 is proportional to |a| 2 , the factor |a| 2 drops out.
We first discuss the implementation in diabatic representation, in which the norm P 0 of the MSAC wave . In this way, the limit x k is set such that P 0 captures the decaying tails of the MSAC resonances in the classically-forbidden regions to within r k semi-classical 1/e decay lengths outside the classical turning points. Here we use r k = 3, which is large enough for P 0 to capture the entire resonance norm, and small enough to not include substantial probability from the oscillatory wave-function tails. The exact value of r k is not important. Note that we define the boundaries via the upper adiabatic potential, V u (x), in both the diabatic and adiabatic representations. The limits x l and x k are visualized in Fig. 2.
The amplitudes a 1 and a 2 of the oscillatory wavefunction tails are found by first locating a position x p where the (classically allowed) tail of ψ 1 (x) takes an extremal value close to the positive limit of the spatial integration range, x max [see circle in the inset of Fig. 2 Using three adjacent carrier points of each ψ i (x) with x p at the center, we compute the wave numbers Due to the positionindependence of the diabatic internal states, {|1 , |2 }, the fluxes in the two wave-function components just add up. The 1/e lifetime τ nad,F C = 1/Γ nad,F C , obtained from the time-independent MSAC wave function in the diabatic picture, then is given by It is noted that k 1 ≈ k 2 in all cases, whereas the ratio a 2 /a 1 increases with V .
We also obtain the lifetimes in the adiabatic representation, in which the wave-function computation is numerically more stable. The norm integral P 0 is computed with the same boundary x k as in the diabatic representation, We first find a peak location x p of the (classically-allowed) tail of ψ d (x) near the integration limit, x max [see circle in the inset of Fig. 2 (b)]. In the flux calculation in the adiabatic representation, the x-dependence of the adiabatic internal-state basis {|u (x), |d (x)} must be considered. Therefore, the adiabatic two-component wave function (ψ d (x), ψ u (x)) is transformed into diabatic representation, (ψ 1 (x), ψ 2 (x)) at three adjacent x-values centered at x p . The decay rate Γ ad,F C and the 1/e lifetime τ ad,F C of the time-independent MSAC wave function in the adiabatic picture are then computed from (ψ 1 (x), ψ 2 (x)) at x p , using equations from the previous paragraph.

MSAC lifetimes from the Breit-Wigner formula
In an alternative, quite different method, we also obtain the MSAC lifetimes from the Breit-Wigner formula (BW) [42]. In the asymptotic regions, time-independent real-valued solutions on the classically allowed potentials are locally of the form ψ(x) = a cos(k(x)x + φ(W )), with an energy-dependent phase φ. The asymptotic solution is a superposition of incident and back-scattered waves of respective forms exp(−ikx) and exp(i(kx + 2δ), with the usual scattering phase shift δ [42]. It is thus seen that the phase φ in the time-independent solution equals the scattering phase, δ = φ. According to the BW formula, the decay rate of a MSAC, at the center of the scattering resonance, is given by Γ BW = 2(dW/dφ), where the derivative is taken at a fixed location x B well outside the classically allowed range of the bound component of the MSAC wave function. Here we pick a location close to x max ; the exact value of x B is not important. The phase is then obtained from the classically-allowed tails of the wave functions ψ 1 (x) or ψ d (x) at x B , in the diabatic and adiabatic representations, respectively, using where the integer m is continually adjusted as a function of W for continuity of φ(W ). We subtract a background phase φ 0 (W ) that arises from the phase shift of the non-resonant solutions away from the MSACs and that is computed from where the potential V * (x) = −x/2 in the diabatic and V * (x) = V d (x) in the adiabatic representation. Note that for vanishing coupling, V = 0, the phase would be that of an Airy-function solution [41]). The BW decay rates and lifetimes then become where * = nad and * = ad for the diabatic and adiabatic representation, respectively. In summary of this subsection, we obtain four values for the nonadiabatic decay times of MSACs from solutions of time-independent two-component Schrödinger equations in diabatic and adiabatic representation, namely τ nad,F C , τ nad,BW , τ ad,F C , and τ ad,BW . As expected and shown below, these generally agree very well with each other, with the values from the adiabatic picture being more accurate due to the vanishing of the Aand B-coupling terms at large |x|.
Here, we also have B ud (x) = −B du (x) for all x. The TDSE in diabatic representation follows from Eq. 3.
As initial conditions for the MSAC wave functions at time t = 0 in the diabatic and the adiabatic representations, we use the respective time-independent solutions obtained in Sec. II B 3. The MSAC wave functions from Sec. II B 3 exhibit oscillatory tails near the boundaries of the integration grid, as seen in Fig. 2. To avoid numerical instability, the MSAC wave functions entered as initial states are set to zero between their outermost nodes and the respective spatial integration boundaries, ±x max .
At the core of the TDSE method is to absorb the outgoing flux and to eliminate reflections from the boundaries [43]. The wave-function norms then drop exponentially, thereby revealing the decay time of the MSAC entered as initial state. The absorption is implemented by padding all diagonal potentials with imaginary absorbing layers near the spatial integration boundaries at ±x max . The absorbing layers rise smoothly from zero at locations well-outside the classical turning points, ±x l , to a maximal value at ±x max . The utilized time-propagation method is a Crank-Nicolson scheme [44] that is similar to schemes used in our recent work on tractor atom interferometry [45] and Rydberg-ion molecules [26], where nonadiabatic transitions were quantitatively described. More details on the method can be found there. In the present work, the time-dependent computations are performed with a spatial-grid step size of ∆x = 10 −3 , the same as in the time-independent methods described in Sec. II B, and a time-step size of ∆t = 10 −3 (all in scaled units). We have checked that a reduction of ∆t does not significantly affect the lifetimes found for the MSAC wave functions. The TDSE computations in the diabatic and adiabatic representations yield MSAC lifetimes denoted τ nad,T DSE and τ ad,T DSE , respectively.

A. Comparison of methods
In Fig. 3 (a), we first present a comparison of results for a moderately adiabatic case, V = 1.528, for ν = 0 to 11. The log-scale plot shows excellent agreement of lifetime data from all six methods over the entire range of ν, over which the lifetime drops by about a factor of 30. Among the methods, we consider the adiabatic wavefunction flux results, τ ad,F C , to be the most accurate and precise for the following reasons. The adiabatic analysis is less prone to numerical inaccuracy in the classicallyforbidden tails of the wave functions, because the adiabatic couplings A and B drop off rapidly with increasing |x| (see Figs. 1 and 2, and arguments presented in Sec. II B). This reduces the amplitude of the classicallyforbidden tails, thereby alleviating their trend towards exponential divergence. Further, the flux method is insensitive to background-phase effects, which affects the BW method at low V (see Sec. III B).
To exhibit small deviations of the results of the other five methods from τ ad,F C , in Fig. 3 (b) we show the ratios τ * /τ ad,F C , with * denoting the other methods. Im-portantly, the values for τ deviate by no more than 11% from τ ad,F C . The four results from the TIDSE agree to within 2% from each other, for V = 1.528, with small deviations attributed to numerical inaccuracy and to the systematic inaccuracy of the BW method at low V (see Sec. III B). The computations based on the TDSE deviate by up to 11% from τ ad,F C . This may be due to the susceptibility of the time-dependent computations to imperfections of the absorbing-wall implementation, such as less-than-perfect absorption of the outgoing flux and spurious reflections. Indeed, for ν ≤ 4, where the absorbing walls are the farthest away from the high-amplitude regions of the MSAC wave functions, the TDSE lifetime results deviate by less than about 5% from the TIDSE results. It is also noted that the diabatic and adiabatic TDSE calculations differ by less than about 1% from each other for all ν-values. This indicates that numerical issues, such as spatial-step or time-step sizes, introduce about the same, %-level of uncertainty in the TDSE and the TIDSE calculations.
Overall, the close agreement across the six methods in Fig. 3 proves the fundamental validity of all methods used. The quite good agreement between the TDSE and TIDSE calculations provides a particularly high level of validation, as the methods of how to extract the lifetimes from the TDISE and TDSE computations are quite different, yet both approaches yield very similar results.

B. Lifetimes vs adiabaticity
A main outcome of the work are the MSAC lifetimes over a wide range of the adiabaticity V and the vibrational quantum number ν. We have performed computations for a set of V -values ranging from V = 0.306 (least adiabatic) to V = 2.75 (most adiabatic). To assist with the interpretation of various regimes, we define the quality factor, Q, of the resonances as the angular frequency in harmonic approximation times the state's norm divided by the time derivative of the norm, or Q = τ ω = τ /(2 √ V ). There, ω is evaluated from the adiabatic potential V u (x) near x = 0. Note Q is unit-less and the same in scaled and physical units.
In Fig. 4 (a) we show the Q-values for the MSAC vibrational ground state, ν = 0, versus V . Noting that the number of oscillations after which the survival probability drops below 50% is approximately Q/9, it is seen that the ground-state MSACs may be considered only barely oscillatory for 0.3 < V 0.6, as in these cases it only takes a few oscillation periods or less for half of the groundstate MSAC population to decay. For V ≈ 1 it already takes a few tens of oscillation periods before the groundstate MSACs are half decayed. However, as V rises above about a value of 2, the ground-state MSACs quickly become highly stable against nonadiabatic decay. At the largest V -value tested, V = 2.75, it takes > 10 9 oscillation periods for half of the ground-state MSAC population to nonadiabatically decay [see Fig. 4 (a)]. The rapid The wide range of MSAC level damping is further visualized in Fig. 4 (b), where we show four examples of the wave-function scattering phases, φ(W )−φ 0 (W ), that are used for the calculation of BW lifetimes according to Sec. II B 5. At the MSAC energies, W ν , the phases exhibit rises in steps of π. The energy widths of the rises drop from a large fraction of the level spacing at V = 0.306 to too narrow to be visible at V = 2.444. Fig. 4 (b) reiterates the vast range of nonadiabatic damping behavior that is seen over the range 0.3 < V < 2.444, a range over which V varies by about one and Q by about ten orders of magnitude. Fig. 4 (b) also shows that at the lowest V -values the resonances are wide enough and the slopes at the resonances, d(φ−φ 0 )/dW , are small enough that background trends and cross talk between neighboring MSACs will affect the d(φ − φ 0 )/dW -readings at the resonance centers, W ν . This makes lifetimes from the BW formula inaccurate at low V , as seen below. Lifetimes from flux calculations are not susceptible to this type of inaccuracy.
In Fig. 5 we show lifetime results from TIDSE computations for ten values of V for MSACs within an energy range of about 3.8 s.u. from the potential minima of V u . (The computationally more intensive TDSE computations were performed only for the intermediate case of Fig. 5 (a) demonstrates good agreement between the TIDSE methods over a wide range of conditions. For all ν-values studied, the MSAC lifetimes increase by six to ten orders of magnitude, as V is increased from 0.306 to 2.75. In the following, we discuss the dependence of the lifetimes on ν in several regimes of V .
In the nonadiabatic regime, V 0.6, the lifetime barely depends on the vibrational quantum number, ν, and for the least-adiabatic case, V = 0.306, the lifetime actually increases with ν. This behavior, which may seem counter-intuitive at first, reflects the fact that for V ≪ 1 the anti-crossing half width, x w = 2V , which is an estimate for the reach of the A-and B-couplings, only is a fraction of the spatial extent of the MSAC wave function onṼ u , as seen above in Fig. 1 (a). As a result, for V ≪ 1 the spatial extent of the interaction range that causes the nonadiabatic decay, measured relative to the wave-function extent, decreases with increasing ν, leading to an increase in lifetime with increasing ν. This mechanism becomes more transparent in an analysis based on Fermi's Golden Rule (see Sec. IV A). Arguing semi-classically, one may say that at the lowest V -values studied the lifetime increases with ν because with increasing ν the MSAC particle spends less of its time in the anti-crossing region. It is noted that increasing ν for the purpose of increasing the MSAC lifetime is not a useful concept to generate long-lived MSACs (for atom trapping, for instance), because of the generally very low Q-values at V 0.6 [see Fig. 4 (a)]. For V 1.2, the MSAC resonances become increasingly adiabatic, with Q-values beginning to range above 100. In the adiabatic regime, the MSAC lifetimes decrease with increasing ν, which is opposite to the trend that is seen in the nonadiabatic regime. The decrease of τ with increasing ν accelerates with increasing V ; at V = 2.75, the largest value studied, the lifetime ratio between ν = 0 and ν = 10 exceeds a factor of 1000. In order to understand this behavior, one may first compare the relative importance of the A-and B-coupling terms in the adiabatic representation. It is found in Sec. IV A that the A-term is quite dominant. As a consequence, at sufficiently large V , the gradient of the trapped wave function, ∂ ∂x ψ u , averaged over the wave-function extent, factors decisively into the nonadiabatic coupling strength. This means that, at the larger V -values, the lifetime should drop with increasing ν, as observed. Noting that wave-function gradient and classical velocity are related, the velocity dependence of the Landau-Zener equation predicts the same trend (see Sec. IV B).
Next, we discuss the deviations between the lifetimes obtained with the TDISE methods. For visibility of small deviations, we display the ratios τ * /τ ad,F C on a fine scale in Fig. 5 (b). The adiabatic and nonadiabatic flux-calculation results agree very well in all regimes. We reiterate that τ ad,F C is still considered to be the most accurate and precise (see Sec. III A). For τ ad,F C 200, which roughly corresponds with V 1, the BW data also agree well. However, for τ ad,F C 200 they yield up to about 20% shorter lifetimes than the flux methods. It is also noted that the two BW results from the diabatic and adiabatic representations still agree very well with each other. The systematic deviation of the BW from the flux-calculation data at low V (nonadiabatic regime) may be attributed to the facts that at low V neighboring BW resonances begin to cross-talk, and that background phase slopes become a significant fraction of the slopes d(φ − φ 0 )/dW at the resonance centers [see Fig. 4 (b)], rendering the BW data less accurate at low V . It is further seen that the numerical noise of the diabatic BW calculations can reach 5% at large V , where the resonances become extremely narrow and the computation of the slopes d(φ − φ 0 )/dW becomes less accurate. Notwithstanding, the overall good agreement, seen on the fine scale in Fig. 5 (b), validates methods and results across the entire V -and ν-regimes studied.

A. Perturbation theory
The adiabatic representation lends itself to a perturbative description of nonadiabatic decay [33,46]. In this approach, we find bound MSACs onṼ u (x) used in Eq. 6, neglecting the nonadiabatic A-and B-couplings (but keeping the diagonal B-terms). These states differ from the true MSACs in that they are infinitely-long-lived, and in that their energies, W ν,F GR , are up to 0.07 scaled units below the true resonance energies, W ν . The energy deviations are most notable at small coupling V , where the off-diagonal non-adiabatic terms are large and cause the largest shifts W ν − W ν,F GR . We denote the wave functions of the coupling-and decay-free approximations of the MSACs as ψ u,ν,F GR (x). The ψ u,ν,F GR (x) are weakly coupled to the continuum of free-particle states on the potentialṼ d (x). The solutions onṼ d (x) are, asymptotically, identical with Airy functions [41]. Factoring in that onṼ d (x) the wave functions extend to both ±∞, as opposed to just one side on a linear potential, we normalize the free states such that the amplitude of their oscillatory tails at large positive x is There, W is the level energy. The solutions ψ d,W,F GR (x) normalized in that way are orthonormal in unit energy, i.e. it is ψ d,W ′ ,F GR |ψ d,W,F GR = δ(W − W ′ ). According to Fermi's Golden Rule (FGR), the transition rate from |ψ u,ν,F GR to |ψ d,W,F GR then is given by Γ F GR = 2π|M | 2 , with a matrix element where we abbreviate ψ u (x) = x|ψ u,ν,F GR and ψ d (x) = x|ψ d,W,F GR . The free-particle energy in the integral equals that of the quasi-bound state, W = W ν,F GR . Also, here all ψ(x) are real, and the integration range is limited by the range of ψ u (x). Since in the present problem B du (x) and d/dx have odd and A du (x) has even parity in x, even ψ u (x) decay into odd solutions ψ d (x) and vice versa. The FGR lifetimes then are The FGR calculation is visualized in Fig. 6 for a small and a large V -value, for the case ν = 3. While the bound and free wave functions, ψ u (x) = x|ψ u,ν=3,F GR and ψ d (x) = x|ψ d,W,F GR , look quite similar in the two cases [see Fig. 6 (a,d)], the coupling matrix elements are very different.
We define the matrix- in Fig. 6 (b,e). For V = 0.306 the m-densities are large, localized to within just the central lobe of ψ u (x), and largely uni-polar, whereas for V = 2.75 the m-densities are weak, spread-out over the entire reach of ψ u (x), highly oscillatory and bi-polar. In both cases, m A (x) is much larger than m B (x) in magnitude, on average. As a result, for small V the nonadiabatic decay is fast and largely driven by couplings localized to within a small, interior region of ψ u (x), whereas for large V the nonadiabatic decay is slow and spread-out over the entire range of ψ u (x). These observations validate statements that we have made in Sec. III B with regard to the ν-dependence of the MSAC lifetimes.
In Fig. 6 (c,f) we show the integrals of m(x), whose asymptotic values, M = M (x max ), give the FGR lifetimes according to Eq. 14. Due to symmetry, the integral M (x) has odd parity about a symmetry point at x = 0 [crosshair in Fig. 6 (c)], and it is M = M (x max ) = 2M (x = 0). For low and moderate values of V , the large amplitudes and the somewhat uni-polar characteristics of m(x) lead to numerically stable results for M and τ F GR . At large V , however, the integral in Eq. 13 is numerically challenging because of the bipolar and highly oscillatory behavior of m Σ (x). It is seen in Fig. 6 (f) that at large V the integral M = m Σ (x)dx comes down to a very small, nearly-vanishing remainder after integration, as evidenced by the fact that M (x) has a near-perfect zero crossing at x = 0, leading to a very small matrix element M . To get converging values for M , at the largest Vvalues studied we had to decrease the spatial step size in the wave-function computations and in the integral for M by a factor of up to about 100 relative to the step size used in the non-perturbative methods. Nevertheless, even at the largest V considered the FGR computations are still quite fast because the wave functions to be computed are scalar.
In Fig. 7 we present the ratio τ F GR /τ ad,F C for all values of V and ν also shown in Fig. 5. As in Figs. 3 (b) and 5 (b), τ ad,F C is used as a reference because the nonperturbative adiabatic wave-function flux calculation is the most accurate and precise. The lifetime ratios are plotted on a log scale covering about two decades, which is fine enough to observe relative deviations as small as about 1% and wide enough to also cover relative deviations for a Landau-Zener model (see Sec. IV B). The τ F GR /τ ad,F C -ratios, plotted in Fig. 7 versus τ ad,F C , follow a quite well-defined trend line at 0.1 to 0.3 below unity, with the lowest deviations occurring in the nonadiabatic and adiabatic limits on the left and right margins of the plot, respectively. At the largest τ ad,F C -values, corresponding to large V -and low ν-values, there is additional numerical noise on the order of ±0.1, caused by the delicate nature of the M -matrix elements at large V (see Fig. 6 and related discussion).
The FGR approach in this work differs from typical applications of FGR in which the wave functions are perturbation-independent and the perturbation has a tunable strength. In contrast, in the present case the perturbation V is fixed for a given set of wave functions ψ u (x) and ψ d (x), and the wave functions themselves depend on the fixed perturbation V . The matrix-element "densities" m(x) have a complex spatial structure and are considered in first order only. The deviations of the FGR from the non-perturbative results are notable, albeit not exceeding about 30%. A practical concern relies in the fact that at large V the spatial step size in the FGR calculation of the matrix element M has to be set very small to achieve convergence, due to the delicate nature of the M -integral at large V (see Fig. 6).

B. Landau-Zener model
For a semi-classical estimate of MSAC lifetimes using the Landau-Zener equation, we use a Landau-Zener tunneling "attempt rate" of twice the vibrational frequency, which gives an attempt rate of R(ν) = (W ν+1 − W ν−1 )/π (scaled units). The LZ coupling equals V and the differential slope of the diabatic potentials equals s = 1, in scaled units. For a fixed particle velocity, v, the Landau-Zener tunneling probability is P LZ = exp(−2πV 2 /(s v)), and the lifetime τ LZ = 1/(RP LZ ), in scaled units. Assuming that a semi-classical picture with a point-particle velocity v suffices to describe the quantum problem of interest, one needs a rule for how to get v. From Fourier transforms of MSAC wave functions in any representation (diabatic or adiabatic), one expects and finds that v could be on the order of √ W ν − V , which also accords with the classical virial theorem for a harmonic oscillator. Further, classically the velocity peaks at v = 2(W ν − V ) at the crossing. For the largest V and lowest ν studied in this work, these v-values produce τ LZvalues that are about 20 orders of magnitude too long. As the exponent in the Landau-Zener tunneling probability is ∝ −1/v, we may surmise that the high-velocity wings in the Fourier transforms of the MSAC wave functions govern the LZ decay rate. Empirically, one finds that v = √ 2W ν , used in the following, overall leads to the best LZ estimates for the MSAC lifetimes (that can still be several orders of magnitude off).
The deviations of τ LZ from quantum calculations are shown in Fig. 7 in terms of τ LZ /τ ad,F C . It is seen that, over our range in V and ν studied, the LZ model may serve as a very rough guideline to predict MSAC lifetimes, as the τ LZ -values stay within a factor of about 20 from τ ad,F C . The inaccuracy of the τ LZ -values accelerates in the adiabatic region (large τ ad,F C ). The strong ν-dependence of τ LZ /τ ad,F C , seen especially in the adiabatic region, reiterates that we have no well-founded rule for the classical velocity v. As such, the poor overall agreement of τ LZ with the quantum results reflects the fact that a semi-classical model applied on a problem in the quantum domain cannot be expected to be accurate.
Considering quantum-classical correspondence, we add that with increasing ν our model system becomes more classical, and with decreasing V the nonadiabatic transitions become relatively well-localized in the spatial region near x = 0. As a result, for V 1, and for V 1 and ν exceeding a V -dependent limit evident from Fig. 7, the τ LZ -values deviate by less than about 50% from the corresponding τ ad,F C -values, and the agreement improves with increasing ν. These observations accord with the expectation that quantum-classical correspondence should occur in the limit of large quantum numbers.

V. CONCLUSION
We have computed nonadiabatic lifetimes of metastable states on symmetric avoided crossings. Among six non-perturbative quantum methods, the results of which generally agree well, a wave-function flux method implemented in the adiabatic representation is the most accurate and precise, with lifetime uncertainties estimated at about 1%. Using the given relations between scaled and physical units, the results are portable to a variety of applications, including Rydberg molecules [20,22] and atom trapping and guiding on dressed potentials [32].
In addition to providing accurate, non-perturbative lifetime data, our comparisons have shown that timedependent perturbation theory in first order, applied to states in the adiabatic representation, with the nonadiabatic coupling terms treated as a perturbation, yields ap-proximate lifetimes that deviate by less than about 30% from the non-perturbative values. Semi-classical estimates based on the Landau-Zener tunneling formula were generally found to be quite inaccurate. This is especially the case for vibrational ground states in the adiabatic (long-lifetime) regime, which are states of paramount relevance in atom trapping and guiding. Expanding on earlier works in atom trapping [33] and Rydberg molecules [26], the non-perturbative methods tested in the present work can be generalized to problems with more than two adiabatic potentials with non-linear spatial dependence and variable mutual couplings.