Time-independent approximations for time-dependent optical potentials

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.


Introduction
Since Ashkin's discovery of the fact that radiation pressure from continuous lasers can be used to trap small particles, almost 50 years ago [1], various types of optical traps have been designed [2]. Different variants of optical traps have found widespread 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) = κ(t)V (x), with V (x) describing the trapping shape and κ(t) some time-dependent modulation, e.g. [8]. There exist, however, also optical potentials that cannot 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 predominantly 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 a e-mail: a.fring@city.ac.uk b e-mail: rebecca.tenney@city.ac.uk Stark potentials of the form V (x) + x E(t), with E(t) describing a laser field and the term x E(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 x E(t), are mixed and terminated after the first iteration. The expressions obtained in this manner are commonly referred to as the strong field approximation [11,[15][16][17].
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 Davydov-Chaban Hamiltonian in presence of time-dependent potential [24], a Bohr Hamiltonian with a time-dependent potential [25], time-dependent Hamiltonians given in terms of linear combinations of SU(1,1) and SU (2) generators [26][27][28], in the inverse construction of time-dependent Hamiltonian [29,30], for systems on noncommutative spaces in time-dependent backgrounds [31], time-dependent non-Hermitian Hamiltonian systems [32][33][34][35] 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 recognise 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 Sect. 2 we briefly recall the Lewis-Riesenfeld method of invariants and describe our proposal of altering the second step in the scheme. In Sect. 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 Sect. 5.

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 for the time-dependent or dressed states |ψ n associated with 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 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) Since all the quantities on the right hand side of (2.5) have been obtained in the previous steps, one is left with a simple integration in time to determine the phase α(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 (2.3) cannot 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 and consider the eigenvalue equation for the full invariant and the unperturbed one separately I (t) |φ n = λ n |φ n , and Assuming that within the perturbation term a small parameter 1 can be identified, we expand the eigenvalues and the eigenfunctions of the unperturbed invariant as The first-order corrections to the eigenvalues and eigenstates of the invariants are then computed in the standard fashion for nondegenerate systems to respectively. For orthonormal functions φ n , we obtain further constraints on the normalisation of contributions in the series n (2.10) Thus, if the zero-order wavefunction is normalised n , we require the higherorder wavefunctions to satisfy the additional constraints Next, we can use these expressions to obtain an approximate solution to the TDSE. Denoting Alternatively, we may also solve the eigenvalue equation (2.3) by using the WKB approximation φ WKB n and compute the phase using that expression (2.14) Assuming that invariant I (t) can be cast into the same form as a time-independent Hamiltonian, with a standard kinetic energy term and a potential V (ξ ), the WKB approximation to first order inh, see, e.g. [36], for the eigenvalue equation (2.15) in the classically allowed region, λ > V (ξ ) and (2.16) in the classically forbidden region λ < V (ξ ), where 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.

Time-dependent potentials with a Stark term
We first demonstrate how to solve TDSE (2.1) for the one-dimensional Stark Hamiltonian involving a time-dependent potential V (x, t) 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 ( When the laser field term involving E(t) dominates the potential term and κ(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 [11,[15][16][17].
In our proposed 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 theory in the second step.

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 Eq. (2.2) for a given Hamiltonian, (3.1) in our case. For this purpose, one usually makes an Ansatz by assuming the invariant to be of a similar form as the Hamiltonian In our case, it involves five unknown time-dependent coefficient functions α(t), β(t), γ (t), δ(t) and ε(t). As usual, we denote the anticommutator by {A, B} := AB + B A. The substitution of (3.2) into (2.2) then yields the following first-order coupled differential equations as constraintṡ Remarkably, despite being overdetermined this system can be solved consistently. We note that the equations in (3.3) and (3.4) almost decouple entirely from each other, being only related by δ. We solve (3.3) first by parameterizing α(t) = σ 2 (t) and integrating twice The auxiliary quantity σ has to satisfy the nonlinear Ermakov-Pinney (EP) [37,38] equation and in addition the electric field has to be parameterised by the solution of the EP-equation σ as The constants c, τ ∈ R result from the integrations. We take here τ > 0. Using the expression for δ from (3.5), we may now also solve the set of equations in (3.4), obtaining with real integration constants c p ,c p and p ∈ R. This means that we cannot 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 ,c p to be complex, hence opening up the treatment to include non-Hermitian PT -symmetric Hamiltonians [39][40][41].
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 potential (3.4) decouple from the remaining ones and since these equations are linear we may solve for potentials that factorize termwise when expanded, that is V ( For another widely used potential, the soft Coulomb potential of the form with k taken to be a real constant, we obtain where we have to restrict A(t) = 1/a(t) = σ −1 .
As mentioned, besides the potential, also the electric field is not entirely unconstrained as they are mutually related via the EP-function σ . 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 [38] that the solutions to (3.6) are where u 1 , u 2 are the two linearly independent solutions of the equation and W = u 1u2 −u 1 u 2 is the corresponding Wronskian. Thus taking the two solutions of (3.12) to be u 1 = A sin(ωt) and u 2 = B cos(ωt) with A, B ∈ R, the solution to EP-equation (3.11) acquires the form The function σ (t) is regular since τ > 0. Therefore, the electric field follows to be where we have chosen the constants c = E 0 and A = √ τ /ω such that E(0) = E 0 . We note that ω = √ τ is a special point at which σ (t) → 1 and also the field becomes timeindependent E(t) → E 0 .
Assembling everything, we have completed the first step in the Lewis-Riesenfeld construction procedure. The invariant acquires the form with σ (t) given by (3.13) and free constants τ , m, ω, c p ,c p and E 0 . The second step, that is to solve the eigenvalue equation (2.3), cannot be carried out exactly for all invariants I (t) of the form in (3.15). We therefore resort to a perturbative approach as outlined in the previous section.

Exact computation
A good indication about the quality of the perturbation theory and the WKB approximation laid 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 (3.1) the potential for instance to be V (x, t) = κ(t)x 2 , κ(t) = 2c κ /σ 4 , we obtain an exactly solvable system that can serve as a benchmark. In this case, expression (3.15) for the invariant simply becomes with α, β, γ , δ, ε as specified in (3.5). The eigenvalue equation is simplified further when eliminating the anticommutator term {x, p} by means of a unitary transformation U = exp(iδx 2 /2α) and the subsequent introduction of the new variable ξ := x/σ . We computê The eigenvalue equation for the transformed, and in this case time-independent, invariant I χ(ξ) = λχ(ξ ) is then solved by where μ ± = ±(E 2 0 m + 4c κ λ)/ √ m(2c κ + mτ ) 3/2 − 1/2 and D ν (z) denotes the parabolic cylinder function. Demanding that the eigenfunctions vanish asymptotically, i.e. lim ξ →±∞ χ(ξ) = 0, imposes the constraint μ ± = n ∈ N 0 and thus quantizes the eigenvalues λ → λ n .
We discard the solution related to D μ − , as its corresponding eigenvalues are not bounded from below. Hence, we are left with the eigenfunctions and eigenvalues 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 for the operator I in (3.16) with Finally, we compute the phase α n (t) in (2.4) by means of (2.5). The right hand side yields 1 so that phase becomes We notice that for ω → √ τ this simply reduces to α n (t) → −λ n /m and the Hamiltonian becomes time independent, so that this choice simply describes the time-independent Schrödinger equation. 1 We used here the integrals for n, s, r ∈ N 0 and (δ,δ) = 0, 1), (0, 0), (1, 1). The function 3F2 (a, b, c; d, f ; z) is the regularized hypergeometric function defined as with (a) k = (a + k) / (a) denoting the Pochhammer symbol.

Perturbative computation
Next, we treat the term V p (x, t) = κ(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 invariant (3.16) as I (t) = I 0 (t) + I p (t) with 26) and the small expansion parameter is identified as ≡ c κ . First, we compute the correction to the eigenvalue of the invariant. Solving the eigenvalue equation (2.7) and computing the expectation values in (2.9), we obtain As we expect, λ n + c κ λ (1) n is precisely λ n in (3.19) expanded up to first order in c κ . Next, we use (2.9) to compute the corrections to the wavefunctions. There are only four terms contributing in the infinite sum. We compute Finally, we evaluate the perturbed expression for the phase α (1) n (t) using Eq. (2.5). Up to first order, we find Notice that for ω → √ τ this simply reduces to α (1) n + c κ λ (1) n )/m. We have now obtained the full perturbative solution to the TDSE as |ψ n (1) as defined in (2.13).

WKB computation
We start by determining the classical turning points ξ ± from the condition λ = V (ξ ). We find , (3.30) so that the WKB quantisation condition yields the exact time-independent eigenvalues (3.32) as found above in (3.19). Next, we specify WKB wavefunction further. Keeping in the classically forbidden regions ξ ∈ (−∞, ξ − ) and ξ ∈ (ξ + , ∞) only the asymptotically decaying parts in ( 2.15) and (2.16), the corresponding WKB wavefunction arê (3.33) andφ 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 ξ ∈ (ξ − , ξ + ) the wavefunctionφ We may compute these expressions by using the explicit expressions for the functions q(ξ ) and p(ξ ). To do so, we use the same abbreviated constants a and b as in (3.21) that convert the potential, eigenvalues and turning points into more compact forms After a lengthy computation, we obtain the WKB wavefunction in the different regions aŝ (3.37) The last remaining constant C 3 may be fixed by the normalisation condition. Converting from theÎ eigenvalue equation back to the I eigenvalue equation with φ = U −1φ and the variable ξ to x/σ , the normalisation condition amounts to Evaluating the integrals in (3.38), we find the n independent constant (3.39) Having found the WKB eigenfunction φ WKB , we can now compute the integrant in (2.14) that yields the WKB approximated Lewis-Riesenfeld phase Thus, we have now obtained a WKB approximated solution ψ WKB n to the time-dependent Schrödinger equation as specified in (2.14). Let us now compare these three solutions.

WKB versus perturbation 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 subsequently compare them. For the exact case, we find the expectation values for the momentum, position and their squares as such that the uncertainty relation becomes where as usual the squared uncertainty is defined as the squared standard deviation A 2 := ψ n | A 2 |ψ n − ψ n | A |ψ n 2 for A = x, p. Since the square root is always greater or equal to 1, the bound in the uncertainty relation x p ≥ 1/2 is always respected.
From the perturbed solution |ψ n (1) we find (3.43) and x p (1) These expressions coincide with the exact expressions expanded up to order one in c k . In Fig. 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 ω approaching √ τ . Unlike the expectation values for position, momenta and their squares the autocorrelation function also captures the influence of the time-dependent phase α(t). We depict this function in Fig. 2. In this case, the overall agreement decreases for larger values of n. Next, we compare directly the wavefunctions obtained three alternative ways. Figure 3 shows an extremely good agreement between the WKB approximation and the exact solution, except near the turning points ξ ± where the WKB approximation is singular. The perturbative solution is in very good agreement with the exact solution for small values of c κ , as expected.
With increasing values of c κ , the perturbative solution starts to deviate stronger in the negative regime for ξ 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 wavefunction. We observe the oscillation of the turning points with time that enter through the function σ (t). As in the time-independent case, extremely good agreement between the WKB approximation and the exact solution,

Goldman-Krivchenko potential with time-dependent perturbation
In trying to identify solvable systems, we have seen in Eq. (3.8) that the value p = −2 is special as in that case the potential becomes the time-independent Goldman-Krivchenko potential [42], being a particular spiked harmonic oscillator [43,44]. This potential may serve also as a benchmark for which we can solve the eigenvalue equation exactly and compare 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 It follows from above that the invariant for this system is with constraint E(t) = E 0 /σ 4 and σ satisfying EP-equation (3.6). 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 1 as the expansion parameter in the perturbative series.

Exact computation
Using the same similarity and variable transformation as in (3.17), we obtain the timeindependent invariant We solve the eigenvalue equationÎ (ζ ) = λ (ζ ) exactly obtaining the solution where L μ ν (z) denotes the generalised Laguerre polynomials, U (ν, μ, z) the confluent hypergeometric function and ν ± : Demanding again that the eigenfunctions vanish asymptotically, i.e. lim ζ →±∞ (ζ ) = 0, imposes ν ± = n ∈ N 0 and thus quantizes λ. We discard the solution related to U , as its corresponding eigenvalues are not bounded from below, leading to the eigenfunctions and eigenvalues Assembling everything, we obtain for the operator I (t) the normalised eigenfunction from U −1 n (x/σ ) as when Reb > −2 and Re(a/σ 2 ) > 0. This completes the second step in the Lewis-Riesenfeld approach. In the third and last step, we determine the phase α by means of (2.5). The right hand side is computed once more to so that the phase acquires the same form as in the previous example Let us now compare these expressions with those obtained in the perturbative computation.

Perturbative computation
We treat the term V p (x, t) = x 2 E(t) with E(t) = E 0 /σ 4 and E 0 1 in the Hamiltonian as a perturbation. Accordingly, we split up invariant (3.16) as The zeroth-order wavefunction φ (0) n is simply φ n (x) in (4.6) with E 0 = 0. From (2.7) and ( 2.9), we compute first two terms in the perturbative series for the eigenvalues As expected, the eigenvalues are time independent and λ (0) n + E 0 λ (1) n corresponds to (4.5) expanded to first order in E 0 . Next, we need to compute the infinite sum in (2.9 ) to determine the corrections to the wavefunctions. In this case, there are only two terms contributing in the infinite sum. We compute In the last step, we compute the perturbed expression for the phase α (1) n (t) using Eq. (2.5). Once more we find up to first order so that we have obtained the full perturbative solution to the TDSE as |ψ n (1) as defined in (2.13).

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 (4.14) In order to achieve convergence, we had to impose the additional constraint b = 2 + 1 with ∈ N 0 . The uncertainty relation becomes with the lower bound x p ≥ 1/2 always well respected. Using the perturbed solutions (4.12) and (4.13), we compute ψ n | x |ψ n (1) = 0, ψ n | x 2 |ψ n (1) = (2n + + 3/2)(mτ − E 0 ) σ 2 τ m 2 , ψ n | p |ψ n (1) = 0, ψ n | p 2 |ψ n (1) = (4n + 2) + 2n + 3/2 2 + 1 The approximated uncertainty relation results to We are now in the position to compare the exact and the perturbative solution. In Fig. 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 ω approaching √ τ . As in the previous example, we also compute the autocorrelation function (3.45) as it captures well the effect from the time-dependent phase α(t). We depict this function in Fig. 6. Once more, the overall agreement decreases for larger values of n.

Conclusions
We have explored the possibility 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 perturbative 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 solvability. The validity of either approach is governed by the validity of the time-independent perturbation theory and the WKB approximation for which explicit expressions can be found in the standard literature.
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 [45].