Time-reversible and norm-conserving high-order integrators for the nonlinear time-dependent Schrödinger equation: Application to local control theory

Nonlinear time-dependent Schrödinger equations contain, by definition, Hamiltonians that depend on the quantum state. Such state-dependent effective Hamiltonians appear in many areas of physics and chemistry. Examples include various nonlinear Schrödinger equations generated by the Dirac–Frenkel time-dependent variational principle,^1–4 e.g., the equations of the multi-configurational time-dependent Hartree method,^5–7 and some numerical methods such as the short-iterative Lanczos algorithm.^8–10 Probably, the best known nonlinear Schrödinger equations, however, are approximate equations for Bose–Einstein condensates,^11,12 in which the Hamiltonian depends on the probability density of the quantum state. The dynamics of a Bose–Einstein condensate is often modeled by solving the celebrated Gross–Pitaevskii equation with a cubic nonlinearity.^13–17 

To solve this equation, several numerical schemes, such as the explicit split-operator algorithm or the time and spatial finite difference methods, have been employed.^18,19 These methods are of low accuracy (in time and/or space) and do not always preserve the geometric properties of the exact solution.²⁰ For example, the Crank–Nicolson finite difference method is geometric but exhibits only second-order convergence with respect to the spatial discretization. To remedy this, the explicit second-order split-operator algorithm,^21–24 commonly used for the linear time-dependent Schrödinger equation, is a great alternative, as it conserves, in some cases, the geometric properties of the exact solution and has spectral accuracy in space. Unfortunately, this algorithm cannot be used for all types of nonlinear time-dependent Schrödinger equations. Indeed, in the case of the Gross–Pitaevskii equation, the algorithm is symmetric and, therefore, time-reversible only because the ordinary differential equation that must be solved when propagating the molecular state with the potential part of the Hamiltonian leaves a nonlinear term invariant in time.¹⁸ We show here that for nonlinear terms of a more general form, this algorithm becomes implicit. If this implicit nature is not taken into account and the explicit version is used, the algorithm loses its time reversibility and efficiency due to its low accuracy that is only of the first order in the time step.

An example of a situation, where a more general nonlinearity appears, is provided by local control theory (LCT). Introduced by Kosloff et al.,²⁵ LCT is a widely used approach to coherent control. In LCT, the pulse is computed on the fly, based on the instantaneous molecular state, in order to increase (or decrease) the expectation value of a specified operator. LCT has been successfully used to control various processes such as energy and population transfer,^25–28 dissociation and association dynamics,^29–32 direction of rotation in molecular rotors,³³ and electron transfer.³⁴ Controlling quantum systems using LCT changes the nature of the time-dependent Schrödinger equation. Because the time dependence of the pulse is determined exclusively by the molecular state, the time-dependent Schrödinger equation becomes autonomous but nonlinear.

The nonlinear nature of LCT is often not acknowledged, and the standard explicit split-operator algorithm²¹ for linear time-dependent Schrödinger equations is used,^29–32,34 instead of its time-reversible, second-order, but implicit alternative. Most previous studies used LCT for applications that required neither high accuracy nor time reversibility and, therefore, could rely on this approximate explicit integrator, which, as we show below, in the context of LCT, indeed, has only first-order accuracy in the time step and is time-irreversible. Such an algorithm, however, would be very inefficient for highly accurate calculations and could not be used at all if exact time reversibility were important. Because this failure of the explicit splitting algorithm in LCT is generic, while its success in the Gross–Pitaevskii equation is rather an exception, it is desirable to develop efficient high-order geometric integrators suitable for a general nonlinear time-dependent Schrödinger equation.

Recently, we presented high-order time-reversible geometric integrators for the nonadiabatic quantum dynamics driven by the linear time-dependent Schrödinger equation with both separable³⁵ and nonseparable³⁶ Hamiltonians. Here, we extend this work to the general nonlinear Schrödinger equation, in order to address the slow convergence and time irreversibility of the explicit split-operator algorithm.

The remainder of the study is organized as follows: In Sec. II, we define the nonlinear time-dependent Schrödinger equation, discuss its geometric properties, and explain how LCT leads to a nonlinear Schrödinger equation. In Sec. III, after demonstrating the loss of geometric properties by Euler methods, we describe how these geometric properties are recovered and accuracy is increased to an arbitrary even order by symmetrically composing the implicit and explicit Euler methods. Then, we describe a general procedure to perform the implicit propagation and derive explicit expressions for the case of LCT. We also show the derivation of the approximate explicit split-operator algorithm for the nonlinear Schrödinger equation, explain how it loses time reversibility, and briefly describe the dynamic Fourier method. Finally, in Sec. IV, we numerically verify the convergence and geometric properties of the integrators by controlling, using LCT, either the population or energy transfer in a two-state two-dimensional model of retinal.³⁷ 

The nonlinear time-dependent Schrödinger equation is the differential equation

describing the time evolution of the state ψ_t driven by the nonlinear Hamiltonian operator $Ĥ(ψt)$⁠, which depends on the state of the system. This dependence on ψ_t is what distinguishes the equation from the linear Schrödinger equation. As the notation in Eq. (1) suggests, we shall always assume that while the operator $Ĥ:ψ↦Ĥ(ψ)ψ$ is nonlinear, for each ψ, the operator $Ĥ(ψ):ϕ↦Ĥ(ψ)ϕ$ is linear. We will also assume that $Ĥ(ψ)$ has real expectation values $⟨Ĥ(ψ)⟩ϕ≔⟨ϕ|Ĥ(ψ)ϕ⟩$ in any state ϕ, which for a linear operator implies that it is Hermitian, i.e., $Ĥ(ψ)†=Ĥ(ψ)$⁠, or, more precisely, that for every ψ, ϕ, and χ,

A paradigm of a nonlinear Schrödinger equation is the Gross–Pitaevskii equation,^13,14,17 in position representation expressed as

where the real coefficient C is positive for a repulsive interaction and negative for an attractive interaction. This equation has a cubic nonlinearity and is useful, e.g., for approximate modeling of the dynamics of a Bose–Einstein condensate.²⁰ Many other examples are provided by the Dirac–Frenkel variational principle,^1,2 which approximates the exact solution of a linear Schrödinger equation with the Hamiltonian $Ĥ$ by an optimal solution of a predefined, restricted form within a certain subset (called the approximation manifold) of the Hilbert space. This optimal solution ψ_t satisfies the equation

where δψ_t is an arbitrary variation in the approximation manifold. Equation (3) is equivalent to the nonlinear Eq. (1) with an effective state-dependent Hamiltonian

where $P^(ψt)$ is the projection operator on the tangent space to the approximation manifold at the point ψ_t.^4,38 (Note that in the very special case, where the projector does not depend on ψ_t, the resulting Schrödinger equation remains linear. This happens in the Galerkin method, in which ψ_t is expanded in a finite, time-independent basis, and the approximation manifold is a vector space.⁴)

With initial condition $|ψt0〉$ and assuming that t ≥ t₀, Eq. (1) has the formal solution $|ψt〉=Û(t,t0;ψ)|ψt0〉$ with the exact evolution operator given by

where the dependence of $Û$ on ψ was added as an argument to emphasize the nonlinear character of Eq. (1). Expression (4) is obtained by solving the differential equation

with the initial condition $Û(t0,t0;ψ)=1$⁠. The Hermitian adjoint of $Û(t,t0;ψ)$ is the operator

where $T¯$ denotes the reverse time-ordering operator.

The nonlinearity of $Ĥ$ leads to the loss of some geometric properties, even if Eq. (1) is solved exactly. Indeed, since the Hamiltonian is nonlinear, the exact evolution operator is also nonlinear.

An operator $Û$ is said to preserve the inner product (or to be unitary) if $⟨Ûψ|Ûϕ⟩=⟨ψ|ϕ⟩$⁠. The exact evolution operator does not preserve the inner product because

if $ψt0≠ϕt0$⁠, where we used the property (5) of the Hermitian adjoint of $Û$ to obtain the third line. The exact nonlinear evolution operator is, therefore, not symplectic because it does not preserve the associated symplectic two-form⁴ 

The nonlinear evolution does not conserve energy

since

where the first and third terms in the intermediate step cancel each other because ψ_t satisfies the nonlinear Schrödinger equation (1). Note, however, that in special cases, such as the Gross–Pitaevskii equation, there exist modified energy functionals that are conserved.²⁰ 

An operator $Û$ is said to conserve the norm ‖ψ‖ ≔ ⟨ψ|ψ⟩^1/2 if $‖Ûψ‖=‖ψ‖$⁠. The exact nonlinear evolution operator $Û≡Û(t,t0;ψ)$ conserves the norm because

where we used relation (5).

In the theory of dynamical systems, an adjoint $Û(t,t0;ψ)*$ of $Û(t,t0;ψ)$ is defined as the inverse of the evolution operator taken [in contrast to the Hermitian adjoint (5)] with a reversed time flow,

An operator $Û(t,t0;ψ)$ is called symmetric if it is equal to its adjoint, i.e., if $Û(t,t0;ψ)=Û(t,t0;ψ)*$⁠. A symmetric operator is also time-reversible because propagating a molecular state $ψt0$ forward to time t and then backward to time t₀ yields $ψt0$ again, i.e.,

For the operator (4), the reverse evolution operator is given by the anti-time-ordered exponential

and therefore, the adjoint is

This shows that the exact evolution operator (4) is symmetric and time-reversible.

Finally, a time evolution is said to be stable if for every ϵ > 0, there is a δ > 0 such that the distance between two states satisfies the condition

Although the exact evolution conserves the norm, because it does not conserve the inner product, we cannot, in general, say anything about preserving the distance,

Moreover, since the sign of the real part of the difference of the inner products can be arbitrary, we cannot deduce anything about stability.

We now show that the local coherent control of the time evolution of a quantum system with an electric field provides another example of a nonlinear Schrödinger equation. The quantum state |ψ_t⟩ of a system interacting with an external time-dependent electric field $E→(t)$ evolves according to the linear time-dependent Schrödinger equation

with a time-dependent Hamiltonian

equal to the sum of the Hamiltonian $Ĥ0$ of the system in the absence of the field and the interaction potential $V^int(t)$⁠. Within the electric-dipole approximation,³⁹ the interaction potential is

where the vector function $μ→(q^)$ of the position operator $q^$ denotes the electric-dipole operator of the system. Direct integration of Eq. (16) with the initial condition $|ψt0〉$ leads to the formal solution $|ψt〉=Û(t,t0)|ψt0〉$ with the exact evolution operator given by the time-ordered exponential

This exact evolution operator has many important geometric properties: it is linear, unitary, symmetric, time-reversible, and stable.^4,40–42 Because it is unitary, the evolution operator conserves the norm as well as the inner product and symplectic structure.⁴ However, since the Hamiltonian is time-dependent, the Schrödinger equation (16) is a nonautonomous differential equation⁴² and, as a consequence, does not conserve energy. For a more detailed presentation and discussion of the above properties, we refer the reader to Ref. 36.

Contrary to Eq. (18), the electric field used in LCT, called control field and denoted by $E→LCT(t)$⁠, is not known explicitly as a function of time. Instead, it is chosen “on the fly” according to the current state |ψ_t⟩ of the system, in order to increase or decrease the expectation value $⟨Ô⟩ψt≔⟨ψt|Ô|ψt⟩$ of a particular operator $Ô$ in the state |ψ_t⟩. More precisely, the control field is computed so that the time derivative of the expectation value,

remains positive or negative at all times. If the operator $Ô$ commutes with the system’s Hamiltonian $Ĥ0$⁠, this goal is achieved by using the field

where λ > 0 is a parameter that scales the intensity of the control field, and the sign in Eq. (21) is chosen according to whether one wants to increase or decrease $⟨Ô⟩ψt$⁠. This claim is proven by inserting the definition (21) of $E→LCT(t)$ into Eq. (20), which yields

for the derivative of the expectation value. This equation confirms that a strictly increasing or strictly decreasing evolution of $⟨Ô⟩ψt$ is guaranteed only if $[Ĥ0,Ô]=0$⁠, largely reducing the choice of operators $Ô$⁠, whose expectation values can be controlled monotonically.

The left-hand side of Eq. (21) suggests that the control field can be viewed either as a function of time or a function of the molecular state [i.e., $E→LCT(t)≡E→LCT(ψt)$]. More precisely, the control field does not depend on time explicitly but only implicitly through the dependence on ψ_t. Therefore, the time-dependent Schrödinger equation changes from a nonautonomous linear to an autonomous nonlinear differential equation.⁴² By acknowledging this nonlinear character, the interaction potential from Eq. (18) becomes

and Eq. (16) becomes an example of a nonlinear time-dependent Schrödinger equation (1) with the Hamiltonian operator $Ĥ(ψ)≔Ĥ0+V^LCT(ψ)$⁠.

Numerical propagation methods for solving the nonlinear equation (1) obtain the state at time t + Δt from the state at time t by using the relation

where Δt denotes the numerical time step and $Ûappr(t+Δt,t;ψ)$ is an approximate nonlinear evolution operator depending on ψ. By construction, all propagation methods converge to the exact solution in the limit Δt → 0. As the exact operator, these approximate evolution operators $Ûappr(t+Δt,t;ψ)$ are, therefore, nonlinear and conserve neither the inner product nor the symplectic form. Moreover, nothing can be said about their stability in general. However, some integrators may lose even the remaining geometric properties of the exact evolution: norm conservation, symmetry, and time reversibility. In this section, we simply state the properties that are lost by different methods; detailed proofs are provided in the  Appendix.

The simplest methods, applicable to both separable and nonseparable and both linear and nonlinear Hamiltonian operators, are the explicit and implicit Euler^4,41 methods, which approximate the exact evolution operator, respectively, as

Both methods are only first-order in the time step and, therefore, very inefficient. Moreover, both Euler methods lose the norm conservation, symmetry, and time reversibility of the exact evolution operator.

Composing the implicit and explicit Euler methods yields, depending on the order of composition, either the implicit midpoint method,

or the trapezoidal rule (or Crank–Nicolson method),

Both methods are second-order, symmetric, and time-reversible regardless of the size of the time step.^36,42 Although the trapezoidal rule conserves the norm of a state evolved with a time-independent linear Hamiltonian,³⁶ it loses this property when the Hamiltonian is time-dependent or nonlinear (which results, as we have seen, in an implicit time dependence). In contrast, the implicit midpoint method remains norm-conserving in all cases.

Because they are symmetric, both methods can be further composed using symmetric composition schemes^4,36,42–46 in order to obtain integrators of an arbitrary even order of accuracy in the time step. Indeed, every symmetric method $Ûp$ of an even order p generates a method $Ûp+2$ of order p + 2 if it is symmetrically composed as

where M is the number of composition steps and γ₁, …, γ_M are real composition coefficients, which satisfy the relations $∑n=1Mγn=1$⁠, γ_M+1−n = γ_n, and $∑n=1Mγnp+1=0,$ guaranteeing the increase in the order of accuracy.⁴² 

The simplest composition methods that are used here are the triple-jump⁴³ and Suzuki’s fractal⁴⁴ with M = 3 and M = 5, respectively. Both are recursive and able to generate integrators of arbitrary even orders of accuracy. Sixth-, eighth-, and tenth-order integrators can also be obtained with nonrecursive composition methods,^45,46 which require fewer composition steps and also minimize their magnitude. These composition methods are, therefore, more efficient than the triple-jump and Suzuki’s fractal and will be referred to as “optimal” compositions. For more details on symmetric compositions, see Ref. 36, where they were applied to a linear Schrödinger equation.

The implicit Euler method requires implicit propagation because its integrator [see Eq. (26)] depends on the result of the propagation, i.e., ψ_t+Δt. In the implicit Euler method, ψ_t+Δt is obtained by solving the nonlinear system

which can be written as f(ψ_t+Δt) = 0 with the nonlinear functional

A nonlinear system f(ψ) = 0 can be solved with the iterative Newton–Raphson method, which computes, until convergence is obtained, the solution ψ^(k+1) at iteration k + 1 from ψ^(k) using the relation

where $Ĵ≔δδψf(ψ)$ is the Jacobian of the nonlinear functional f(ψ).

If the initial guess ψ⁽⁰⁾ is close enough to the exact solution of the implicit propagation, the Newton–Raphson iteration (31) is a contraction mapping and by the fixed-point theorem is guaranteed to converge. We use as the initial guess the result of propagating ψ_t with the explicit Euler method [Eq. (25)]. Note that this initial guess is sufficiently close to the implicit solution only if the time step is small. If the time step is too large, the difference between the explicit and implicit propagations becomes too large for the algorithm to converge, and no solution can be obtained.

Equation (31) requires computing the inverse of the Jacobian, which is an expensive task. It is preferable to avoid this inversion by computing each iteration as

where δψ^(k) solves the linear system

We solve this linear system by the generalized minimal residual method,^47–49 an iterative method based on the Arnoldi process^50,51 (see Sec. I of the supplementary material for a detailed presentation of this algorithm).

The procedure presented for solving the implicit propagation is applicable to any nonlinear system whose Jacobian is known analytically. Therefore, the integrators proposed in Secs. III A and III B can be employed for solving any nonlinear time-dependent Schrödinger equation of the form of Eq. (1), i.e., with the Hamiltonian $Ĥ(ψt)$ depending on the state of the system.

To sum up, each implicit propagation step, given by the evolution operator (26), is performed as follows:

1. Compute the initial guess ψ⁽⁰⁾ using the explicit Euler method [see Eq. (25)]. Choose an error threshold ε and a maximum iteration number m.

2. For k = 0, 1, …, m − 1, Do:

3. Compute δψ^(k) by solving the linear system shown in Eq. (33).

4. Compute a new approximate solution ψ^(k+1) using Eq. (32).

5. If ‖f(ψ^(k+1))‖ ≤ ε, take ψ^(k+1) as the solution of the implicit step.

6. End Do. The algorithm fails when k = m and no approximate solution has been found.

In the case of LCT,

and the Jacobian of the nonlinear functional (30) is

To obtain the third row of Eq. (34), we used $δδψ[V^LCT(ψ)]ψ=V^LCT(ψ)$⁠, where the generalized complex derivative⁵² of the interaction potential is given by the bra vector

The algorithms that we described above apply to Hamiltonians that are not only nonlinear but also nonseparable, i.e., to Hamiltonians $Ĥ$ that cannot be written as a sum $Ĥ=T(p^)+V(q^)$ of a momentum-dependent kinetic term and a position-dependent potential term. If the time-dependent Schrödinger equation is linear and its Hamiltonian is separable, the midpoint method remains implicit, but the split-operator algorithms and their compositions yield explicit high-order integrators satisfying most geometric properties (except for the conservation of energy). In the case of LCT, if $Ĥ0$ is separable, so is the total Hamiltonian, which can be written as $Ĥ(ψ)=T^+V^tot(ψ)$⁠, where $V^tot(ψ)≔V^0+V^LCT(ψ)$ is the sum of the system’s and interaction potential energy operators. It is, therefore, tempting to use the split-operator algorithm, with the hope of obtaining an efficient explicit integrator.

More generally, let us assume that the Hamiltonian operator in the general nonlinear Schrödinger equation (1) can be separated as

The approximate evolution operator is given by

in the TV split-operator algorithm and by

in the VT split-operator algorithm. These integrators are norm-conserving but only first-order and time-irreversible. From their definitions (36) and (37), it follows immediately that the TV and VT algorithms are adjoints³⁵ of each other and require, respectively, explicit and implicit propagations. In analogy to the implicit midpoint algorithm from Sec. III B, a second-order method is obtained by composing the two adjoint methods to obtain the TVT split-operator algorithm,

or the VTV split-operator algorithm if the order of composition is reversed. Both TVT and VTV algorithms are norm-conserving, symmetric, and time-reversible. However, these geometric properties are only acquired if the implicit part, i.e., the propagation with the VT algorithm (38) is performed exactly. This requires solving a nonlinear system, which can be performed using the Newton–Raphson method, as described in Sec. III C. This, however, implies abandoning the explicit nature of the split-operator algorithm, which is one of its main advantages over implicit methods for solving linear Schrödinger equations.

The nonlinearity of Eq. (1) is often not acknowledged. As a consequence, the implicit character of Eqs. (37) and (38) is not taken into account, and explicit alternatives of these algorithms are used. For example, instead of using ψ_t+Δt in the VT algorithm (i.e., performing the implicit propagation exactly), the state $ψt,T^Δt≔e−iT^Δt/ℏψt$ obtained after the kinetic propagation is often used to perform the potential propagation. After composition with the TV algorithm, it yields the approximate explicit TVT split-operator algorithm,

This approximate explicit integrator can be used to successfully perform LCT, but, because it depends on $ψt,T^Δt/2$ instead of ψ_t+Δt/2, it is not time-reversible and achieves only first-order accuracy. Like the approximate explicit TVT algorithm, any other explicit algorithm cannot perform the implicit step of the $ÛVT$ algorithm exactly and, therefore, cannot be time-reversible. Yet, the approximate explicit algorithm does conserve the norm (see the  Appendix for proofs of the geometric properties).

To propagate the molecular state using the algorithms presented in Sec. III, an efficient way of evaluating the action of an operator on the molecular state ψ_t is needed. For this, we use the dynamic Fourier method,^21–24 which is applicable to operators that are sums of products of operators of the form $f(x^)$⁠, where $x^$ denotes either the momentum operator $p^$ or the position operator $q^$⁠. The computation of $f(x^)ψt$ is then performed easily, by pointwise multiplication, in the x-representation, in which $x^$ is diagonal. Whenever required, the representation of ψ_t is changed by the Fourier transform. In the numerical examples below, we used the Fastest Fourier Transform in the West 3 (FFTW3) library⁵³ to perform all of the Fourier transforms. Unlike finite difference methods, the Fourier method shows exponential convergence with respect to the number of grid points (see Figs. S1 and S2 of the supplementary material).

We tested the general integrators for the nonlinear Schrödinger equation, presented in Sec. III, by using them for the local control of a two-dimensional two-state diabatic model of retinal taken from Ref. 37. The model describes the cis–trans photo-induced isomerization of retinal—an ultrafast reaction mediated by a conical intersection and the first event occurring in the biological process of vision. The two vibrational modes of the model are the reaction coordinate θ, an angle describing the torsional motion of the retinal molecule, and a vibronically active coupling mode q_c. In the diabatic representation, the Hamiltonian of the system in the absence of the field,

is separable into a sum of the kinetic energy operator

and the potential energy operator with components

Here (all parameters are in eV units), ω = 0.19 is the vibrational frequency of the coupling mode, m⁻¹ = 4.84 · 10⁻⁴ is the inverse mass of the reaction coordinate, W₁ = 3.6 and W₂ = 1.09 determine the depth of the well in the reaction coordinate for the ground and excited electronic states, respectively, χ₂ = 0.1 is the gradient of the linear perturbation in the excited electronic state, E₂ = 2.48 determines the maximum of the excited electronic state in the reaction coordinate, and ξ = 0.19 is the gradient of the linear coupling between the two electronic states. The two diabatic potential energy surfaces (42) and (43) are displayed in Fig. S3 of the supplementary material.

Note that, to distinguish nuclear, electronic, and molecular operators, in this section, the bold face denotes electronic operators expressed as S × S matrices in the basis of S electronic states and the hat $ˆ$ denotes nuclear operators acting on the Hilbert space of nuclear wavefunctions, i.e., square-integrable functions of D continuous dimensions.

In the simulations, the reaction and coupling coordinates are represented on regular grids consisting of 128 points between θ = −π/2 a.u. and θ = π/2 a.u. and 64 points between q_c = −9 a.u. and q_c = 9 a.u. Figure S1 of the supplementary material confirms that this grid is sufficient by showing that the grid representation of the wavepacket is converged even at the final time t_f. We assume the coordinate independence of the electric-dipole moment operator (Condon approximation) and, therefore, can write $μ→^=μ→1^=ϵ→μ1^$⁠, where $ϵ→$ is a constant unit vector in the direction of $μ→$⁠. In this case, the LCT electric field is aligned with $μ→$⁠, and we can write it as $E→LCT=ϵ→ELCT$⁠. As a consequence, we can drop the vector symbol $⃗$ from $μ→^$ and $E→LCT$ in Eqs. (21) and (23) and consider only the analogous scalar equations satisfied by $μ^$ and E_LCT. In addition, we assume the electric-dipole moment operator to have unit transition elements (⁠$μ^12=μ^21=μ^=1$ a.u.) and zero diagonal elements (⁠$μ^11=μ^22=0$⁠). The calculations presented below aim at simulating the photo-excitation step of the photo-isomerization of the retinal molecule. We, therefore, use as the initial state ψ₀ the ground vibrational state of the harmonic fit of the ground potential energy surface [i.e., a two-dimensional Gaussian wavepacket with q₀ = p₀ = (0, 0) and σ₀ = (0.128, 1) a.u.] with initial populations P₁(0) = 0.999 and P₂(0) = 0.001 of the ground and excited electronic states, respectively. The tiny initial seed population of the excited state is essential for the control because it ensures that Eq. (21) does not stay zero at all times.

Two ways of populating the excited state based on LCT were investigated: the former used as the target observable the population of the excited state described by the projection operator onto the excited state (i.e., $O^=P^2=P21^=P2$⁠), while the latter employed as the target observable the molecular energy described by the unperturbed molecular Hamiltonian (i.e., $O^=H^0$⁠). To show that the monotonic evolution of the target observable $⟨O^⟩ψt$ is guaranteed only if $[O^,H^0]=0$⁠, we also compare the results obtained from control calculations in the presence of nonadiabatic couplings (where $[P^2,H^0]≠0$ and $[H^0,H^0]=0$⁠) and in the absence of nonadiabatic couplings (where both target operators $P^2$ and $H^0$ commute with $H^0$⁠). The control calculations were performed by solving the nonlinear time-dependent Schrödinger equation (1) with the implicit midpoint algorithm combined with the dynamic Fourier method^21–24 for a total time of t_f = 256 a.u. with a time step Δt = 2⁻³ a.u. In addition, intensity parameters λ = 1.430 × 10⁻² and λ = 1.534 × 10⁻¹ were used for the control of excited-state population $P2(t)=⟨P2⟩ψt$ and molecular energy $E0(t)=⟨H^0⟩ψt$⁠, respectively. These parameters were chosen so that the electric fields of the obtained control pulses were similar during the first period.

Figure 1 shows the excited-state population, molecular energy, and obtained control pulse for the control of either the excited-state population (left panels) or molecular energy (right panels). In the figure, the results obtained in the presence and in the absence of nonadiabatic couplings are also compared for each target. The population and energy control schemes result in similar population dynamics, and in both schemes, the population of the excited state reaches 0.99 at time t_f. The carrier frequencies of the obtained control pulses are, as expected, similar and correspond to the electronic transition between the two electronic states of the model. As predicted, when controlling the excited-state population $⟨P2⟩ψt$ in the presence of nonadiabatic couplings given by Eq. (44), the evolution of the population is not monotonic (see the solid line in the inset of the top left panel of Fig. 1) because the control operator does not commute with the molecular Hamiltonian. Compare this with the monotonic increase in the population when the nonadiabatic couplings are zero [dotted line in the same inset; $V^12(qc)=V^21(qc)=0$], and the target operator commutes with the molecular Hamiltonian. In contrast, when controlling the molecular energy $⟨H^0⟩ψt$⁠, its time evolution is always monotonic because the molecular Hamiltonian commutes with itself, whether the nonadiabatic couplings are included or not (see the inset of the middle right-hand panel of Fig. 1). Because increasing the population of the excited state has almost the same effect as increasing the molecular energy, very similar dynamics and control pulses are obtained. Yet, the energy and population controls do not always yield similar results. In the retinal model, when performing energy control, no vibrational energy is pumped into the system because the diagonal terms of the electric-dipole moment operator are all zero by construction (hence $⟨[μ^,T^]⟩ψt=0$⁠). Consequently, only electronic potential energy is added to the system, and the corresponding control pulse is similar to the one obtained from the population control.

To verify the orders of convergence predicted in Sec. III B, we performed convergence analysis of control simulations using various integrators. Simulations with each integrator were repeated with different time steps, and the resulting wavefunctions at the final time t_f were compared. As a measure of the convergence error, we used the L₂-norm $‖ψtf(Δt)−ψtf(Δt/2)‖$⁠, where ψ_t(Δt) denotes the state at time t obtained after propagation with the time step Δt. Figure 2 displays the convergence behavior of both Euler methods, approximate explicit TVT split-operator algorithm, trapezoidal rule, and proposed implicit midpoint method as well as its symmetric compositions, when controlling the excited-state population. Note that all integrators have their predicted orders of convergence. The approximate explicit TVT split-operator algorithm is, for the reasons mentioned in Sec. III E, only first-order and not second-order as one might naively expect. For the convergence of other simulations, we refer the reader to Figs. S4–S6 of the supplementary material. Together, these results confirm that both population and energy control follow the predicted order of convergence regardless of the presence of nonadiabatic couplings.

Because the higher-order methods require more work to perform each step, a higher order of convergence may not guarantee higher efficiency. Therefore, we evaluated the efficiency of each method directly by measuring the computational cost needed to reach a prescribed convergence error. Figure 3 shows the convergence error as a function of the central processing unit (CPU) time and confirms that, except for very crude calculations, higher-order integrators are more efficient than any of the first- and second-order methods. For example, to reach errors below a rather high threshold of 3 × 10⁻⁴, the fourth-order integrator obtained with the Suzuki composition scheme is already more efficient than any of the first- or second-order algorithms. The efficiency gain increases further when highly accurate results are desired. Indeed, for an error of 10⁻⁹, the eighth-order optimal method is 48 times faster than the basic, second-order implicit midpoint method and ∼400 000 times faster than the approximate explicit TVT split-operator algorithm (for which, due to its inefficiency, the speedup had to be estimated by linear extrapolation). High accuracy is hard to achieve with the explicit methods because both the explicit Euler and approximate explicit TVT split-operator algorithms have only first-order convergence. Note also that the cost of implicit methods is not a monotonous function of the error because the Newton–Raphson method needs more iterations to converge for larger than for smaller time steps. Indeed, for time steps (or errors) larger than a critical value, the CPU time might in fact increase with further increasing time step (or error). The efficiency plots of other control simulations (see Figs. S7–S9 of the supplementary material) confirm that the increase in efficiency persists regardless of the control target (energy or population) and presence or absence of nonadiabatic couplings.

Figure 4 analyzes how the time reversibility and norm conservation depend on the time step. The figure confirms that all proposed integrators are exactly time-reversible and norm-conserving regardless of the time step (the slow increase in the error with decreasing time step is due to the accumulation of numerical roundoff errors because smaller time steps require a larger number of steps to reach the same final time t_f). In contrast, Fig. 4 shows that an unrealistically small time step would be required for the trapezoidal rule and both Euler methods to conserve the norm exactly, and for the explicit split-operator algorithm and both Euler methods to be exactly time-reversible. Figures S10–S12 of the supplementary material confirm that neither the chosen control objective nor the nonadiabatic couplings influence the geometric properties of the integrators.

We also checked how the conservation of geometric properties by the integrators depends on time t for a fixed time step. For these simulations, we used a greater final time t_f = 2048 a.u. and an intentionally large time step Δt = 2⁻² a.u. The grid was modified to 256 points between θ = ±3π/2 a.u. and 64 points between q_c = ±9 a.u., ensuring that the grid representation of the wavefunction at the final time t_f was converged (see Fig. S2 of the supplementary material). Figure 5 displays the time evolution of the geometric properties for the elementary integrators (i.e., the trapezoidal rule, implicit midpoint, approximate explicit split-operator, and both Euler methods). Figure 5(a) confirms that the implicit midpoint method and the approximate explicit TVT split-operator algorithm conserve the norm exactly (i.e., to machine precision) even though a large time step was used. In contrast, the trapezoidal rule and both Euler methods do not conserve the norm. Figure 5(b) shows that only the trapezoidal rule and the implicit midpoint method are time-reversible. However, due to the nonlinearity of the Schrödinger equation (1) and the accumulation of roundoff errors, the time reversibility of these integrators slowly deteriorates as time increases. (For a more detailed analysis of this gradual loss of time reversibility, we refer the reader to Sec. V of the supplementary material.) Figures 5(c)–5(e) confirm that even the implicit midpoint method does not conserve the inner product, distance between two states, and total energy; this is not surprising because, due to nonlinearity, the exact solution does not conserve these properties either. Figures S13–S15 of the supplementary material also confirm that neither the chosen control objective nor the nonadiabatic couplings influence the time evolution of the geometric properties of these integrators.

We presented high-order time-reversible integrators for the nonlinear time-dependent Schrödinger equation and demonstrated their efficiency and geometric properties on the problem of local control of quantum systems. The basic time-reversible integrator is an adaptation of the implicit midpoint method to the nonlinear Schrödinger equation and is obtained by composing the explicit and implicit Euler methods. It is norm-conserving, symmetric, time-reversible, and of second order of accuracy in the time step. Because it is symmetric, the implicit midpoint method can be composed using symmetric composition methods to obtain integrators of an arbitrary even order of accuracy. These higher-order integrators conserve all the properties of the original second-order method.

In contrast, the explicit TVT split-operator algorithm is generally only an approximate adaptation of the standard second-order TVT split-operator algorithm to the nonlinear Schrödinger equation, which results from LCT. Because this integrator is not implicit, it is only of first-order accuracy in the time step and loses time reversibility while still conserving the norm. The trapezoidal rule, another popular algorithm for solving the Schrödinger equation, remains symmetric and time-reversible, but does not conserve the norm of the wavefunction propagated with a nonlinear Schrödinger equation.

Although we applied the proposed algorithms only to the special case of LCT, they should be useful for any nonlinear time-dependent Schrödinger equation if high accuracy, norm conservation, and time reversibility of the solution are desired.

See the supplementary material for a detailed description of the generalized minimal residual method, convergence of the wavefunction with respect to the grid density, a plot of the potential energy surfaces, an additional analysis of the convergence, efficiency, and geometric properties of various integrators, and a detailed analysis of the gradual loss of time reversibility of the implicit midpoint method.

The authors thank Seonghoon Choi for useful discussions and acknowledge financial support from the Swiss National Science Foundation within the National Center of Competence in Research “Molecular Ultrafast Science and Technology” (MUST) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 683069 – MOLEQULE).

The data that support the findings of this study are available within this article and its supplementary material.

Here, we verify which geometric properties of the exact evolution are preserved by various integrators. The analysis generalizes the analysis from Appendix A of Ref. 36 for the linear to the nonlinear Schrödinger equation. To simplify the proofs, wherever it is not ambiguous, we shall use the abbreviated notation $Ûappr(ψ)≔Ûappr(t+Δt,t;ψ)$ for the evolution operator for a single time step and ϵ ≔ Δt/ℏ for the time step divided by Planck’s constant.

In general, the norm is conserved if and only if

which follows from a derivation analogous to Eq. (9) for the exact operator $Û(ψ)$⁠. As in the linear case, neither Euler method conserves the norm because

and

The trapezoidal rule, norm-conserving in the linear case, loses this property for nonlinear Hamiltonians since

and the last nontrivial expression does not reduce to 1 because $Ĥ(ψt)≠Ĥ(ψt+Δt)$⁠. In contrast, both the implicit midpoint and approximate explicit TVT split-operator algorithms conserve the norm even in the nonlinear setting because