Time-independent approximations for time-dependent optical potentials

TL;DR

This paper proposes an approximate method for solving time-dependent quantum systems with optical potentials by combining invariant-based approaches with perturbation or WKB approximations, demonstrated on solvable models.

Contribution

It introduces a hybrid approach that applies approximate solutions to the invariant eigenvalue problem within the Lewis-Riesenfeld framework for optical potentials.

Findings

The method yields accurate results compared to exact solutions.

Perturbative and WKB approximations are effective in this context.

The approach broadens the applicability of invariant-based techniques.

Abstract

We explore the possibility of modifying the Lewis-Riesenfeld method of invariants developed originally to find exact solutions for time-dependent quantum mechanical systems for the situation in which an exact invariant can be constructed, but the subsequently resulting time-independent eigenvalue system is not solvable exactly. We propose to carry out this step in an approximate fashion, such as employing standard time-independent perturbation theory or the WKB approximation, and subsequently feeding the resulting approximated expressions back into the time-dependent scheme. We illustrate the quality of this approach by contrasting an exactly solvable solution to one obtained with a perturbatively carried out second step for two types of explicitly time-dependent optical potentials.

Equations

$$i\hslash {\partial }_t|{\psi }_n⟩=H(t)|{\psi }_n⟩,$$

$$\frac{dI(t)}{dt}={\partial }_{t}I(t)+\frac{1}{i\hslash }[I(t),H(t)]=0.$$

$$I(t)|{\varphi }_n⟩={\lambda }_n|{\varphi }_n⟩,$$

$$|{\psi }_n⟩=e^{i{\alpha }_n(t)}|{\varphi }_n⟩$$

$$\frac{d\alpha (t)}{dt}=\frac{1}{\hslash }⟨{\varphi }_n|i\hslash {\partial }_t-H(t)|{\varphi }_n⟩.$$

$$I(t)=I_0(t)+ϵI_p(t),$$

$$I(t)|{\varphi }_n⟩={\lambda }_n|{\varphi }_n⟩,\text{and}{I}_0(t)|{\varphi }_n^{(0)}⟩={\lambda }_n^{(0)}|{\varphi }_n^{(0)}⟩.$$

$${\lambda }_n={\lambda }_n^{(0)}+ϵ{\lambda }_n^{(1)}+{ϵ}^2{\lambda }_n^{(2)}+\mathcal{O}({ϵ}^3),\text{and}|{\varphi }_n⟩=|{\varphi }_n^{(0)}⟩+ϵ|{\varphi }_n^{(1)}⟩+{ϵ}^2|{\varphi }_n^{(2)}⟩+\mathcal{O}({ϵ}^3),$$

$${\lambda }_n^{(1)}=⟨{\varphi }_n^{(0)}|I_p|{\varphi }_n^{(0)}⟩,\text{and}|{\varphi }_n^{(1)}⟩={\sum }_{k\ne n}\frac{⟨{\varphi }_k^{(0)}|I_p|{\varphi }_n^{(0)}⟩}{{\lambda }_n^{(0)}-{\lambda }_k^{(0)}}|{\varphi }_k^{(0)}⟩,$$

$$1=⟨{\varphi }_n|{\varphi }_n⟩=⟨{\varphi }_n^{(0)}|{\varphi }_n^{(0)}⟩$$

$$\sum _{k=0}^{\ell }⟨{\varphi }_n^{(\ell -k)}|{\varphi }_n^{(k)}⟩=0.$$

$${|{\psi }_n⟩}^{(1)}=e^{i{\alpha }_n^{(1)}(t)}{|{\varphi }_n⟩}^{(1)},\text{and}{\alpha }^{(1)}(t)=\frac{1}{\hslash }\int dt{}^{(1)}⟨{\varphi }_n|i\hslash {\partial }_t-H(t){|{\varphi }_n⟩}^{(1)}\text{.}$$

$$|{\psi }_n^{\text{WKB}}⟩=e^{i{\alpha }_n^{\text{WKB}}(t)}|{\varphi }_n^{\text{WKB}}⟩,\text{and}{\alpha }^{\text{WKB}}(t)=\frac{1}{\hslash }\int dt⟨{\varphi }_n^{\text{WKB}}|i\hslash {\partial }_t-H(t)|{\varphi }_n^{\text{WKB}}⟩.$$

$${\hat{\varphi }}^{\text{WKB}}(\xi )=\frac{A}{\sqrt{p(\xi )}}{e}^{\frac{i}{\hslash }{\int }^{\xi }p(z)dz}+\frac{B}{\sqrt{p(\xi )}}{e}^{-\frac{i}{\hslash }{\int }^{\xi }p(z)zx}$$

$${\hat{\varphi }}^{\text{WKB}}(\xi )=\frac{C}{\sqrt{q(\xi )}}{e}^{\frac{1}{\hslash }{\int }^{\xi }q(z)dz}+\frac{D}{\sqrt{q(\xi )}}{e}^{-\frac{1}{\hslash }{\int }^{\xi }q(z)dz}$$

$$p(\xi ):=\sqrt{2[\lambda -V(\xi )]}\text{and}q(\xi ):=\sqrt{2[V(\xi )-\lambda ]}.$$

$$H(t)=\frac{p^2}{2m}+\frac{m{\omega }^2}{2}{x}^2+V(x,t)+xE(t).$$

$$I(t)=\frac{1}{2}\left[\alpha (t)p^2+\beta (t)V(x)+\gamma (t)x+\delta (t)\{x,p\}+\epsilon (t)x^2\right].$$

$$\dot{\alpha }$$

$$\beta$$

$$\alpha ={\sigma }^2,\gamma =2m{\sigma }^{2}E(t),\delta =-m\sigma \dot{\sigma },\epsilon =m^2{\dot{\sigma }}^2+m^2\frac{\tau }{{\sigma }^2}.$$

$$\ddot{\sigma }+{\omega }^2\sigma =\frac{\tau }{{\sigma }^3},$$

$$E(t)=\frac{c}{{\sigma }^3}.$$

$${\beta }_p=m{\sigma }^2{\kappa }_p,V_p=c_p{x}^p,{\kappa }_p=\frac{{\tilde{c}}_p}{{\sigma }^{2+p}},$$

$$V_{\text{Gauss}}(x,t)=A(t)\left(e^{-\lambda (t)x^2}-1\right)=\sum _{n=1}^{\infty }{\kappa }_n(t)V_n(x),$$

$$V_{2n}=x^{2n},{\kappa }_n=\frac{(-1{)}^n}{n!}\frac{1}{{\sigma }^{2+2n}},$$

$$V_{\text{sCoulomb}}(x,t)=A(t)\frac{1}{\sqrt{x^2+k^2{a}^2(t)}}=\sum _{n=1}^{\infty }{\kappa }_n(t)V_n(x),$$

$$V_{2n}=x^{2n},{\kappa }_n=\frac{(-1{)}^n(2n)!}{(2n)!!(2n)!!}\frac{1}{{\sigma }^{2+2n}{k}^{1+2n}},$$

$$\sigma =\sqrt{u_1^2+\tau \frac{u_2^2}{W^2}},$$

$$\ddot{u}+{\omega }^{2}u=0,$$

$$\sigma (t)=\frac{1}{\sqrt{2}A\omega }\sqrt{\tau +A^4{\omega }^2+\left(\tau -A^4{\omega }^2\right)\cos (2\omega t)}.$$

Full text

††conference: Time-independent approximations for time-dependent optical potentials

1 Introduction

Since Ashkin’s discovery of the fact that radiation pressure from continuous lasers can be used to trap small particles, almost fifty years ago [1], various types of optical traps have been designed [2]. Different variants of optical traps have found widespead applications to fixate particles [3], atoms [4], molecules [5], and even living cells such as viruses and bacteria [6, 7]. Once the objects are fixed in the potential their properties can be investigated in a very controlled and otherwise inaccessible manner. The general underlying fundamental quantum mechanical description is provided by the time-dependent Schrödinger equation (TDSE) involving explicitly time-dependent Hamiltonians in which the optical potentials typically take on the generic form $V(x,t)=\kappa(t)V(x)$, with $V(x)$ describing the trapping shape and $\kappa(t)$ some time-dependent modulation, e.g. [8]. There exist, however, also optical potentials that can not be factorised in this manner, such as for instance periodic optical lattices [9].

In general, there are only very few exact solutions known to the TDSE and for many physical applications one relies predominatly on approximation methods. Most approximation methods for time-dependent systems are carried out on the level of the time-evolution operator, such as the adiabatic and the sudden approximation [10]. In more concrete settings less general approximation methods have been developed to account for the specifics of the system. For instance, a very active field of research involving time-dependent Hamiltonians is the area of strong laser fields [11]. The systems considered in this context involve Stark potentials of the form $V(x)+xE(t)$, with $E(t)$ describing a laser field and the term $xE(t)$ dominating or being comparable in strength to the potential $V(x)$. For these scenarios approximation schemes have been developed as a mixture of perturbative expansions based on the Du-Hamel formula carried out for time-evolution operator [12, 13, 14]. When iterated, these expressions give rise to various versions of the Born series. In the high intensity regime the two perturbative series, one in $V(x)$ and the other in $xE(t)$, are mixed and terminated after the first iteration. The expressions obtained in this manner are commonly referred to as the strong field approximation [15, 16, 17, 11].

There are, however, methods purposefully designed to solve the TDSE exactly. The Lewis-Riesenfeld method of invariants [18] is one of them. This scheme has been applied successfully to many models and scenarios, such as the harmonic oscillator with time-dependent mass and frequency in one [19] and two dimensions [20, 21], the damped harmonic oscillator [22], a time-dependent Coulomb potential [23], a time-dependent Hamiltonians given in terms of linear combinations of SU(1,1) and SU(2) generators [24, 25, 26], in the inverse construction of time-dependent Hamiltonian [27, 28], for systems on noncommutative spaces in time-dependent backgrounds [29], time-dependent non-Hermitian Hamiltonian systems [30, 31, 32, 33] and other specific systems.

Even though many exact solutions were found, in comparison they are still rare and for more involved potentials the method appears to be nonapplicable. This is partly due to the attempt to complete the Lewis and Riesenfeld method in its entirety. Here we propose to exploit a particular feature of the scheme and modify it when it can only be completed partially. In essence the method consists of three steps, of which often only the first step, that is the construction of the invariants, can be completed. Instead of abandoning the scheme one should recognize that with the completion of the first step one has achieved a major simplification and has transformed the system from a time-dependent first order differential equation to a time-independent eigenvalue equation. Usually one demands in the next step the solvability for this time-independent system, which again is a rare property. For many concrete optical potentials this seems to be entirely unachievable. We explore here the possibility of not terminating the procedure at this stage, but instead to use perturbative theory or the WKB approximation and subsequently feed the approximated expressions back into the Lewis-Riesenfeld approach. Naturally one may also use other approximation methods at this point.

We study here two concrete examples of optical potentials for which an exact method may be obtained and then compare it to the result for which the second step of the Lewis and Riesenfeld method has only been carried out perturbatively. Our manuscript is organised as follows: In section 2 we briefly recall the Lewis-Riesenfeld method of invariants and describe our proposal of altering the second step in the scheme. In section 3 and 4 we investigate in detail two concrete classes of optical potentials. The potentials are chosen in such a way that the full scheme can be completed, hence allowing a comparison with the approximated approach. Our conclusion and an outlook into future problems is provided in section 5.

2 An approximate Lewis-Riesenfeld method of invariants

We start by recalling the key steps of the method of invariants and then describe how they can be modified in an appropriate fashion. The scheme was introduced originally by Lewis and Riesenfeld [18], for the purpose of solving the TDSE

[TABLE]

for the time-dependent or dressed states $\left|\psi_{n}\right\rangle$ associated to the explicitly time-dependent Hamiltonian $H(t)$.

The Lewis-Riesenfeld method of invariants is made up of three main stages: The initial step in this approach consists of constructing a time-dependent invariant $I(t)$ from the evolution equation

[TABLE]

Often this step can be completed and an exact form for the invariant $I(t)$ can be found. In the next step one needs to solve the corresponding eigenvalue system of the invariant $I(t)$

[TABLE]

for time-independent eigenvalues $\lambda_{n}$ and for the time-dependent states $\left|\phi_{n}\right\rangle$. Provided the Hamiltonian $H(t)$ is Hermitian also the invariant $I(t)$ is Hermitian and therefore the eigenvalues $\lambda_{n}$ are guaranteed to be real. The virtue of this equation, compared to the TDSE in (1), is that one has reduced the original evolutionary problem in form of a first order differential equation to an eigenvalue equation in which $t$ simply plays the role of a parameter. Hence one just needs to solve a time-independent eigenvalue problem. To complete this step the system in (3) needs to be solvable. It is this requirement one can weaken and employ *time-independent *approximation methods to complete step two.

The final third step relates the eigenstates in (3) with the complete solution of the TDSE. It was shown in [18] that the states

[TABLE]

satisfy the TDSE (1) provided that the real function $\alpha(t)$ in (4) obeys

[TABLE]

Since all the quantities on the right hand side of (5) have been obtained in the previous steps, one is left with a simple integration in time to determine the phase $\alpha(t)$. These key equations serve mainly for reference purposes and we refer the reader to [18] for more details.

For many systems we might succeed in carrying out the first step in the procedure and construct an explicit expression for the invariant $I(t)$. However, the process stalls often in the second step and for most Hamiltonians the eigenvalue equation for the invariants $I(t)$ in (3) can not be solved exactly. In this case we can, however, use approximation methods. For instance, when the potential can be separated into two terms, with one being dominating the other in absolute value, we can set up a standard time-independent perturbation theory. For this purpose let us briefly recall the main formulae in this approach.

We are splitting the invariant as

[TABLE]

and consider the eigenvalue equation for the full invariant and the unperturbed one sparately

[TABLE]

Assuming that within the perturbation term a small parameter $\epsilon\ll 1$ can be identified, we expand the eigenvalues and the eigenfunctions of the unperturbed invariant as

[TABLE]

with $\lambda_{n}^{(k)}=1/k!\left.d\lambda_{n}/d\epsilon^{k}\right|_{\epsilon=0}$, $\left|\phi_{n}^{(k)}\right\rangle=1/k!\left.d\phi_{n}/d\epsilon^{k}\right|_{\epsilon=0}$. The first order corrections to the eigenvalues and eigenstates of the invariants are then computed in the standard fashion to

[TABLE]

respectively. For orthonormal functions $\phi_{n}$, we obtain further constraints on the normalization of contributions in the series

[TABLE]

Thus if the zero order wavefunction is normalized to $1=\left\langle\phi_{n}^{(0)}\right.\left|\phi_{n}^{(0)}\right\rangle$, we require the higher order wave functions to satisfy the additional constraints

[TABLE]

Next we can use these expressions to obtain an approximate solution to the TDSE. Denoting $\left|\phi_{n}\right\rangle^{(1)}:=\left|\phi_{n}^{(0)}\right\rangle+\epsilon\left|\phi_{n}^{(1)}\right\rangle$ we obtain

[TABLE]

Alternatively, we may also solve the eigenvalue equation (3) by using the WKB approximation $\left|\phi_{n}^{\text{WKB}}\right\rangle$ and compute the phase using that expression

[TABLE]

Assuming that invariant $I(t)$ can be cast into the same form as a time-independent Hamiltonian, with a stndard kinetic energy term and a potential $V(\xi)$, the WKB approximation to first order in $\hbar$, see e.g. [34], for the eigenvalue equation (3) reads

[TABLE]

in the classically allowed region, $\lambda>V(\xi)$ and

[TABLE]

in the classically forbidden region $\lambda<V(\xi)$, where

[TABLE]

The constants $A$, $B$, $C$, $D$ need to be determined by the appropriate asymptotic WKB matching and the normalisation conditions.

Let us now apply this approximation scheme to some concrete systems.

3 Time-dependent potentials with a Stark term

We first demonstrate how to solve the TDSE (1) for the one-dimensional Stark Hamiltonian involving a time-dependent potential $V(x,t)$

[TABLE]

In order to cover optical potentials of the form $V(x,t)$ in our treatment, we are slightly more general than in the standard Stark Hamiltonian where the potential is just depending on $x$ and allow for an explicit time-dependence in the potential $V(x,t)$ as well as in an electric or laser field $E(t)$. At first we assume that the potential factorizes as $V(x,t)=\kappa(t)V(x)$. When the laser field term involving $E(t)$ dominates the potential term and $\kappa(t)=$const several well known and successful approaches have been developed. For instance, the strong field approximation is a mixture of perturbative expansions based on the Du-Hamel formula carried out on the level of the time-evolution operator [15, 16, 17, 11].

In our poposed approach we assume that the first step in the Lewis and Riesenfeld can be carried out and resort to an approximation in form of perturbation thory in the second step.

3.1 Construction of time-independent invariants

In order to carry out the first step in the Lewis-Riesenfeld approach to solve time-dependent systems we need to construct the invariant $I(t)$ by solving equation (2) for a given Hamiltonian, (18) in our case. For this purpose one usually makes an Ansatz by assuming the invariant to be of a similar form as the Hamiltonian

[TABLE]

In our case it involves five unknown time-dependent coefficient functions $\alpha(t)$, $\beta(t)$, $\gamma(t)$, $\delta(t)$ and $\varepsilon(t)$. As ususal we denote the anti-commutator by $\{A,B\}:=AB+BA$. The substitution of (19) into (2) then yields the following first order coupled differential equations as constraints

[TABLE]

Remarkably, despite being overdetermined this system can be solved consistently. We note that the equations in (20) and (21) almost decouple entirely from each other, being only related by $\delta$. We solve (20) first by parameterizing $\alpha(t)=\sigma^{2}(t)$ and integrating twice

[TABLE]

The auxiliary quantity $\sigma$ has to satisfy the nonlinear Ermakov-Pinney (EP) [35, 36] equation

[TABLE]

and in addition the electric field has to be parameterised by the solution of the EP-equation $\sigma$ as

[TABLE]

The constants $c,\tau\in\mathbb{R}$ result from the integrations. We take here $\tau>0$. Using the expression for $\delta$ from (22), we may now also solve the set of equations in (21), obtaining

[TABLE]

with real integration constants $c_{p},\tilde{c}_{p}$ and $p\in\mathbb{R}$. This means that we can not choose the electric field $E(t)$ in our Hamiltonian and the potential $V(x,t)$ entirely a priori and independently from each other. Notice that we may extend the analysis by allowing the constants $c_{p},\tilde{c}_{p}$ to be complex, hence opening up the treatment to include non-Hermitian $\mathcal{PT}$-symmetric Hamiltonians [37, 38, 39].

First we notice that the only time-independent potential is obtained for $p=-2$, so that the potential part in $H(t)$ becomes the solvable Goldman-Krivchenko potential. Crucially, the constraining equations involving the potential (21)  decouple from the remaining ones and since these equations are linear we may solve for potentials that factorize termwise when expanded, that is $V(x,t)=\sum\nolimits_{p}$ $\kappa_{p}(t)V_{p}(x)$. For instance, for a time-dependent Gaussian potential of the form

[TABLE]

we obtain

[TABLE]

where have to restrict $A(t)=\lambda(t)=\sigma^{-2}$. For another widely used potential, the soft Coulomb potential of the form

[TABLE]

with $k$ taken to be a real constant, we obtain

[TABLE]

where we have to restrict $A(t)=1/a(t)=\sigma^{-1}$.

As mentioned, besides the potential, also the electric field is not entirely unconstrained as they are mutually related via the EP-function $\sigma$. However, as we shall demonstrate the solutions of the EP-equation are such that it will still allow for a large class of interesting fields, notably periodic, to be investigated in an exactly solvable manner. It was found by Pinney [36] that the solutions to (23) are

[TABLE]

where $u_{1}$, $u_{2}$ are the two linearly independent solutions of the equation

[TABLE]

and $W=u_{1}\dot{u}_{2}-\dot{u}_{1}u_{2}$ is the corresponding Wronskian. Thus taking the two solutions of (29) to be $u_{1}=A\sin(\omega t)$ and $u_{2}=B\cos(\omega t)$ with $A$, $B\in\mathbb{R}$, the solution to the EP-equation (28) acquires the form

[TABLE]

The function $\sigma(t)$ is regular since $\tau>0$. Therefore the electric field follows to be

[TABLE]

where we have chosen the constants $c=E_{0}$ and $A=\sqrt{\tau}/\omega$ such that $E(0)=E_{0}$. We note that $\omega=\sqrt{\tau}$ is a special point at which $\sigma(t)\rightarrow 1$ and also the field becomes time-independent $E(t)\rightarrow E_{0}$.

Assembling everything we have completed the first step in the Lewis-Riesenfeld construction procedure. The invariant acquires the form

[TABLE]

with $\sigma(t)$ given by (30) and free constants $\tau$, $m$, $\omega$, $c_{p}$, $\tilde{c}_{p}$ and $E_{0}$.

The second step, that is to solve the eigenvalue equation (3), can not be carried out exactly for all invariants $I(t)$ of the form in (32). We therefore resort to a perturbative approach as outlined in the previous section.

3.2 Testing the semi-exact solutions

3.2.1 Exact computation

A good indication about the quality of the perturbation theory and the WKB approximation layed out above can be obtained by comparing both approximations to an exact expression. For most cases this is of course not possible, but taking in (18) the potential for instance to be $V(x,t)=\kappa(t)x^{2}\,$, $\kappa(t)=2c_{\kappa}/\sigma^{4}$, we obtain an exactly solvable system that can serve as a benchmark. In this case the expression (32) for the invariant simply becomes

[TABLE]

with $\alpha$, $\beta$, $\gamma$, $\delta$, $\varepsilon$ as specified in (22). The eigenvalue equation is simplified further when eliminating the anticommutator term $\{x,p\}$ by means of a unitarity tranformation $U=\exp(i\delta x^{2}/2\alpha)$ and the subsequent introduction of the new variable $\xi:=x/\sigma$. We compute

[TABLE]

The eigenvalue equation for the transformed, and in this case time-independent, invariant $\hat{I}\chi(\xi)=\lambda\chi(\xi)$ is then solved by

[TABLE]

where $\mu_{\pm}=\pm(E_{0}^{2}m+4c_{\kappa}\lambda)/\sqrt{m}(2c_{\kappa}+m\tau)^{3/2}-1/2$ and $D_{\nu}(z)$ denotes the parabolic cylinder function. Demanding that the eigenfunctions vanish asymptotically, i.e. $\lim_{\xi\rightarrow\pm\infty}$ $\chi(\xi)=0$, imposes the constraint $\mu_{\pm}=n\in\mathbb{N}_{0}$ and thus quantizes the eigenvalues $\lambda\rightarrow\lambda_{n}$. We discard the solution related to $D_{\mu_{-}}$, as its corresponding eigenvalues are not bounded from below. Hence, we are left with the eigenfuctions and eigenvalues

[TABLE]

The eigenvalues are indeed time-independent as we expect in the context of the Lewis-Riesenfeld approach. Assembling the above and using the orthonormality property of the parabolic cylinder function $\int\nolimits_{-\infty}^{\infty}D_{n}(x)D_{m}(x)dx=n!\sqrt{2\pi}\delta_{nm}$, we obtain the normalized eigenfunction $\phi_{n}=U^{-1}\chi_{n}$

[TABLE]

for the operator $I$ in (33) with

[TABLE]

Finally we compute the phase $\alpha_{n}(t)$ in (4) by means of (5). The right hand side yields111We used here the integrals

$\displaystyle\int\limits_{-\infty}^{\infty}x^{2n}D_{2s+\delta}(x)D_{2r+\bar{\delta}}(x)$ $\displaystyle=$ $\displaystyle(-1)^{s}2^{s+r-n+\frac{1+\delta-\bar{\delta}}{2})}\sqrt{\pi}\Gamma\left(\frac{1}{2}+s+\delta\right)\Gamma\left(2n+1+\delta\right)$

$\displaystyle\times_{3}\tilde{F}_{2}\left(-s,n+1,n+\frac{1}{2}+\delta;\frac{1}{2}+\delta,n-r+1+\frac{\delta-\bar{\delta}}{2};1\right)$

for $n,s,r\in\mathbb{N}_{0}$ and $(\delta,\bar{\delta})=(0,1),(0,0),(1,1)$. The function ${}_{3}\tilde{F}_{2}\left(a,b,c;d,f;z\right)$ is the regularized hypergeometric function defined as

${}_{3}\tilde{F}_{2}\left(a,b,c;d,f;z\right)=\frac{1}{\Gamma\left(d\right)\Gamma\left(f\right)}\sum\limits_{k=0}^{\infty}\frac{(a)_{k}(b)_{k}(c)_{k}}{(d)_{k}(f)_{k}}\frac{z^{k}}{k!},$

(40) with $(a)_{k}=\Gamma\left(a+k\right)/\Gamma\left(a\right)$ denoting the Pochhammer symbol.

[TABLE]

so that phase becomes

[TABLE]

We notice that for $\omega\rightarrow\sqrt{\tau}$ this simply reduces to $\alpha_{n}(t)\rightarrow-\lambda_{n}/m$ and the Hamiltonian becomes time-independent, so that this choice simply describes the time-independent Schrödinger equation.

3.2.2 Perturbative computation

Next we treat the term $V_{p}(x,t)=\kappa(t)x^{2}$ in the Hamiltonian as a perturbation, so that we may view the system as being in the strong field approximation. Accordingly we split up the invariant (33) as $I(t)=I_{0}(t)+\epsilon I_{p}(t)$ with

[TABLE]

and the small expansion parameter is identified as $\epsilon\equiv c_{\kappa}$. First we compute the correction to the eigenvalue of the invariant. Solving the eigenvalue equation (7) and computing the expectation values in (9) we obtain

[TABLE]

As we expect, $\lambda_{n}^{(0)}+$ $c_{\kappa}\lambda_{n}^{(1)}$ is precisely $\lambda_{n}$ in (36) expanded up to first order in $c_{\kappa}$. Next we use (9) to compute the corrections to the wavefunctions. There are only four terms contributing in the infinite sum. We compute

[TABLE]

Finally we evaluate the perturbed expression for the phase $\alpha_{n}^{(1)}(t)$ using equation (5). Up to first order we find

[TABLE]

Notice that for $\omega\rightarrow\sqrt{\tau}$ this simply reduces to $\alpha_{n}^{(1)}(t)=-t(\lambda_{n}^{(0)}+c_{\kappa}\lambda_{n}^{(1)})/m$. We have now obtained the full perturbative solution to the TDSE as $\left|\psi_{n}\right\rangle^{(1)}$ as defined in (13).

3.2.3 WKB computation

We start by determining the classical turning points $\xi_{\pm}$ from the condition $\lambda=V(\xi)$. We find

[TABLE]

so that the WKB quantisation condition

[TABLE]

yields the *exact *time-independent eigenvalues

[TABLE]

as found above in (36). Next we specify WKB wavefunction further. Keeping in the classically forbidden regions $\xi\in(-\infty,\xi_{-})$ and $\xi\in(\xi_{+},\infty)$ only the asymptotically decaying parts in (15) and (16), the corresponding WKB wavefunction are

[TABLE]

and

[TABLE]

respectively. At this point $C_{3}$ is the only undetermined constant left. Carrying out the appropriate WKB matching we obtain for the classically allowed region $\xi\in(\xi_{-},\xi_{+})$ the wavefunction

[TABLE]

We may compute these expressions by using the explicit expressions for the functions $q(\xi)$ and $p(\xi)$. To do so we use the same abbreviated comnstants $a$ and $b$ as in (38) that convert the potential, eigenvalues and turning points into more compact forms

[TABLE]

After a lengthy computation we obtain the WKB wavefunction in the different regions as

[TABLE]

The last remaining constant $C_{3}$ may be fixed by the normalization condition. Converting from the $\hat{I}$ eigenvalue equation back to the $I$ eigenvalue equation with $\phi=U^{-1}\hat{\phi}$ and the variable $\xi$ to $x/\sigma$, the normalisation condition amounts to

[TABLE]

Evaluating the integrals in (55), we find the $n$ independent constant

[TABLE]

Having found the WKB eigenfunction $\phi^{\text{WKB}}$, we can now compute the integrant in (14) that yields the WKB approximated Lewis-Riesenfeld phase

[TABLE]

Thus we have now obtained a WKB approximated solution $\psi_{n}^{\text{WKB}}$ to the time-dependent Schrödinger equation as specified in (14). Let us now compare these three solutions.

3.2.4 WKB versus pertubation theory versus exact solution

In order to obtain an idea about the quality of these approximations let us compute some physical quantities in an exact and perturbative manner and subequently compare them. For the exact case we find the expectation values for the momentum, position and their squares as

[TABLE]

such that the uncertainty relation becomes

[TABLE]

where as usual the squared uncertainty is defined as the squared standard deviation $\Delta A^{2}:=\left\langle\psi_{n}\right|A^{2}\left|\psi_{n}\right\rangle-\left\langle\psi_{n}\right|A\left|\psi_{n}\right\rangle^{2}$ for $A=x,p$. Since the square root is always greater or equal to $1$, the bound in the uncertainty relation $\Delta x\Delta p\geq 1/2$ is always respected.

From the perturbed solution $\left|\psi_{n}\right\rangle^{(1)}$ we find

[TABLE]

and

[TABLE]

These expresssions coincide with the exact expressions expanded up to order one in $c_{k}$. In figure 1 we compare the time-dependent expectation values for $x$, $x^{2}$, $p$, $p^{2}$ computed in an exact way with those computed in a perturbative fashion. In general the agreement is very good for small values of $c_{k}$. Overall the agreement is increasing for large values of $n$ as well as $m$ and for $\omega$ approaching $\sqrt{\tau}$.

A further useful quantity to compute that illustrates the quality of the perturbative approach is the autocorrelation function

[TABLE]

Unlike the expectation values for position, momenta and their squares the autocorrelation function also captures the influence of the time-dependent phase $\alpha(t)$. We depict this function in figure 2. In this case the overall agreement decreases for larger values of $n$.

Next we compare directly the wave functions obtained three alternative ways. Figure 3 shows an extremely good agreement between the WKB approximation and the exact solution, except near the turning points $\xi_{\pm}$ where the WKB approximation is singular. The pertubative solution is in very good agreement with the exact solution for small values of $c_{\kappa}$, as expected. With increasing values of $c_{\kappa}$ the perturbative solution starts to deviate stronger in the negative regime for $\xi$ and large values of $n$.

Let us next see how these properties are inherited in the time-dependent system. Figure 4 displays the real part of the full time-dependent wave function. We observe the oscillation of the turning points with time that enter through the function $\sigma(t)$. As in the time-independent case, extremely good agreement between the WKB approximation and the exact solution, except near the turning points $x_{\pm}$. The perturbative solution slightly overshoots at the maxima and minima, especially in the negative time regime. The discrepancy becomes worse for larger values of $n$, which we do not show here.

4 Goldman-Krivchenko potential with time-dependent perturbation

In trying to identify solvable systems we have seen in equation (25) that the value $p=-2$ is special as in that case the potential becomes the time-independent Goldman-Krivchenko potential [40], being a particular spiked harmonic oscillator [41, 42]. This potential may serve also as a benchmark for which we can solve the eigenvalue equation exactly and campare it to the perturbative solution. Hence we take this potential as our unperturbed system and perturb it by dropping the Stark term and replacing it by $x^{2}E(t)$. Thus we consider the time-dependent Hamiltonian

[TABLE]

It follows from above that the invariant for this system is

[TABLE]

with constraint $E(t)=E_{0}/\sigma^{4}$ and $\sigma$ satisfying the EP-equation (23). Also in this case we may complete the remaining steps in the Lewis-Riesenfeld approach and hence compare the exact and the perturbative solution. We identify $E_{0}\ll 1$ as the expansion parameter in the perturbative series.

4.1 Testing the approximate solution

4.1.1 Exact computation

Using the same similarity and variable transformation as in (34) we obtain the time-independent invariant

[TABLE]

We solve the eigenvalue equation $\hat{I}\Xi(\zeta)=\lambda\Xi(\zeta)$ exactly obtaining the solution

[TABLE]

where $L_{\nu}^{\mu}(z)$ denotes the generalized Laguerre polynomials, $U\left(\nu,\mu,z\right)$ the confluent hypergeometric function and $\nu_{\pm}:=\pm\left(2+b-2\lambda/a\right)/4$, $a=\sqrt{\tau m^{2}+2mE_{0}}$, $b:=\sqrt{1+4m^{2}\Omega^{2}}$. Demanding again that the eigenfunctions vanish asymptotically, i.e. $\lim_{\zeta\rightarrow\pm\infty}$ $\Xi(\zeta)=0$, imposes $\nu_{\pm}=n\in\mathbb{N}_{0}$ and thus quantizes $\lambda$. We discard the solution related to $U$, as its corresponding eigenvalues are not bounded from below, leading to the eigenfuctions and eigenvalues

[TABLE]

Assembling everything we obtain for the operator $I(t)$ the normalized eigenfunction from $U^{-1}\Xi_{n}(x/\sigma)$ as

[TABLE]

when $\mathop{\mathrm{R}e}b>-2$ and $\mathop{\mathrm{R}e}(a/\sigma^{2})>0$. This completes the second step in the Lewis-Riesenfeld approach. In the third and last step we determine the phase $\alpha$ by means of (5). The right hand side is computed once more to

[TABLE]

so that the phase acquires the same form as in the previous example

[TABLE]

Let us now compare these expressions with those obtained in the perturbative computation.

4.1.2 Perturbative computation

We treat the term $V_{p}(x,t)=x^{2}E(t)$ with $E(t)=E_{0}/\sigma^{4}$ and $E_{0}\ll 1$ in the Hamiltonian as a perturbation. Accordingly we split up the invariant (33) as $I(t)=I_{0}(t)+E_{0}I_{p}(t)$ with

[TABLE]

The zeroth order wavefunction $\left|\phi_{n}^{(0)}\right\rangle$ is simply $\phi_{n}(x)$ in (68) with $E_{0}=0$. From (7) and (9) we compute first two terms in the perturbative series for the eigenvalues

[TABLE]

As expected the eigenvalues are time-independent and $\lambda_{n}^{(0)}+E_{0}\lambda_{n}^{(1)}$ corresponds to (67) expanded to first order in $E_{0}$. Next we need to compute the infinite sum in (9) to determine the corrections to the wavefunctions. In this case there are only two terms contributing in the infinite sum. We compute

[TABLE]

In the last step we compute the perturbed expression for the phase $\alpha_{n}^{(1)}(t)$ using equation (5). Once more we find up to first order

[TABLE]

so that we have obtained the full perturbative solution to the TDSE as $\left|\psi_{n}\right\rangle^{(1)}$ as defined in (13).

4.1.3 Exact versus perturbative solutions

As previously, we compute several physical quantities to compare the exact and the perturbative solution. The momentum, position, squared momentum and squared position expectation values are computed to

[TABLE]

In order to achieve convergence we had to impose the additional constraint $b=2\ell+1$ with $\ell\in\mathbb{N}_{0}$. The uncertainty relation becomes

[TABLE]

with the lower bound $\Delta x\Delta p\geq 1/2$ always well respected.

Using the perturbed solutions (74) and (75) we compute

[TABLE]

The approximated uncertainty relation results to

[TABLE]

We are now in the position to compare the exact and the pertubative solution. In figure 5 we compare the time-dependent expectation values for $x^{2}$ and $p^{2}$ computed in an exact way with those computed in a perturbative fashion. As in the previous example, the agreement is very good for small values of the expansion parameter, $E_{0}$ in this case. Overall the agreement is increasing for large values of $n$ as well as $m$ and for $\omega$ approaching $\sqrt{\tau}$.

As in the previous example we also compute the autocorrelation function (62) as it captures well the effect from the time-dependent phase $\alpha(t)$. We depict this function in figure 6. Once more, the overall agreement decreases for larger values of $n$.

5 Conclusions

We have explored the possibilty of a modified approximated Lewis and Riesenfeld method by solving the time-independent eigenvalue equation in the second step by means of standard time-independent perturbation theory. We have tested the quality of this approach for two classes of optical potentials by comparing the exact solutions obtained from the completely exact solution of the Lewis-Riesenfeld approach to the approximated ones, the perturbative approach and the WKB approximation. We computed some standard expectation values and the autocorrelation functions in two alternative ways. For the pertubative approach we found in general good agreement which is naturally improved in quality for smaller values of the expansion parameters. The WKB approximation is not limited to these small parameters and only deviates significantly at the turning points. Our semi-exactly solvable approach significantly widens the scope of the Lewis-Riesenfeld method and allows to tackle more complicated physical situations that are not possible to treat when insisting on full exact solvabilty

While in this paper we were mainly concerned with a proof of concept related to the modification of the Lewis-Riesenfeld method, it would naturally be very interesting now to apply the scheme to study more interesting and challenging physical phenomena. To this end one may carry out further approximations. Here we are still assuming that the first step in the approach, i.e. the computation of the invariant, can be carried out in an exact manner. Naturally one may also weaken this requirement and work with an approximated invariant in the second step. We will present this possibility elsewhere [43].

Acknowledgments: RT is supported by a City, University of London Research Fellowship. We thank Carla Figueira de Morisson Faria for discussions and comments.