Analytical Results for the Dynamics of Parabolic Level-Crossing Model

We study the dynamics of a two-level crossing model with a parabolic separation of the diabatic energies. The solutions are expressed in terms of the tri-confluent Heun equations --- the generalization of the confluent hypergeometric equations. We obtain analytical approximations for the state populations in terms of Airy and Bessel functions. Applicable expressions are derived for a large part of the parameter space. We also provide simple formulas which connect local solution in different time regimes. The validity of the analytical approximations is shown by comparing them to numerical simulations.


I. INTRODUCTION
Level crossing models are crucial for the understanding of non-adiabatic transitions in physics, chemistry and biology [1]. The best known and probably most widely studied levelcrossing model is the Landau-Zener model, which was formulated by Landau in 1932 for analyzing atomic collisions in both near-sudden [2] and near-adiabatic limits [3], and was subsequently solved by Zener using parabolic cylinder functions [4]. At around the same time, Stückelberg derived a sophisticated tunneling formula based on analytical continuation of the semi-classical WKB solutions across the Stokes lines [5], and Majorana derived the transition probability formula independently using integral representation of the survival amplitude in connection with the dynamics of a spin-1/2 in a time-varying magnetic field [6]. The Landau-Zener model assumes a constant coupling between bare states in the diabatic basis and a linearly varying separation of diabatic energies [7]. The tunneling probability exp(−ζ) in the Landau-Zener model only depends on a single dimensionless parameter ζ ≡ 2π f 2 /( |(F 1 − F 2 )V|), where f is the coupling matrix element in the diabatic basis, F 1 and F 2 are the slopes of the intersecting diabatic potential curves, and V is the velocity of the perturbation variable, e.g., the relative collision velocity [8]. Since Landau, Zener, Stückelberg and Majorana's pioneering works, the linear two-state model has been applied to atomic [9,10] and molecular collisions [11], atoms in intense laser fields [12,13], spin tunneling in molecular nano-magnets [14], tunneling of Bose-Einstein condensates in accelerated optical lattices [15], and optical tunneling in waveguide arrays [16].
Although the Landau-Zener model has achieved great success over the last century, there are indeed cases where the assumption of linear crossing between the diabatic states breaks down. To be more precise, one may employ the Stückelberg tunneling formula to the individual tunneling events [5,17], as long as any pair of Landau-Zener crossings are well-separated from each other. In other words, for cases in which the crossing points merge together as a result of external electric or magnetic fields, the Landau-Zener linearization fails, and the linear-dependence of diabatic energies should be replaced by a parabolic one [18]. The parabolic model was first introduced * Corresponding author. Email: yangbrookchen@yahoo.co.uk by Bikhovskii, Nikitin and Ovchinnikova in 1965 in the context of slow atomic collisions [19] and later re-evaluated by Delos and Thorson [20,21] and Crothers [22][23][24] in various limits. Two decades later, Shimshoni and Gefen incorporated environment-induced dissipation and dephasing into the parabolic model [25]. Suominen derived analytical approximations for the final state populations [26], and applied the results to the dynamics of cold atoms in magnetic traps [27]. In the same period, Zhu and Nakamura derived an exact formula for the scattering matrices in terms of a convergent infinite series, in which the coefficients satisfy a five-term recursion relation [28][29][30][31][32]. Nakamura and co-workers applied the results to laser assisted surface ion neutralization [33] and the laser-controlled photochromism in functional molecules [34]. Over the last decade, Lehto incorporated super-parabolic level-glancing effects into the parabolic model [35,36], and studied complete population inversion due to phase-jump couplings [37]. Zhang and co-workers described the population dynamics of driven dipolar molecules in the parabolic level-glancing model in terms of the confluent Heun functions [38]. Recently, the parabolic model has found applications in topological systems, including the interband tunneling of fermionic atoms near the merging transition of Dirac cones in tunable honeycomb optical lattices [39], and the interband tunneling of two-dimensional electrons near the topological transition in type II Weyl semimetals [40]. The dynamics of the parabolic level-crossing model, probably unknown to most physicists, may be written in terms of the tri-confluent Heun function [41], which is derived from the general Heun function [42] by the coalescence of three finite regular singularities with infinity. Unlike the parabolic cylinder function appearing in the Landau-Zener model, comprehensive information on the asymptotic behavior of confluent Heun functions, in general, does not exist. The main difficulty is due to the fact that the Heun functions do not possess any integral representations in terms of simpler special functions. Hence, Majorana's integral representation method is not applicable for determining the transition probability at infinity. Nevertheless, in the subsequent sections, we will show that valuable analytical approximations which only involves ordinary special functions can still be obtained in a large part of parameter space.
The paper is organized as follows: In Section II, we introduce the parabolic level-crossing model in the context of a two-level atom driven by a classical laser field, in which the arXiv:1906.11499v1 [quant-ph] 27 Jun 2019 laser detuning varies quadratically with time, and the Rabi frequency at resonance is time-independent. We express the final state population in terms of a single Stokes multiplier. In Section III, we derive analytical approximations for the transition amplitudes in both short-and long-time regimes. The validity of the analytical approximations is tested via comparisons with numerical results. We also discuss a way to connect analytical solutions in different time regimes. Finally, in Section IV, we conclude our studies, and discuss extensions of the results to other level-crossing models.

III. ANALYTICAL APPROXIMATIONS FOR THE TRANSITION AMPLITUDE
In the last section, we have shown that the dynamics of a two-level atom dipole-interacting with an off-resonant classical electric field with constant amplitude and parabolic detuning can be solved in terms of the tri-confluent Heun functions. We discussed the relationship between the final transition probability and the Stoke multipliers which connect asymptotic expansions of the tri-confluent Heun functions. However, due to mathematical difficulties involved, it is better to develop analytical approximations to the transition amplitudes, rather than to solve the connection problem rigorously.
To begin with, let us rewrite Eq. (6), the differential equation which governs U 1 , as where τ ≡ t + α/β. From Eq. (17), we see that the sign of α/β does not alter the nature of the equation. Hence, without loss of generality, we may assume that α/β is a positive number. Eq. (17) can also be rewritten as For large α/β, we compare the last two terms in Eq. (18). The short time regime is defined by β 2 τ 4 /16 α 2 τ 2 /8, i.e., |t| t * ≈ α/β; whereas the long time regime is defined by β 2 τ 4 /16 α 2 τ 2 /8, i.e., |t| t * . In the short time regime, the quadratic term is dominant in Eq. (18). Hence, the transition dynamics can be described by the Landau-Zener formula, and solved via the parabolic cylinder functions. But for small α/β, we should compare the third and the last terms in Eq. (18). Hence, the short time regime is defined by β 2 τ 4 /16 |βτ|/2, i.e., |τ| 2|β| −1/3 ; whereas the long time regime is defined by |τ| 2|β| −1/3 . In the short time regime, the linear term is dominant in Eq. (18). Hence, the transition dynamics is solved via the Airy functions, which is different from that obtained from the Landau-Zener formula. In the long time regime where the quartic term is dominant in Eq. (18), the transition dynamics can be solved via the Bessel functions.
In the following subsections, we discuss in detail the analytical approximations for the transition amplitude in different time regimes.
A. Dynamics in the short-time regime for small α/β In this subsection, we analyze the dynamics for the transition amplitude in the time regime |t| t * , so that β 2 τ 4 /16 α 2 τ 2 /8. We require |α| √ 2β 2/3 , so that α 2 τ 2 /8 |βτ|/2 is satisfied. We may rewrite Eq. (18) as where the right-hand side of Eq. (19) is regarded as a known function of τ. Let us denote z = ( iβ 2 ) and where . From Eq. (20), we see that the homogeneous equationŪ 1 = zŪ 1 is solved by the Airy functions Ai(z) and Bi(z). To be specific, we define the cubic root of iβ 2 as e ±i π 6 ( |β| 2 ) 1 3 for β = ±|β|, so that arg z ∈ ( π 6 , 7π 6 ) for β > 0 and arg z ∈ (− 7π 6 , − π 6 ) for β < 0. The particular solution of Eq. (20) is given by where W {Ai(ζ), Bi(ζ)} = 1/π is the Wronskian of Ai(ζ) and Bi(ζ). As the homogeneous solutionŪ 1 (z) is a linearly combination of Ai(z) and Bi(z), we need to evaluate integrals of the form z n y 2 dz and z n y 1 y 2 dz, where y refers to Ai(z) or Bi(z) respectively. After integration, Eq. (20) is solved by (details can be found in Appendix D) whereŪ 1 (z) ≡ c 1 Ai(z)+c 1 Bi(z) with c 1 and c 2 being arbitrary constants, and where a ≡ z + + z − and b ≡ z + z − . In Fig. 1, we compare the analytical approximation of the transition probability with the numerical solution. The result shows that the transition probability is well-described by Eqs. (22) and (23) in the short-time regime for small α/β.
B. Dynamics in the short-time regime for large α/β In this subsection, we analyze the dynamics for the transition amplitude in the time regime |t| t * , so that β 2 τ 4 /16 α 2 τ 2 /8. We require |α| √ 2β 2/3 , so that α 2 τ 2 /8 |βτ|/2 is satisfied. In this regime, it is better to use the original equation which governs U 1 . From Eq. (6), we obtain where the right-hand side of Eq. (24) is regarded as a known function of t. Let us denote z = e −iπ/4 α 1/2 t, then we obtain where a ≡ −i| f | 2 /α − 1/2 and z * ≡ e −iπ/4 α 3/2 /β. The homogeneous equationŪ 1 = (a + 1 4 z 2 )Ū 1 is solved by the parabolic cylinder functions U(a, ±z) and U(−a, ±iz). We choose the pair of linearly independent solutions as {U(a, z), U(−a, ∓iz}, where the plus and minus signs are for α > 0 and α < 0 re-spectively. The particular solution of Eq. (25) is given by where W {U(a, ζ), U(−a, ∓iζ)} = ±ie ∓iπ( a 2 + 1 4 ) is the Wronskian of U(a, ζ) and U(−a, ∓iζ). As the homogeneous solution U 1 (z) is a linearly combination of U(a, z) and U(−a, ∓iz), we need to evaluate integrals of the form z n y 2 dz and z n y 1 y 2 dz, where y refers to U(a, ζ) or U(−a, ∓iζ) respectively. After integration, Eq. (25) can still be solved by (details can be found in Appendix D)  Fig. 2b, we compare the transition probability obtained from |Ū 1 (z)| 2 with the exact solution, and the result shows that the final transition probability is well-described by the Landau-Zener formula in the wholetime range for large α/β.

D. Analytical Approximations to the Connection Problem
In the last subsections, we have studied the dynamics for the transition amplitude in both short-and long-time regimes. In this subsection, we focus on the problem of connecting different local solutions of the transition dynamics. As we showed in Fig. 2b, the transition dynamics for large α/β can be welldescribed by the Landau-Zener formula. Hence, we concentrate on the transition dynamics for small α/β. Without loss of generality, we only consider the parabolic glancing case with α = 0. Then Eq. (17) may be written as where λ ≡ β/4. Substitution of U 1 ≡ ρe iS into Eq. (33) yields the following set of differential equations 1 ρ where the initial conditions are ρ(−∞) = 0 andρ(−∞) = 1. In the long time limit, we expect that the transition probability ρ 2 = |U 1 | 2 converges toward a stationary value. Hence, we assumeρ ≈ 0 andṠ 2 = | f | 2 + λ 2 t 4 in Eq. (34a), which yields where λ 2 ρ f t −5 in the long-time limit. We will use this solution in later discussions.

IV. CONCLUSION
To summarize, we studied the transition dynamics of the parabolic model -a two-state system subject to a quadratically detuning over an infinite time interval. The solutions are expressed in terms of the tri-confluent Heun functions, which are the generalizations of the conventional confluent hypergeometric functions. Instead of rigorously solving the Stokes multipliers which connect the asymptotic solutions of the transition amplitudes, we derived concise analytical approximations to the transition amplitudes in both short-and long-time regimes, and provided practical formulas for connecting local solutions in different regimes. We gave applicable estimation of the critical times that separate different time regimes. The transition dynamics is shown to be well-described by the ana-lytical formulas in the whole-time range by comparison with exact results.
In future works, we would like to extend our study to superlinear model, in which the laser detuning is a polynomial in time with cubic or higher degrees, and the sub-linear model in which the laser detuning is a rational function in time. We also want to determine the critical times that separate the short-and long-time regimes rigorously, instead of only giving a rough estimate of the magnitude.

ACKNOWLEDGMENTS
The Authors would like to thank the Science and Technology Development Fund of the Macau SAR for providing support, FDCT 023/2017/A1.

Appendix A: Equations of Motion for a Two-level Atom
Dipole-Interacting with a Classical Driving Field In this appendix, we provide the background knowledge of the two-level atom description of light-matter integrations for readers. When the two target atomic levels are nearly resonant with the driving field, while on the same time the other atomic levels are detuned far off resonance, we may regard the system as two discrete non-degenerate states, e.g., the ground state |g and the first excited state |e , then the Hamiltonian of the atom may be written asĤ atom = E e |e e| + E g |g g|. For the case when the atom interacts with an external electric field under dipole approximation, the interaction Hamiltonian becomeŝ H int = −d · E(r cm ), wherer cm is the position operator for the center of mass. When the De Broglie wavelength of the atom is small compared to the interatomic spacing, the center of mass position may be treated classically [45].
In the following, we discuss the coherent excitation of the two-level atom under the driving of an external electric field. Let us denote the state of the two-level atom as |ψ(t) = c e (t)|e + c g (t)|g , the Schrödinger equations for the two wave amplitudes have the form whered ≡ qx is the transition electric dipole moment of the atom. For atom that possesses inversion symmetry, the energy eigenstates are symmetric or anti-symmetric, and hence the expectation values of the dipole moment vanishes, e|d|e = g|d|g = 0. Let us denote the energies of the two atomic levels as E e = 2 ω 0 and E g = − 2 ω 0 , where ω 0 ≡ E e − E g is the energy difference between the two atomic levels. For a monochromatic wave, the electric field may be written as where E 0 is the complex electric field amplitude at the position r cm , and e p is the polarization vector of the incident electric field. If we denote the transition dipole moment e|d|g as d ge , the Schrödinger equations for the two wave amplitudes may be simplified as In the reference frame rotating about the z-axis with frequency ω(t), the state of the two-level atom is |ψ (t) = exp( i τ 0 ω(τ)dtS z )|ψ(t) . Then the relationships between the wave amplitudes in the rotated and unrotated frame are where ∆ ≡ ω − ω 0 is the laser detuning. Applying the rotating wave approximation, we neglect the fast-oscillating terms like exp{±2i t 0 ω(τ)dτ} in the Schrödinger equations, and obtain where f ≡ − E * 0 2 d ge · e * p is the Rabi frequency for the transition dipole moment. For the special case that ω = ω 0 , the laser detuning ∆ vanishes, and the coupled-mode equations become iȧ e = f a g and iȧ g = f * a e , which are equivalent tö a e −ḟ fȧ e + | f | 2 a e = 0,ä g −ḟ * When the atom is initially at the ground state, we have c g (0) = 1 and c e (0) = 0, or equivalently a g (0) = 1 and a e (0) = 0. For the special case that the Rabi frequency is time-independent, Eq. (A5) is solved by a e = sin f t and a g = cos f t, or equivalently c e = e − i 2 t 0 ω(τ)dτ sin f t and c g = e i 2 t 0 ω(τ)dτ cos f t in the lab frame, which results in the occupation probabilities |c e | 2 = sin 2 f t and |c g | 2 = cos 2 f t.

Integrals of products of parabolic cylinder functions
In this sub-appendix, we evaluate indefinite integrals of the form I n ≡ z n y 2 dz and J n ≡ z n y 1 y 2 dz in terms of parabolic cylinder functions and their first derivatives, where y 1 and y 2 are both solutions of the parabolic cylinder differential equation y = (a + 1 4 z 2 )y. If we write I n ≡ 1 2 R n − zR n y 2 − R n yy + R n y 2 , we obtain R n − (4a + z 2 )R n − zR = 2z n . We may expand R n in a series as ∞ k=0 c k z k , where the coefficients c k satisfies the following three-term recursion relation (k + 2)(k + 1)kc k+2 − 4akc k − (k − 1)c k−2 =        2, k = n + 1; 0, otherwise.

Integrals of products of Bessel functions
In this sub-appendix, we evaluate indefinite integrals of the form I n ≡ τ n y 2 dτ and J n ≡ τ n y 1 y 2 dτ in terms of Bessel functions and their first derivatives, where y 1 and y 2 are both solutions of the differential equation y = −λ 2 τ 4 y. Using the ansatz I n = 1 2 R n + λ 2 τ 4 R n y 2 −R n yy +R n y 2 , we obtain R n + 4λ 2 τ 4 R n + 8λ 2 τ 3 R n = 2τ n , where R n satisfies the following recursion relation R n = 2τ n+3 − 4λ 2 (n + 5)R n+6 (n + 1)(n + 2)(n + 3) .