Analytical Results for the Dynamics of Parabolic Level-Crossing Model

TL;DR

This paper derives analytical solutions for a two-level parabolic crossing model using special functions, providing approximations valid over extensive parameter ranges and validated by numerical comparisons.

Contribution

It introduces new analytical approximations for the model's dynamics using Heun, Airy, and Bessel functions, expanding understanding of parabolic level-crossing phenomena.

Findings

Analytical solutions expressed via special functions.

Approximate formulas valid across large parameter spaces.

Validation through numerical simulations.

Abstract

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.

Equations

$$i\frac{d{a}_1}{dt}=-\frac{\Delta }{2}{a}_1+f{a}_2,i\frac{d{a}_2}{dt}=f^{\ast }{a}_1+\frac{\Delta }{2}{a}_2,$$

$$i\frac{d{C}_1}{dt}=f{e}^{-i{\int }_0^t\Delta ds}{C}_2,i\frac{d{C}_2}{dt}=f^{\ast }{e}^{i{\int }_0^t\Delta ds}{C}_1,$$

$$\frac{d^2{C}_1}{d{t}^2}+\left(i\Delta -\frac{\dot{f}}{f}\right)\frac{d{C}_1}{dt}+|f{|}^2{C}_1=0,$$

$$\frac{d^2{C}_2}{d{t}^2}-\left(i\Delta +\frac{{\dot{f}}^{\ast }}{f^{\ast }}\right)\frac{d{C}_2}{dt}+|f{|}^2{C}_2=0,$$

$$\frac{d^2{U}_1}{d{t}^2}+J(t)U_1=0,$$

$$J(t)=|f{|}^2-\frac{1}{2}\frac{d}{dt}\left(i\Delta -\frac{\dot{f}}{f}\right)-\frac{1}{4}{\left(i\Delta -\frac{\dot{f}}{f}\right)}^2.$$

$$\frac{d^2{U}_1}{d{t}^2}+\left(|f{|}^2-\frac{i\alpha }{2}+\frac{{\alpha }^2{t}^2}{4}\right)U_1=0.$$

$$\frac{d^2{U}_1}{d{t}^2}+\left(|f{|}^2-\frac{i\alpha }{2}-\frac{i\beta t}{2}+\frac{{\alpha }^2{t}^2}{4}+\frac{\alpha \beta t^3}{4}+\frac{{\beta }^2{t}^4}{16}\right)U_1=0.$$

$$\frac{d^2{U}_1}{d{z}^2}+\left(\mu -\frac{{\xi }^2}{4}+\nu z-\frac{3}{2}\xi z^2-\frac{9}{4}{z}^4\right)U_1=0,$$

$$\mu$$

$$\frac{d^2{V}_1}{d{z}^2}+(\xi +3z^2)\frac{d{V}_1}{dz}+(\mu +(\nu -3)z)V_1=0,$$

$$T_1(\mu ,\nu ,\xi ;z)=e^{-\frac{1}{2}(z^3+\xi z)}{z}^{\frac{\nu }{3}-1}\sum _{k\ge 0}{a}_k(\mu ,\nu ,\xi )z^{-k},$$

$$T_2(\mu ,\nu ,\xi ;z)$$

$$=e^{\frac{1}{2}(z^3+\xi z)}{z}^{-\frac{\nu }{3}-1}\sum _{k\ge 0}(-1{)}^k{a}_k(\mu ,-\nu ,\xi )z^{-k},$$

$$3(k+3)a_{k+3}+\mu a_{k+2}+\xi (k+2-\nu /3)a_{k+1}$$

$$+(k+1-\nu /3)(k+2-\nu /3)a_k=0.$$

$$U_1(t\to -\infty )$$

$$\approx A_1(\beta /6{)}^{-\frac{2}{3}}{e}^{-i(\frac{\beta }{12}{t}^3+\frac{\alpha }{4}{t}^2-\frac{1}{6}\frac{{\alpha }^3}{{\beta }^2}-\frac{\pi }{3})}{t}^{-2}.$$

$$C_1(t\to -\infty )$$

$${\dot{C}}_1(t\to -\infty )$$

$$T_2(\mu ,\nu ,\xi ;z)=\mathcal{L}{T}_1(\mu ,\nu ,\xi ;z)-\mathcal{L}C(\mu ,\nu ,\xi )T_1(\mu ,\nu ,\xi ;z),$$

$$U_1(t\to \infty )$$

$$-C({\omega }^4\mu ,-\nu ,{\omega }^2\xi )e^{i(\frac{\beta }{12}{t}^3+\frac{\alpha }{4}{t}^2-\frac{1}{6}\frac{{\alpha }^3}{{\beta }^2})}].$$

$$C_1(t\to \infty )\approx -A_{1}C({\omega }^4\mu ,-\nu ,{\omega }^2\xi )e^{-i\frac{{\alpha }^3}{6{\beta }^2}},$$

$$\frac{d^2{U}_1}{d{\tau }^2}+\left[|f{|}^2-\frac{i\beta \tau }{2}+\frac{{\beta }^2}{16}{\left({\tau }^2-\frac{{\alpha }^2}{{\beta }^2}\right)}^2\right]U_1=0,$$

$$\frac{d^2{U}_1}{d{\tau }^2}+\left(|f{|}^2+\frac{{\alpha }^4}{16{\beta }^2}-\frac{i\beta \tau }{2}-\frac{{\alpha }^2{\tau }^2}{8}+\frac{{\beta }^2{\tau }^4}{16}\right)U_1=0.$$

$$\frac{d^2{U}_1}{d{\tau }^2}+\left(|f{|}^2+\frac{{\alpha }^4}{16{\beta }^2}-\frac{i\beta \tau }{2}\right)U_1=\left(\frac{{\alpha }^2{\tau }^2}{8}-\frac{{\beta }^2{\tau }^4}{16}\right)U_1,$$

$$\frac{d^2{U}_1}{d{z}^2}-z{U}_1$$

$$\equiv {\pi }_4(z)U_1,$$

$$U_{1p}(z)$$

$$+\mathrm{Bi}(z){\int }_{z_0}^z\frac{\mathrm{Ai}(\zeta ){\pi }_4(\zeta ){\bar{U}}_1(\zeta )}{\mathcal{W}\{\mathrm{Ai}(\zeta ),\mathrm{Bi}(\zeta )\}}d\zeta ,$$

$$U_1(z)\approx \left(1-\frac{1}{2}\frac{dG(z)}{dz}\right){\bar{U}}_1(z)+G(z){\bar{U}}_1^{\prime }(z),$$

$$G(z)$$

$$+\left(\frac{4}{9}-\frac{2ab}{3}\right)z+b^2-\frac{6a}{7}]+\frac{{\alpha }^4}{16{\beta }^2},$$

$$\frac{d^2{U}_1}{d{t}^2}+\left(|f{|}^2-\frac{i\alpha }{2}+\frac{{\alpha }^2{t}^2}{4}\right)U_1\approx \left(\frac{i\beta t}{2}-\frac{\alpha \beta t^3}{4}\right)U_1,$$

$$\frac{d^2{U}_1}{d{z}^2}-\left(a+\frac{z^2}{4}\right)U_1\approx \left(-\frac{z}{2z^{\ast }}+\frac{z^3}{4z^{\ast }}\right)U_1\equiv {\pi }_3(z)U_1,$$

Full text

I introduction

Level crossing models are crucial for the understanding of non-adiabatic transitions in physics, chemistry and biology nakamura2012nonadiabatic . The best known and probably most widely studied level-crossing model is the Landau-Zener model, which was formulated by Landau in 1932 for analyzing atomic collisions in both near-sudden landau1932theoryI and near-adiabatic limits landau1932theory , and was subsequently solved by Zener using parabolic cylinder functions zener1932non . 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 stuckelberg1932theory , 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 majorana1932atomi . The Landau-Zener model assumes a constant coupling between bare states in the diabatic basis and a linearly varying separation of diabatic energies wittig2005landau . The tunneling probability $\exp(-\zeta)$ in the Landau-Zener model only depends on a single dimensionless parameter $\zeta\equiv 2\pi f^{2}/(\hbar|(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 di2005majorana . Since Landau, Zener, Stückelberg and Majorana’s pioneering works, the linear two-state model has been applied to atomic smirnov2003physics ; nikitin2012theory and molecular collisions child2010molecular , atoms in intense laser fields delone1985atoms ; kazantsev1990mechanical , spin tunneling in molecular nano-magnets wernsdorfer1999quantum , tunneling of Bose-Einstein condensates in accelerated optical lattices morsch2001bloch , and optical tunneling in waveguide arrays khomeriki2005nonadiabatic .

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 stuckelberg1932theory ; shevchenko2010landau , 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 garraway1995wave . The parabolic model was first introduced by Bikhovskii, Nikitin and Ovchinnikova in 1965 in the context of slow atomic collisions bikhovskii1965probability and later re-evaluated by Delos and Thorson delos1972solution ; delos1972studies and Crothers crothers1975stueckelberg ; crothers1976perturbed ; crothers1977stueckelberg in various limits. Two decades later, Shimshoni and Gefen incorporated environment-induced dissipation and dephasing into the parabolic model shimshoni1991onset . Suominen derived analytical approximations for the final state populations suominen1992parabolic , and applied the results to the dynamics of cold atoms in magnetic traps pakarinen2000atomic . 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 zhu1992two ; zhu1992twoII ; zhu1993two ; zhu1994two ; zhu2001nonadiabatic . Nakamura and co-workers applied the results to laser assisted surface ion neutralization teranishi1997semiclassical and the laser-controlled photochromism in functional molecules tamura2006laser . Over the last decade, Lehto incorporated super-parabolic level-glancing effects into the parabolic model lehto2012superparabolic ; lehto2013zhu , and studied complete population inversion due to phase-jump couplings lehto2016two . 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 zhang2016analytic . 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 fuchs2012interband , and the interband tunneling of two-dimensional electrons near the topological transition in type II Weyl semimetals malla2019high .

The dynamics of the parabolic level-crossing model, probably unknown to most physicists, may be written in terms of the tri-confluent Heun function ronveaux1995heun , which is derived from the general Heun function heun1888theorie 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 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.

II The Parabolic Model

Within the rotating-wave approximation, the wave amplitudes of a two-level atom dipole-interacting with a classical electric field are governed by the following coupled equations (details are provided in Appendix A)

[TABLE]

where $\Delta$ is the laser detuning, and $f$ is the Rabi frequency at resonance. After applying the transformations $a_{1}=C_{1}e^{\frac{i}{2}\int_{0}^{t}\Delta ds}$ and $a_{2}=C_{2}e^{-\frac{i}{2}\int_{0}^{t}\Delta ds}$, the coupled equations becomes

[TABLE]

or equivalently

[TABLE]

where $C_{1}$ and $C_{2}$ satisfy the normalization condition $|C_{1}|^{2}+|C_{2}|^{2}=1$. The problem of non-adiabatic transition is to determine the transition probability $|C_{1}(\infty)|^{2}$, subjected to the conditions $|C_{1}(-\infty)|=0$ and $|C_{2}(-\infty)|=1$. Using the change of variable $C_{1}=U_{1}\exp\{-\frac{1}{2}\int_{0}^{t}pds\}$, where $p\equiv i\Delta-\dot{f}/f$, we obtain the Schrödinger form of Eq. (3a)

[TABLE]

In the conventional Landau-Zener model, the laser detuning varies linearly with time, $\nu\equiv\alpha t$, and the the Rabi frequency at resonance is time-independent, which yields

[TABLE]

Eq. (5) becomes the parabolic cylinder equation $U_{1}^{\prime\prime}+(n+\frac{1}{2}-\frac{1}{4}z^{2})U_{1}=0$ via the change of variables $z\equiv e^{-i\pi/4}\alpha^{1/2}t$ and $n\equiv i|f|^{2}/\alpha$. The transition probability $|C_{1}(\infty)|^{2}$ can be obtained from the asymptotic expansions of the parabolic cylinder function (details are provided in Appendix B).

In contrast to the conventional Landau-Zener model, the laser detuning in the parabolic model varies quadratically with time, $\Delta\equiv\alpha t+\frac{1}{2}\beta t^{2}$, and the Rabi frequency at resonance is time-independent, which yields

[TABLE]

Eq. (6) may be transformed into the canonical form of tri-confluent Heun equation. Let us perform the transformation $z\equiv h^{-1}(t+\alpha/\beta)$ with $h^{6}=-\frac{9}{4}\lambda^{-2}$ and $\lambda\equiv\beta/4$, then Eq. (6) becomes the second canonical form of the tri-confluent Heun equation (THE2 equation) ronveaux1995heun

[TABLE]

which may be transformed into the first canonical form of the tri-confluent Heun equation (THE1 equation) ronveaux1995heun

[TABLE]

via the transformation $U_{1}=e^{-\frac{1}{2}(z^{3}+\xi z)}V_{1}$. To be precise and for later convenience, we define $h\equiv e^{i\pi/6}(3/2\lambda)^{1/3}$ for $\lambda>0$, and define $\mathcal{L}U_{1}(\mu,\nu,\xi;z)=U_{1}(\omega^{4}\mu,\omega^{3}\nu,\omega^{2}\xi;\omega z)$ for $\omega=e^{i\pi/3}$.

Similar to Ziner’s approach to the linear level-crossing model zener1932non , the transition probability may be derived from the asymptotic expansions of the tri-confluent Heun function at different sectors in the complex plane. The THE2 equation has two independent solutions $T_{1}(\mu,\nu,\xi;z)$ and $T_{2}(\mu,\nu,\xi;z)$, where $T_{1}(\mu,\nu,\xi;z)$ has the following asymptotic expansion in the sector $|\arg{z}|<\frac{\pi}{2}$ ronveaux1995heun

[TABLE]

and $T_{2}(\mu,\nu,\xi;z)$ has the following asymptotic expansion in the sector $\frac{\pi}{2}<\arg z<\frac{3\pi}{2}$ ronveaux1995heun

[TABLE]

where $a_{0}(\mu,\nu,\xi)=1$, $a_{1}(\mu,\nu,\xi)=-\mu/3$, $a_{2}(\mu,\nu,\xi)=\frac{1}{18}(\mu^{2}+\xi(\nu-3))$ and $a_{k}(\mu,\nu,\xi)$ satisfies the following four-term recursion relation ronveaux1995heun

[TABLE]

For $t\rightarrow-\infty$, we have $z=-h^{-1}|t+\alpha/\beta|$ and $\arg z=\frac{5\pi}{6}$. Hence, we may use the solution $T_{2}(\mu,\nu,\xi;z)=T_{1}(\mu,-\nu,\xi;-z)$, so that

[TABLE]

Using the relation $C_{1}=U_{1}\exp\{-\frac{i}{2}\int_{0}^{t}(\alpha s+\frac{\beta}{2}s^{2})ds\}$, we obtain

[TABLE]

The constant $A_{1}$ is determined by $|\dot{C}_{1}|=|f\dot{C}_{2}|=|f|$, which yields $A_{1}=\frac{|f|}{3}(\frac{\beta}{6})^{-\frac{1}{3}}$. The large $|t|$ solutions Eqs. (13a) and (13b) can also be obtained from the method of direct integration (details are provided in Appendix C). For $t\rightarrow\infty$, we have $z=h^{-1}(t+\alpha/\beta)$ and $\arg z=-\pi/6$. Hence, Eq. (II), the asymptotic expansion for $T_{2}(\mu,\nu,\xi;z)$ may not be used. In order to evaluate $T_{2}(\mu,\nu,\xi;z)$ for $t\rightarrow\infty$, we have to use the connection formula ronveaux1995heun

[TABLE]

where $\mathcal{L}C(\mu,\nu,\xi)\equiv C(\omega^{4}\mu,-\nu,\omega^{2}\xi)$. $C(\mu,\nu,\xi)$ is the Stokes multiplier which connects the asymptotic expansions of $U_{1}(z)$ at different sectors, and is an entire function of $\mu$, $\nu$ and $\xi$ sibuya1975global . For $t\rightarrow\infty$, we have $\arg(h^{-1}(t+\alpha/\beta))=-\frac{\pi}{6}$ and $\arg(\omega h^{-1}(t+\alpha/\beta))=\frac{\pi}{6}$. Hence, we may use the asymptotic expansion of $T_{1}(\mu,\nu,\xi;z)$ and obtain

[TABLE]

Using the relation $C_{1}=U_{1}e^{-i(\frac{\beta}{12}t^{3}+\frac{\alpha}{4}t^{2})}$, we obtain

[TABLE]

which yields $|C_{1}(\infty)|^{2}=|A_{1}C(\omega^{4}\mu,-\nu,\omega^{2}\xi)|^{2}$. Hence, the final transition probability $|C_{1}(\infty)|^{2}$ depends only on the Stokes multiplier $C(\mu,\nu,\xi)$. However, it is in general not an easy task to obtain exact formulas for the Stokes multipliers. Although Zhu and Nakamura derived an exact formula for the Stokes multipliers in terms of a sophisticated infinite series generated by a five-term recursion relation zhu1992stokes , a compact formula for the final transition probability which is similar to the Landau-Zener formula has not yet existed. In the following section, we derive concise and explicit expressions for the transition dynamics in the parabolic model, and provide an analytical approximation for connecting solutions of the transition amplitude in different time regimes.

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

[TABLE]

where $\tau\equiv t+\alpha/\beta$. From Eq. (17), we see that the sign of $\alpha/\beta$ does not alter the nature of the equation. Hence, without loss of generality, we may assume that $\alpha/\beta$ is a positive number. Eq. (17) can also be rewritten as

[TABLE]

For large $\alpha/\beta$, we compare the last two terms in Eq. (18). The short time regime is defined by $\beta^{2}\tau^{4}/16\lessapprox\alpha^{2}\tau^{2}/8$, i.e., $|t|\lessapprox t^{*}\approx\alpha/\beta$; whereas the long time regime is defined by $\beta^{2}\tau^{4}/16\gtrapprox\alpha^{2}\tau^{2}/8$, i.e., $|t|\gtrapprox 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 $\alpha/\beta$, we should compare the third and the last terms in Eq. (18). Hence, the short time regime is defined by $\beta^{2}\tau^{4}/16\lessapprox|\beta\tau|/2$, i.e., $|\tau|\lessapprox 2|\beta|^{-1/3}$; whereas the long time regime is defined by $|\tau|\gtrapprox 2|\beta|^{-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.

III.1 Dynamics in the short-time regime for small $\alpha/\beta$

In this subsection, we analyze the dynamics for the transition amplitude in the time regime $|t|\lessapprox t^{*}$, so that $\beta^{2}\tau^{4}/16\lessapprox\alpha^{2}\tau^{2}/8$. We require $|\alpha|\lessapprox\sqrt{2}\beta^{2/3}$, so that $\alpha^{2}\tau^{2}/8\lessapprox|\beta\tau|/2$ is satisfied. We may rewrite Eq. (18) as

[TABLE]

where the right-hand side of Eq. (19) is regarded as a known function of $\tau$. Let us denote $z=(\frac{i\beta}{2})^{\frac{1}{3}}[\tau+\frac{2i}{\beta}(|f|^{2}+\frac{\alpha^{4}}{16\beta^{2}})]$, then we obtain $\tau=(\frac{i\beta}{2})^{-\frac{1}{3}}z-\frac{2i}{\beta}(|f|^{2}+\frac{\alpha^{4}}{16\beta^{2}})$ and

[TABLE]

where $z_{\pm}\equiv(\frac{i\beta}{2})^{1/3}(2i\gamma\pm 1)\frac{\alpha}{\beta}$ and $\gamma\equiv\frac{1}{\alpha}(|f|^{2}+\frac{\alpha^{4}}{16\beta^{2}})$. From Eq. (III.1), we see that the homogeneous equation $\bar{U}^{\prime\prime}_{1}=z\bar{U}_{1}$ is solved by the Airy functions $\operatorname{Ai}(z)$ and $\operatorname{Bi}(z)$. To be specific, we define the cubic root of $\frac{i\beta}{2}$ as $e^{\pm i\frac{\pi}{6}}(\frac{|\beta|}{2})^{\frac{1}{3}}$ for $\beta=\pm|\beta|$, so that $\arg z\in(\frac{\pi}{6},\frac{7\pi}{6})$ for $\beta>0$ and $\arg z\in(-\frac{7\pi}{6},-\frac{\pi}{6})$ for $\beta<0$. The particular solution of Eq. (III.1) is given by

[TABLE]

where $\mathscr{W}\{\operatorname{Ai}(\zeta),\operatorname{Bi}(\zeta)\}=1/\pi$ is the Wronskian of $\operatorname{Ai}(\zeta)$ and $\operatorname{Bi}(\zeta)$. As the homogeneous solution $\bar{U}_{1}(z)$ is a linearly combination of $\operatorname{Ai}(z)$ and $\operatorname{Bi}(z)$, we need to evaluate integrals of the form $\int z^{n}y^{2}dz$ and $\int z^{n}y_{1}y_{2}dz$, where $y$ refers to $\operatorname{Ai}(z)$ or $\operatorname{Bi}(z)$ respectively. After integration, Eq. (III.1) is solved by (details can be found in Appendix D)

[TABLE]

where $\bar{U}_{1}(z)\equiv c_{1}\operatorname{Ai}(z)+c_{1}\operatorname{Bi}(z)$ with $c_{1}$ and $c_{2}$ being arbitrary constants, and

[TABLE]

where $a\equiv z_{+}+z_{-}$ and $b\equiv 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 (III.1) in the short-time regime for small $\alpha/\beta$.

III.2 Dynamics in the short-time regime for large $\alpha/\beta$

In this subsection, we analyze the dynamics for the transition amplitude in the time regime $|t|\lessapprox t^{*}$, so that $\beta^{2}\tau^{4}/16\lessapprox\alpha^{2}\tau^{2}/8$. We require $|\alpha|\gtrapprox\sqrt{2}\beta^{2/3}$, so that $\alpha^{2}\tau^{2}/8\gtrapprox|\beta\tau|/2$ is satisfied. In this regime, it is better to use the original equation which governs $U_{1}$. From Eq. (6), we obtain

[TABLE]

where the right-hand side of Eq. (24) is regarded as a known function of $t$. Let us denote $z=e^{-i\pi/4}\alpha^{1/2}t$, then we obtain

[TABLE]

where $a\equiv-i|f|^{2}/\alpha-1/2$ and $z^{*}\equiv e^{-i\pi/4}\alpha^{3/2}/\beta$. The homogeneous equation $\bar{U}^{\prime\prime}_{1}=(a+\frac{1}{4}z^{2})\bar{U}_{1}$ is solved by the parabolic cylinder functions $U(a,\pm z)$ and $U(-a,\pm iz)$. We choose the pair of linearly independent solutions as $\{U(a,z),U(-a,\mp iz\}$, where the plus and minus signs are for $\alpha>0$ and $\alpha<0$ respectively. The particular solution of Eq. (25) is given by

[TABLE]

where $\mathscr{W}\{U(a,\zeta),U(-a,\mp i\zeta)\}=\pm ie^{\mp i\pi(\frac{a}{2}+\frac{1}{4})}$ is the Wronskian of $U(a,\zeta)$ and $U(-a,\mp i\zeta)$. As the homogeneous solution $\bar{U}_{1}(z)$ is a linearly combination of $U(a,z)$ and $U(-a,\mp iz)$, we need to evaluate integrals of the form $\int z^{n}y^{2}dz$ and $\int z^{n}y_{1}y_{2}dz$, where $y$ refers to $U(a,\zeta)$ or $U(-a,\mp i\zeta)$ respectively. After integration, Eq. (25) can still be solved by (details can be found in Appendix D)

[TABLE]

where $\bar{U}_{1}(z)\equiv c_{1}U(a,z)+c_{2}U(-a,\mp iz)$ with $c_{1}$ and $c_{2}$ being arbitrary constants, and $G(z)=-(6+8a-z^{2})/(6z^{*})$. As we can see from Fig. 2a, the Landau-Zener solution $|\bar{U}_{1}(z)|^{2}=|c_{1}U(a,z)+c_{2}U(-a,\mp iz)|^{2}$ coincides with the exact result of the transition probability obtained from Eq. (17), but the analytical approximation Eq. (27) does not provide better result than the Landau-Zener solution. In Fig. 2b, we compare the transition probability obtained from $|\bar{U}_{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 whole-time range for large $\alpha/\beta$.

III.3 Dynamics in the long-time regime

In this subsection, we analyze the dynamics for the transition amplitude in the time regime $|t|\gtrapprox t^{*}$, so that $\beta^{2}\tau^{4}/16\gtrapprox\alpha^{2}\tau^{2}/8$. We may rewrite Eq. (18) as

[TABLE]

where the right-hand side of Eq. (III.3) is regarded as a known function of $\tau$. The homogeneous equation $\bar{U}_{1}^{\prime\prime}=-\frac{1}{16}\beta^{2}\tau^{4}\bar{U}_{1}$ is solved by the Bessel functions of order $1/6$, that is, $w_{1}\equiv\sqrt{\tau}J_{1/6}(\beta\tau^{3}/12)$ and $w_{2}\equiv\sqrt{\tau}J_{-1/6}(\beta\tau^{3}/12)$. The particular solution of Eq. (III.3) is given by

[TABLE]

where $\mathscr{W}\{w_{1}(\tau),w_{2}(\tau)\}=-3/\pi$ is the Wronskian of $w_{1}(\tau)$ and $w_{2}(\tau)$. As the homogeneous solution $\bar{U}_{1}$ is a linearly combination of $w_{1}$ and $w_{2}$, we need to evaluate integrals of the form $I_{n}\equiv\int\tau^{n}w^{2}d\tau$ and $J_{n}\equiv\int\tau^{n}w_{1}w_{2}d\tau$, where $w$ refers to $w_{1}$ or $w_{2}$ respectively. After integration, Eq (III.3) is solved by (details can be found in Appendix D)

[TABLE]

where $\bar{U}_{1}(\tau)=c_{1}\sqrt{\tau}J_{1/6}\left(\beta\tau^{3}/12\right)+c_{2}\sqrt{\tau}J_{-1/6}\left(\beta\tau^{3}/12\right)$ with $c_{1}$ and $c_{2}$ being arbitrary constants, and

[TABLE]

Here ${}_{p}F_{q}(a_{1},\cdots,a_{p};b_{1},\cdots,b_{q};z)$ is the generalized hypergeometric function of order $p$, $q$, and $G(\tau)$ has the following asymptotic expansion in the limit $\tau\rightarrow\infty$

[TABLE]

As we can see from Fig. 3, the analytical approximation of the transition probability $|U_{1}|^{2}$, which is determined by Eq. (30) with $G(\tau)\approx\alpha^{2}/\beta^{2}\tau^{-1}$, agrees with the exact result in the long-time regime.

III.4 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 $\alpha/\beta$ can be well-described by the Landau-Zener formula. Hence, we concentrate on the transition dynamics for small $\alpha/\beta$. Without loss of generality, we only consider the parabolic glancing case with $\alpha=0$. Then Eq. (17) may be written as

[TABLE]

where $\lambda\equiv\beta/4$. Substitution of $U_{1}\equiv\rho e^{iS}$ into Eq. (33) yields the following set of differential equations

[TABLE]

where the initial conditions are $\rho(-\infty)=0$ and $\dot{\rho}(-\infty)=1$. In the long time limit, we expect that the transition probability $\rho^{2}=|U_{1}|^{2}$ converges toward a stationary value. Hence, we assume $\ddot{\rho}\approx 0$ and $\dot{S}^{2}=|f|^{2}+\lambda^{2}t^{4}$ in Eq. (34a), which yields $S(t)\approx\lambda t^{3}/3+S_{0}$ and $\frac{d}{dt}\left(\rho^{2}\sqrt{|f|^{2}+\lambda^{2}t^{4}}\right)=2\lambda t\rho^{2}$. A direct computation yields

[TABLE]

where $\rho\approx\rho_{f}(1-\frac{|f|^{2}}{2\lambda^{2}}t^{-4})$ and $\dot{\rho}\approx-\frac{2|f|^{2}}{\lambda^{2}}\rho_{f}t^{-5}$ in the long-time limit. We will use this solution in later discussions.

As we have discussed in the last subsections, the transition dynamics for the parabolic-glancing case can be well-approximated by the Airy and Bessel functions in the short- and long-time regimes. Hence, we only need to connect the local solutions at two critical times $t_{1}$ and $t_{2}$ — for $t\in(-\infty,t_{1}]$, we have $U_{1}^{(I)}=c_{1}^{(I)}w_{1}(t)+c_{2}^{(I)}w_{2}(t)$; for $t\in[t_{1},t_{2}]$, we have $U_{1}^{(II)}=c_{1}^{(II)}\operatorname{Ai}(z)+c_{2}^{(II)}\operatorname{Bi}(z)$ with $z\equiv(\frac{i\beta}{2})^{1/3}(t+\frac{2i}{\beta}|f|^{2})$; and for $t\in[t_{2},\infty)$, we have $U_{1}^{(III)}=c_{1}^{(III)}w_{1}(t)+c_{2}^{(III)}w_{2}(t)$, where $w_{1}\equiv\sqrt{t}J_{1/6}(\beta t^{3}/12)$ and $w_{2}\equiv\sqrt{t}J_{-1/6}(\beta t^{3}/12)$. Then the connection formulas for the transition amplitudes are

[TABLE]

A direct computation yields

[TABLE]

where $\mathscr{W}\{\operatorname{Ai}(z(t)),\operatorname{Bi}(z(t))\}=(i\beta/2)^{1/3}/\pi$ is the Wronskian of $\operatorname{Ai}(z(t))$ and $\operatorname{Bi}(z(t))$, $\dot{w}_{1}(t)=\lambda t^{5/2}J_{-5/6}(\lambda t^{3}/3)$, $\dot{w}_{2}(t)=\lambda t^{5/2}J_{5/6}(\lambda t^{3}/3)$, and $\mathscr{W}\{w_{1}(t),w_{2}(t)\}=-3/\pi$ is the Wronskian of $w_{1}(t)$ and $w_{2}(t)$. Hence, the dynamics of the transition amplitude is determined up to two critical times $t_{1}$ and $t_{2}$. As we can see from Eq. (33), the long-time limit is determined by $\lambda^{2}t^{4}\gg 2\lambda t$, which yields $t\gg(2/\lambda)^{1/3}$; whereas the short-time limit is determined by $\lambda^{2}t^{4}\ll 2\lambda t$, which yields $t\ll(2/\lambda)^{1/3}$. Hence, we may approximate the critical times $t_{1}$ and $t_{2}$ by $-(2/\lambda)^{1/3}$ and $(2/\lambda)^{1/3}$ respectively. In Fig. 4a, we compare the analytical approximation $|U_{1}|^{2}=|U_{1}^{(I)}|^{2}\cup|U_{1}^{(II)}|^{2}\cup|U_{1}^{(III)}|^{2}$ of the transition probability with the numerical solution of Eq. (33). The result shows that the transition probability is well-described by $|U_{1}^{(I)}|^{2}\cup|U_{1}^{(II)}|^{2}$ for $t\in(-\infty,t_{2})$. But for $t\in[t_{2},\infty)$, the transition probability, which approaches a non-zero limiting value in the long-time limit, may not be well-described by $|U_{1}^{(III)}|^{2}$, as the Bessel function approaches zero at $t\rightarrow\infty$, e.g., $w_{1}(t)=\sqrt{t}J_{1/6}(\frac{\beta}{12}t^{3})\approx\sqrt{\frac{6}{\pi\lambda}}t^{-1}\cos(\frac{\beta}{12}t^{3}-\frac{\pi}{3})$. In order to fix this problem, we modify the transition amplitude $U_{1}^{(III)}$ by replacing $c_{1}^{(III)}w_{1}+c_{2}^{(III)}w_{2}$ by $\rho e^{iS}+c_{1}^{(III)}w_{1}+c_{2}^{(III)}w_{2}$, where $S\approx\lambda t^{3}/3+S_{0}$ and $\rho$ is determined by Eq. (35). A direct computation yields

[TABLE]

where $\dot{S}^{2}=|f|^{2}+\lambda^{2}t^{4}$ and $\tilde{U}_{1}^{(II)}\equiv U_{1}^{(II)}-\rho_{f}F\cos(\lambda t^{3}/3+S_{0})$. For $t\rightarrow\infty$, as $w_{1}$ and $w_{2}$ decreases as $t^{-1}$, the final transition probability $|U_{1}(\infty)|^{2}$ is given by $\rho_{f}^{2}$, which is determined by Eq. (38a). In Fig. 4b, we compare the modified analytical approximation of the transition probability with the numerical solution of Eq. (33). The result shows that the transition probability $|U_{1}|^{2}$ is well-described by $|U_{1}^{(I)}|^{2}\cup|U_{1}^{(II)}|^{2}\cup|U_{1}^{(III)}|^{2}$ in the whole time range after the replacement $c_{1}^{(III)}w_{1}+c_{2}^{(III)}w_{2}\rightarrow\rho e^{iS}+c_{1}^{(III)}w_{1}+c_{2}^{(III)}w_{2}$.

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 analytical formulas in the whole-time range by comparison with exact results.

In future works, we would like to extend our study to super-linear 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.

Acknowledgements.

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\rangle$ and the first excited state $|e\rangle$, then the Hamiltonian of the atom may be written as $\hat{H}_{atom}=E_{e}|e\rangle\langle e|+E_{g}|g\rangle\langle g|$. For the case when the atom interacts with an external electric field under dipole approximation, the interaction Hamiltonian becomes $\hat{H}_{int}=-\hat{\mathbf{d}}\cdot\mathbf{E}(\hat{\mathbf{r}}_{cm})$, where $\hat{\mathbf{r}}_{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 garrison2008quantum .

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 $|\psi(t)\rangle=c_{e}(t)|e\rangle+c_{g}(t)|g\rangle$, the Schrödinger equations for the two wave amplitudes have the form

[TABLE]

where $\hat{\mathbf{d}}\equiv q\hat{\mathbf{x}}$ 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, $\langle e|\hat{\mathbf{d}}|e\rangle=\langle g|\hat{\mathbf{d}}|g\rangle=0$. Let us denote the energies of the two atomic levels as $E_{e}=\frac{\hbar}{2}\omega_{0}$ and $E_{g}=-\frac{\hbar}{2}\omega_{0}$, where $\hbar\omega_{0}\equiv E_{e}-E_{g}$ is the energy difference between the two atomic levels. For a monochromatic wave, the electric field may be written as $\mathbf{E}=\frac{1}{2}\{E_{0}\exp[i\int_{0}^{t}\omega(\tau)d\tau]\mathbf{e}_{p}+h.c.\}$, where $E_{0}$ is the complex electric field amplitude at the position $\mathbf{r}_{cm}$, and $\mathbf{e}_{p}$ is the polarization vector of the incident electric field. If we denote the transition dipole moment $\langle e|\hat{\mathbf{d}}|g\rangle$ as $\mathbf{d}_{ge}$, the Schrödinger equations for the two wave amplitudes may be simplified as

[TABLE]

In the reference frame rotating about the $z$-axis with frequency $\omega(t)$, the state of the two-level atom is $|\psi^{\prime}(t)\rangle=\exp(\frac{i}{\hbar}\int_{0}^{\tau}\omega(\tau)dtS_{z})|\psi(t)\rangle$. Then the relationships between the wave amplitudes in the rotated and unrotated frame are

[TABLE]

Substitution of Eq. (41) into Eqs. (40a) - (40b) yields

[TABLE]

where $\Delta\equiv\omega-\omega_{0}$ is the laser detuning. Applying the rotating wave approximation, we neglect the fast-oscillating terms like $\exp\{\pm 2i\int_{0}^{t}\omega(\tau)d\tau\}$ in the Schrödinger equations, and obtain

[TABLE]

where $f\equiv-\frac{E_{0}^{*}}{2\hbar}\mathbf{d}_{ge}\cdot\mathbf{e}_{p}^{*}$ is the Rabi frequency for the transition dipole moment. For the special case that $\omega=\omega_{0}$, the laser detuning $\Delta$ vanishes, and the coupled-mode equations become $i\dot{a}_{e}=fa_{g}$ and $i\dot{a}_{g}=f^{*}a_{e}$, which are equivalent to

[TABLE]

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. (43) is solved by $a_{e}=\sin ft$ and $a_{g}=\cos ft$, or equivalently $c_{e}=e^{-\frac{i}{2}\int_{0}^{t}\omega(\tau)d\tau}\sin ft$ and $c_{g}=e^{\frac{i}{2}\int_{0}^{t}\omega(\tau)d\tau}\cos ft$ in the lab frame, which results in the occupation probabilities $|c_{e}|^{2}=\sin^{2}ft$ and $|c_{g}|^{2}=\cos^{2}ft$.

Appendix B The Landau-Zener Model

In this appendix, we provide supplementary information on the large variable expansion of the transition probability in the Landau-Zener model. In the Landau-Zener model, the wave amplitude after the change of variable $U_{1}=C_{1}\exp\{\frac{i\alpha}{2}\int_{0}^{t}sds\}$ is governed by the Schrödinger equation $\ddot{U}_{1}+J(t)U_{1}=0$, where $J(t)=|f|^{2}-\frac{i\alpha}{2}+\frac{\alpha^{2}t^{2}}{4}$. After the change of variable $z=e^{-i\pi/4}\alpha^{1/2}t$ and $n=i|f|^{2}/\alpha$, the Schrödinger equation becomes the parabolic cylinder equation $U_{1}^{\prime\prime}+(n+\frac{1}{2}-\frac{1}{4}z^{2})U_{1}=0$, and is solved by the parabolic cylinder function $D_{n}(z)$, which has the following asymptotic expansion for large $|z|$

[TABLE]

For $\alpha>0$, we have $z=e^{-i\pi/4}\sqrt{|\alpha|}t$, whereas for $\alpha<0$, we have $z=e^{i\pi/4}\sqrt{|\alpha|}t$. For $t\rightarrow-\infty$, we have $z=e^{\pm i3\pi/4}R$ with $R\equiv\sqrt{|\alpha|}|t|$, which yields

[TABLE]

Hence, $|C_{1}(-\infty)|=|A_{\pm}|e^{\mp 3\pi\gamma/2}$ and $|C_{2}(-\infty)|=0$, which violates the initial conditions $C_{1}(-\infty)=0$ and $|C_{2}(-\infty)|=1$. We now check the linearly independent solutions $D_{-n-1}(\mp iz)=D_{-n-1}(e^{\pm i\pi/4}R)$. From Eq. (44a), we obtain

[TABLE]

which yields $|C_{1}(-\infty)|=0$ and $|C_{2}(-\infty)|=|A_{\pm}||\gamma|^{-1/2}e^{|\gamma|\pi/4}$. Hence, the choice of parameters $A_{\pm}=|\gamma|^{1/2}e^{-|\gamma|\pi/4}$ satisfies the initial conditions $C_{1}(-\infty)=0$ and $|C_{2}(-\infty)|=1$. For $t\rightarrow+\infty$, we have $\mp iz=e^{\mp 3i\pi/4}R$, which yields

[TABLE]

A direct computation yields $|C_{2}(\infty)|^{2}=e^{-2\pi|\gamma|}$ and $|C_{1}(\infty)|^{2}=1-e^{-2\pi|\gamma|}$.

Appendix C Large variable solution of $y^{\prime\prime}\pm i\Delta y^{\prime}+|f|^{2}y=0$

In this appendix, we consider a slightly general form of Eqs. (3a) and (3b), and analyze the large variable solution of $y^{\prime\prime}\pm i\Delta y^{\prime}+|f|^{2}y=0$, where $\Delta\equiv\alpha t+\frac{1}{2}\beta t^{2}+\cdots$ is a polynomial of time. Let us denote $y=UV$, where $V$ obeys $\pm i\Delta\dot{V}+|f|^{2}V=0$, and is solved by $V=V_{0}\exp\{\mp i|f|^{2}\int_{t_{0}}^{t}\Delta^{-1}dt^{\prime}\}$. Substitution of $y=UV$ into $y^{\prime\prime}\pm i\Delta y^{\prime}+|f|^{2}y=0$ yields

[TABLE]

For large $|t|$, Eq. (48) is approximated by $\ddot{U}\pm i\nu\dot{U}=0$, which is solved by $U\approx U_{1}\int_{t_{0}}^{t}\exp\{\mp i\int_{t_{0}}^{t^{\prime}}\Delta dt^{\prime\prime}\}dt^{\prime}+U_{0}$. For the Landau-Zener case that $\Delta=\alpha t$, we obtain $V=V_{0}\left(t/t_{0}\right)^{\mp i\gamma}$, $U\approx U_{1}e^{\pm i\alpha t_{0}^{2}/2}\sqrt{\frac{2}{\pm i\alpha}}\int_{z_{0}}^{z}e^{-z^{\prime 2}}dz^{\prime}+U_{0}$, where $\gamma\equiv|f|^{2}/\alpha$ and $z\equiv\sqrt{\pm\frac{i\alpha}{2}}t$. Substitution of the result into $y=UV$ yields

[TABLE]

where the complementary error function $\operatorname{erfc}(z)$ has the following asymptotic expansion

[TABLE]

Hence, for large $|t|$, we obtain

[TABLE]

As a result, the wave amplitudes $C_{1}$ and $C_{2}$ in Eqs. (3a) and (3b) subjected to the conditions $C_{1}(-\infty)=0$ and $|C_{2}(-\infty)|=1$ are asymptotically determined by

[TABLE]

Hence, we obtain $|C_{2}(\infty)|/|C_{2}(-\infty)|=e^{-\pi\gamma}$. We now analyze the parabolic model with $\Delta=\alpha t+\frac{1}{2}\beta t^{2}$. Following similar procedures described previously, we obtain

[TABLE]

where $U_{2}\equiv U_{1}\exp\{\pm i(\frac{\alpha}{2}t_{0}^{2}+\frac{\beta}{6}t_{0}^{3})\}$. We now analyze the large $|t|$ behavior of the integral $\int_{t}^{\infty}e^{-\pi_{n}(t^{\prime})}dt^{\prime}$, where $\pi_{n}(t)$ is an $n$-order polynomial of $t$, e.g., $\pi_{3}(t)=\pm i(\frac{\alpha}{2}t^{2}+\frac{\beta}{6}t^{3})$. A direct computation gives

[TABLE]

which yields

[TABLE]

Hence, we obtain

[TABLE]

where $c_{1}=-V_{0}U_{2}$ and $c_{2}=V_{0}\{U_{0}+U_{2}\int_{t_{0}}^{\infty}\exp[\mp i(\frac{\alpha}{2}t^{2}+\frac{\beta}{6}t^{3})]dt\}$. Hence, we obtain

[TABLE]

where $|C_{2}(\infty)|/|C_{2}(-\infty)|=1$, which is not consistent with the exact solution of the final transition probability. Hence, we conclude that the asymptotic solutions alone cannot guarantee a correct final transition probability, and a more careful treatment of the connection problem is needed for the parabolic model.

Appendix D Integrals involving products of Airy, parabolic cylinder, and Bessel functions

To begin with, let us derive some general results regarding integrals of the form $I_{n}\equiv\int z^{n}y^{2}dz$ and $J_{n}\equiv\int z^{n}y_{1}y_{2}dz$, where $y_{1}$ and $y_{2}$ are solutions of the differential equation $y^{\prime\prime}=fy$. If we write $I_{n}=Py^{2}+Qyy^{\prime}+Ry^{\prime 2}$, we obtain $J_{n}=Py_{1}y_{2}+\frac{1}{2}Q(y_{1}y_{2}^{\prime}+y_{2}y_{1}^{\prime})+Ry_{1}^{\prime}y_{2}^{\prime}$, where $P$ and $Q$ satisfy $P=\frac{1}{2}R^{\prime\prime}-fR$ and $Q=-R^{\prime}$, and $R$ is a solution of the third-order differential equation $R^{\prime\prime\prime}-4fR^{\prime}-2f^{\prime}R=2z^{n}$. As a consequence, $\mathcal{R}\equiv(n+3)(n+2)(n+1)R/2-z^{n+3}$ is a solution of the third-order differential equation $\mathcal{R}^{\prime\prime\prime}-4f\mathcal{R}^{\prime}-2f^{\prime}\mathcal{R}=2z^{n+3}(2(n+3)f/z+f^{\prime})$. Let us denote $y_{1}$ and $y_{2}$ as two independent solutions of $y^{\prime\prime}=fy$ and $\mathscr{W}\{y_{1},y_{2}\}$ being the Wronskian of $y_{1}$ and $y_{2}$. A direct computations shows that $\mathscr{W}$ is a constant. Then for any linear combinations of $y_{1}$ and $y_{2}$, i.e., $y\equiv\alpha y_{1}+\beta y_{2}$, we have

[TABLE]

D.1 Integrals of products of Airy functions

In this sub-appendix, we evaluate indefinite integrals of the form $I_{n}\equiv\int z^{n}y^{2}dz$ and $J_{n}\equiv\int z^{n}y_{1}y_{2}dz$ in terms of Airy functions and their first derivatives, where $y_{1}$ and $y_{2}$ are solutions of the Airy differential equation $y^{\prime\prime}=zy$. A direct computation shows the following recursion relation

[TABLE]

where $R_{n}$ satisfies the recursion relation

[TABLE]

which yields $R_{3}=-(z^{3}+3)/7$ and $R_{4}=-(z^{4}+4z)/9$. Hence, a direct computation yields

[TABLE]

Let us denote $y_{1}=\operatorname{Ai}(z)$, $y_{2}=\operatorname{Bi}(z)$, and $y=\alpha y_{1}+\beta y_{2}$. From Eq. (D), we obtain

[TABLE]

D.2 Integrals of products of parabolic cylinder functions

In this sub-appendix, we evaluate indefinite integrals of the form $I_{n}\equiv\int z^{n}y^{2}dz$ and $J_{n}\equiv\int 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^{\prime\prime}=(a+\frac{1}{4}z^{2})y$. If we write $I_{n}\equiv\left(\frac{1}{2}R_{n}^{\prime\prime}-zR_{n}\right)y^{2}-R_{n}^{\prime}yy^{\prime}+R_{n}y^{\prime 2}$, we obtain $R_{n}^{\prime\prime\prime}-(4a+z^{2})R_{n}^{\prime}-zR=2z^{n}$. We may expand $R_{n}$ in a series as $\sum_{k=0}^{\infty}c_{k}z^{k}$, where the coefficients $c_{k}$ satisfies the following three-term recursion relation

[TABLE]

A direct computation yields

[TABLE]

which corresponds to $R_{1}=-2$ and $R_{3}=(16a-2z^{2})/3$. Let us denote $y_{1}=U(a,z)$, $y_{2}=U(-a,\mp iz)$, and $y=\alpha y_{1}+\beta y_{2}$. From Eq. (D), we obtain

[TABLE]

D.3 Integrals of products of Bessel functions

In this sub-appendix, we evaluate indefinite integrals of the form $I_{n}\equiv\int\tau^{n}y^{2}d\tau$ and $J_{n}\equiv\int\tau^{n}y_{1}y_{2}d\tau$ in terms of Bessel functions and their first derivatives, where $y_{1}$ and $y_{2}$ are both solutions of the differential equation $y^{\prime\prime}=-\lambda^{2}\tau^{4}y$. Using the ansatz $I_{n}=\left(\frac{1}{2}R_{n}^{\prime\prime}+\lambda^{2}\tau^{4}R_{n}\right)y^{2}-R_{n}^{\prime}yy^{\prime}+R_{n}y^{\prime 2}$, we obtain $R_{n}^{\prime\prime\prime}+4\lambda^{2}\tau^{4}R_{n}^{\prime}+8\lambda^{2}\tau^{3}R_{n}=2\tau^{n}$, where $R_{n}$ satisfies the following recursion relation

[TABLE]

A direct computation yields

[TABLE]

which corresponds to $R_{3}=1/(4\lambda^{2})$, $R_{4}=\tau/(6\lambda^{2})$ and $R_{5}=\tau^{2}/(8\lambda^{2})$. We need to evaluate $I_{n}$ with $n=0,1,2$. We may expand $R_{n}$ in a series as $\sum_{k=0}^{\infty}c_{k}\tau^{k}$, where the coefficients $c_{k}$ satisfy the following two-term recursion relation

[TABLE]

A direct computation yields

[TABLE]

where ${}_{p}F_{q}(a_{1},\cdots,a_{p};b_{1},\cdots,b_{q};z)$ is the generalized hypergeometric function of order $p$, $q$. To be specific, let us denote $y_{1}=\sqrt{\tau}J_{1/6}(\lambda\tau^{3}/3)$, $y_{2}=\sqrt{\tau}J_{-1/6}(\lambda\tau^{3}/3)$, and $y=\alpha y_{1}+\beta y_{2}$. From Eq. (D), we obtain ($n=0,1,2$)

[TABLE]

In particular, for $z\rightarrow\infty$ and $p=q-1$, the generalized hypergeometric function ${}_{p}F_{q}(\mathbf{a};\mathbf{b};-z)$ has the following asymptotic expansion in $|\arg(z)|<\pi$

[TABLE]

where $\kappa\equiv q-p+1$, $\nu\equiv\sum_{l}a_{l}-\sum_{l}b_{l}+(q-p)/2$. For $R_{n}$ with $n=0$ or $2$, we have $\kappa=2$, $\nu=-(n+5)/3$, and

[TABLE]

where the plus and minus signs are for $n=0$ and $n=2$ respectively. For $n=0$, we have $z^{-1}\ll z^{-5/6}$; whereas for $n=2$, we have $z^{-7/6}\ll z^{-1}$. Hence, a direct computation yields $R_{2}(\tau)\approx\frac{8}{\beta^{2}\tau}$ and

[TABLE]

For $n=1$, since ${}_{2}F_{3}(1,1;\frac{8}{6},\frac{9}{6},\frac{10}{6};-z)$ is a linear combination of $z^{-1}$ and $z^{-1}\ln z$ near $z=\infty$, $R_{1}(\tau)$ is a linear combination of $\tau^{-2}$ and $\tau^{-2}\ln\tau$. Hence, we have $R_{2}(\tau)\gg R_{1}(\tau)\gg R_{0}(\tau)$ near $\tau=\infty$.