Time-reversible and norm-conserving high-order integrators for the nonlinear time-dependent Schr\"{o}dinger equation: Application to local control theory

The explicit split-operator algorithm has been extensively used for solving not only linear but also nonlinear time-dependent Schr\"{o}dinger equations. When applied to the nonlinear Gross-Pitaevskii equation, the method remains time-reversible, norm-conserving, and retains its second-order accuracy in the time step. However, this algorithm is not suitable for all types of nonlinear Schr\"{o}dinger equations. Indeed, we demonstrate that local control theory, a technique for the quantum control of a molecular state, translates into a nonlinear Schr\"{o}dinger equation with a more general nonlinearity, for which the explicit split-operator algorithm loses time reversibility and efficiency (because it has only first-order accuracy). Similarly, the trapezoidal rule (the Crank-Nicolson method), while time-reversible, does not conserve the norm of the state propagated by a nonlinear Schr\"{o}dinger equation. To overcome these issues, we present high-order geometric integrators suitable for general time-dependent nonlinear Schr\"{o}dinger equations and also applicable to nonseparable Hamiltonians. These integrators, based on the symmetric compositions of the implicit midpoint method, are both norm-conserving and time-reversible. The geometric properties of the integrators are proven analytically and demonstrated numerically on the local control of a two-dimensional model of retinal. For highly accurate calculations, the higher-order integrators are more efficient. For example, for a wavefunction error of $10^{-9}$, using the eighth-order algorithm yields a $48$-fold speedup over the second-order implicit midpoint method and trapezoidal rule, and $400000$-fold speedup over the explicit split-operator algorithm.


I. INTRODUCTION
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.][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.][15][16][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,19hese methods are of low accuracy (in time and/or space) and do not always preserve the geometric properties of the exact solution. 20For 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][22][23][24] commonly used for the lina) Electronic mail: julien.roulet@epfl.chb) Electronic mail: jiri.vanicek@epfl.chear 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, timereversible only because the ordinary differential equation that must be solved when propagating the molecular state with the potential part of the Hamiltonian leaves the nonlinear term invariant in time. 18We show here that for nonlinear terms of 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., 25 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) an expectation value of a specified operator.LCT has been successfully used to control various processes such as energy and population transfer, [25][26][27][28] dissociation and association dynamics, [29][30][31][32] direction of rotation in molecular rotors 33 , and electron transfer. 34Controlling quantum systems using LCT changes the nature of the timedependent 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.
A paradigm of a nonlinear Schrödinger equation is the Gross-Pitaevskii equation, 13,14,17 in position representation expressed as ∇ 2 ψ t (q) + V (q)ψ t (q) + C|ψ t (q)| 2 ψ t (q), 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. 20Many other examples are provided by the Dirac-Frenkel variational principle, 1,2 which approximates the exact solution of a linear Schrödinger equation with 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 Ĥ(ψ t ) = P (ψ t ) Ĥ, where P (ψ t ) is the projection operator on the tangent space to the approximation manifold at the point ψ t . 4,38Note 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, timeindependent basis and the approximation manifold is a vector space. 4)

A. Geometric properties of the exact evolution operator
With initial condition |ψ t0 and assuming that t ≥ t 0 , Eq. ( 1) has the formal solution |ψ t = Û (t, t 0 ; ψ)|ψ 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 initial condition Û (t 0 , t 0 ; ψ) = 1.The Hermitian adjoint of Û (t, t 0 ; ψ) is the operator (5) 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 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. 20n operator Û is said to conserve the norm ψ := ψ|ψ 1/2 if Û ψ = ψ .The exact nonlinear evolution operator Û ≡ Û (t, t 0 ; ψ) conserves the norm because where we used relation (5).
In the theory of dynamical systems, an adjoint Û (t, t 0 ; ψ) * of Û (t, t 0 ; ψ) is defined as the inverse of the evolution operator taken [in contrast to the Hermitian adjoint (5)] with a reversed time flow: 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.

B. Nonlinear character of local control theory
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.
equal to the sum of the Hamiltonian Ĥ0 of the system in the absence of the field and the interaction potential Vint (t).Within the electric-dipole approximation, 39 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 initial condition |ψ t0 leads to the formal solution |ψ t = Û (t, t 0 )|ψ t0 with the exact evolution operator given by the time-ordered exponential 1][42] Because it is unitary, the evolution operator conserves the norm as well as the inner product and symplectic structure. 4However, since the Hamiltonian is time-dependent, the Schrödinger equation ( 16) is a nonautonomous differential equation, 42 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 which 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 we can control monotonically.
The left-hand side of Eq. (21) suggests that the control field can be either viewed 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. 42By acknowledging this nonlinear character, the interaction potential from Eq. (18) becomes and Eq. ( 16) becomes an example of a nonlinear timedependent Schrödinger equation (1) with Hamiltonian operator Ĥ(ψ) := Ĥ0 + VLCT (ψ).

III. GEOMETRIC INTEGRATORS FOR THE NONLINEAR SCHR ÖDINGER EQUATION
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 Appendix A.

A. Loss of geometric properties by Euler methods
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.

B. Recovery of geometric properties and increasing accuracy by composition
Composing the implicit and explicit Euler methods yields, depending on the order of composition, either the implicit midpoint method Ûmid (t + ∆t, t; ψ t+∆t/2 ) := Ûexpl (t + ∆t, t + ∆t/2; ψ t+∆t/2 ) or the trapezoidal rule (or Crank-Nicolson method) Ûtrap (t + ∆t, t; ψ) Both methods are second-order, symmetric, and timereversible regardless of the size of the time step. 36,42Although the trapezoidal rule conserves the norm of a state evolved with a time-independent linear Hamiltonian, 36 it loses this property when the Hamiltonian is timedependent 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][43][44][45][46] in order to obtain integrators of 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, γ 1 , . . ., γ M are real composition coefficients which satisfy the relations and a more-involved third condition 42 guaranteeing the increase in the order of accuracy.
The simplest composition methods that are used here are the triple-jump 43 and Suzuki's fractal 44 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.

C. Solving the implicit step in a general nonlinear 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 ψ (0) is close enough to the exact solution of the implicit propagation, the Newton-Raphson iteration ( 31) is a contraction mapping and by the fixedpoint 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][48][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 a 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 ψ (0) using the explicit Euler method [see Eq. ( 25)].Choose an error threshold ε and a maximum iteration number m.
6. End Do.The algorithm fails when k = m and no approximate solution has been found.

D. Solving the implicit step in LCT
In the case of LCT, ÛLCT,impl (ψ and the Jacobian of the nonlinear functional (30) is To obtain the third row of Eq. ( 34), we used , where the generalized complex derivative 52 of the interaction potential is given by the bra vector E. Approximate application of the explicit split-operator algorithm to the nonlinear Schrödinger equation The algorithms that we described above apply to Hamiltonians that are not only nonlinear but also nonseparable, i.e., to Hamiltonians Ĥ which cannot be written as a sum Ĥ = T (p)+V (q) of a momentum-dependent kinetic term and 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 + Vtot (ψ), where Vtot (ψ) := V0 + VLCT (ψ) is the sum of the system's and interaction potential energy operators.It is, therefore, tempting to use the splitoperator 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 ÛTV (t + ∆t; t, ψ t ) := e −i T ∆t/ e −i Vtot(ψt)∆t/ (36)   in the TV split-operator algorithm and by ÛVT (t+∆t; t, ψ t+∆t ) := e −i Vtot(ψt+∆t)∆t/ e −i T ∆t/ (37) in the VT split-operator algorithm.These integrators are norm-conserving but only first-order and timeirreversible.From their definitions ( 36) and ( 37), it follows immediately that the TV and VT algorithms are adjoints 35 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 ÛTVT (t + ∆t; t, ψ t+∆t/2 ) := ÛTV (t + ∆t; t + ∆t/2, ψ t+∆t/2 ) × ÛVT (t + ∆t/2; t, ψ t+∆t/2 ), (38)   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)- (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 −i T ∆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 splitoperator algorithm Ûexpl TVT (t + ∆t; t, ψ t, T ∆t/2 ) := ÛTV (t + ∆t; t + ∆t/2, ψ t, T ∆t/2 ) × ÛVT (t + ∆t/2; t, ψ t, T ∆t/2 ).(39)   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 ÛV T algorithm exactly and, therefore, cannot be time-reversible.Yet, the approximate explicit algorithm does conserve norm (see Appendix A for proofs of the geometric properties).

F. Dynamic Fourier method
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][22][23][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 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 Fourier transform.In the numerical examples below, we used the Fastest Fourier Transform in the West 3 (FFTW3) library 53 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).

IV. NUMERICAL EXAMPLES
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 potential energy operator with components Here (all parameters are in eV units), ω = 0.19 is the vibrational frequency of the coupling mode, m −1 = 4.84 • 10 −4 is the inverse mass of the reaction coordinate, W 1 = 3.6 and W 2 = 1.09 determine the depth of the well in the reaction coordinate for the ground and excited electronic states, respectively, χ 2 = 0.1 is the gradient of the linear perturbation in the excited electronic state, E 2 = 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 = E LCT .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 to simulate the photo-excitation step of the photo-isomerization of the retinal molecule.We therefore use as initial state ψ 0 the ground vibrational state of the harmonic fit of the ground potential energy surface [i.e., a two-dimensional Gaussian wavepacket with q 0 = p 0 = (0, 0), and σ 0 = (0.128, 1) a.u.] with initial populations P 1 (0) = 0.999 and P 2 (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., Ô = P2 = P 2 1 = P 2 ), while the latter employed as the target observable the molecular energy described by the unperturbed molecular Hamiltonian (i.e., Ô = Ĥ0 ).To show that the monotonic evolution of the target observable Ô ψt is guaranteed only if [ Ô, Ĥ0 ] = 0, we also compare the results obtained from control calculations in the presence of nonadiabatic couplings (where [ P2 , Ĥ0 ] = 0 and [ Ĥ0 , Ĥ0 ] = 0) and in the absence of nonadiabatic couplings (where both target operators P2 and Ĥ0 commute with Ĥ0 ).The control calculations were performed by solving the nonlinear timedependent Schrödinger equation ( 1) with the implicit midpoint algorithm combined with the dynamic Fourier method [21][22][23][24] for a total time t f = 256 a.u. with a time step ∆t = 2 −3 a.u.In addition, intensity parameters λ = 1.430 × 10 −2 and λ = 1.534 × 10 −1 were used for the control of excited-state population P 2 (t) = P 2 ψt and molecular energy E 0 (t) = Ĥ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 P 2 ψ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 of the population when the nonadiabatic couplings are zero [dotted line in the same inset; V12 (q c ) = V21 (q c ) = 0] and the target operator commutes with the molecular Hamiltonian.In contrast, when controlling the molecular energy Ĥ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 2 -norm ψ t f (∆t) − ψ t f (∆t/2) , where ψ t (∆t) denotes the state at time t obtained after propagation with time step ∆t. Figure 2 displays the convergence behavior of both Euler methods, approximate explicit TVT split-operator algorithm, trapezoidal rule, and the proposed implicit midpoint method as well as its symmetric compositions, when controlling the excited state popu- lation.Notice 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 naïvely expect.For the convergence of other sim-ulations, 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, higherorder integrators are more efficient than any of the firstand second-order methods.For example, to reach errors below a rather high threshold of 3 × 10 −4 , the fourthorder 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 −9 , the eighth-order optimal method is 48 times faster than the basic, second-order implicit midpoint method and approximately 400000 times faster than the approximate explicit TVT split-operator algorithm (for which, due to its inefficiency, the speedup had FIG. 4. Norm conservation (left) and time reversibility (right) of various integrators at the final time t f as a function of the time step ∆t used for the local control of population in the presence of nonadiabatic couplings.Time reversibility is measured by the distance between the initial state ψ0 and a "forward-backward" propagated state ψ0 := Û (0, t) Û (t, 0)ψ0 [Eq.(11)] and line labels are the same as in Fig. 2.
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 splitoperator algorithms have only first-order convergence.Notice 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 timereversible and norm-conserving regardless of the time step (the slow increase of the error with decreasing time step is due to the accumulation of numerical roundoff errors because a smaller time steps requires a larger number FIG. 5. Geometric properties of various integrators used for the local population control in the presence of nonadiabatic couplings.Panel (a) shows that only the implicit midpoint and approximate explicit split-operator methods conserve the norm, while panel (b) demonstrates that only the implicit midpoint method and the trapezoidal rule are timereversible.(Reversibility is measured as in Fig. 4).Bottom three panels show that no method conserves (c) the inner product, (d) distance between two states (which would imply stability), or (e) total energy Etot(t) := E0(t) + VLCT(ψt) ψ t because even the exact evolution operator does not preserve these properties.State φ0 is ψ0 displaced along the reaction coordinate, i.e., a Gaussian wavepacket with parameters q0 = (0.1, 0), p0 = (0, 0), and σ0 = (0.128, 0) a.u.The time step ∆t = 2 −2 a.u. was used for all calculations and line labels are the same as in Fig. 2. Note that only a few points of the Euler methods are visible in some of the plots because the results of the Euler methods leave the range of these plots very rapidly.
of steps to reach the same final time t f ).In contrast, the figure shows that an unrealistically small time step would be required for the trapezoidal rule and both Euler methods to conserve norm exactly and for the explicit splitoperator algorithm 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 depend 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 −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).The top panel 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.The second panel 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.)The bottom three panels of Fig. 5 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.

V. CONCLUSION
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 of 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.
The data that support the findings of this study are contained in the paper and supplementary material.

SUPPLEMENTARY MATERIAL
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.
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 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 (where in the last step we used the commutativity of the middle two factors in the previous expression) and Ûexpl TVT (ψ t, T ∆t/2 ) † Ûexpl TVT (ψ t, T ∆t/2 ) = e i T /2 e i Vtot(ψ t, T ∆t/2 ) e i T /2 × e −i T /2 e −i Vtot(ψ t, T ∆t/2 ) e −i T /2 = 1. (A5)

Symmetry and time reversibility
Neither Euler method is symmetric or time-reversible because they are adjoints of each other: The approximate explicit TVT split-operator algorithm is also time-irreversible because forward propagation is not cancelled by backward propagation: Ûexpl TVT (t; t + ∆t, ψ V T ) Ûexpl TVT (t + ∆t; t, ψ t, T ∆t/2 ) = e i T /2 e i Vtot(ψ V T ) e i T /2 × e −i T /2 e −i Vtot(ψ t, T ∆t/2 ) e −i T /2 = 1, (A7) where ψ V T := e i T /2 ψ t+∆t = e −i Vtot(ψ t, T ∆t/2 ) ψ t, T ∆t/2 denotes the state obtained after forward propagation of ψ t, T ∆t/2 with the potential evolution operator and Vtot (ψ V T ) = Vtot (ψ t, T ∆t/2 ) was used to obtain the last line.It is clear from Eq. (A7) that if the nonlinear term is as in the Gross-Pitaevskii equation, i.e., if V tot (ψ, q) = V 0 (q) + C|ψ(q)| 2 with C a real constant, the approximate explicit TVT split-operator algorithm becomes time-reversible because Vtot (ψ V T ) = Vtot (ψ t, T ∆t/2 ) in that particular case; the reason is that wavefunctions ψ V T (q) and ψ t, T ∆t/2 (q) in position representation only differ by a position-dependent phase factor.However, the equality Vtot (ψ V T ) = Vtot (ψ t, T ∆t/2 ) does not hold for other nonlinear Hamiltonians, such as the one in LCT, which contains a more general nonlinearity.
In contrast, both the implicit midpoint method and trapezoidal rule are always time-reversible because they are symmetric: Ûmid (t + ∆t, t; ψ t+∆t/2 ) * = Ûmid (t, t + ∆t; ψ t+∆t/2 ) −1 and where R m is an m × m upper triangular matrix.Once this factorization is performed, the minimization problem simplifies to y m = argmin y βe 1 − Hm y = argmin y ḡm − Rm y , where ḡm := Q m β 1 = (γ 1 , . . ., γ m+1 ) T and can be decomposed as ḡm = (g m , γ m+1 ) T .The solution of the minimization problem is now trivially obtained by solving the linear triangular system R m y m = g m , which yields Note that it can be proven that the residue after m iterations of the algorithm is given by the last element of ḡm , i.e., r m = |γ m+1 |. 1 Instead of constructing Hm and computing its decomposition (4) after m iterations of Arnoldi's method, both can be performed in an iterative manner, at each iteration j of Arnoldi's method.This makes it possible to obtain ḡj and therefore to check if r j ≤ at each iteration step j, without even solving the minimization problem.
To do so, Hj+1 is obtained by concatenation of Hj and the elements obtained in steps 4 and 7 of Arnoldi's method, i.e., Then, multiplying Hj+1 on the left by the orthogonal matrix obtained at the previous iteration yields an almost upper triangular matrix of dimension (j + 2) × (j + 1).Note that η and ν denote the arbitrary elements j + 1, j + 1 and j + 2, j + 1 of the matrix Rj+1 .From there, Rj+1 , which is upper triangular, is obtained by rotating Rj+1 : where is a rotation matrix of dimension (j + 2) × (j + 2).The complex coefficients of Ω j are selected so that multiplying Ω j with Rj+1 yields the upper triangular matrix Rj+1 .From Eq. ( 10), Q j+1 is obtained by applying the rotation matrix Ω j to Q j as

Figures S1 and S2
show the convergence of the wavefunction with respect to the grid density at the final times t f = 256 a.u. and t f = 2048 a.u., respectively.For both figures, the goal of the calculation was to increase the population of the excited state in the presence of nonadiabatic couplings and the approximate explicit TVT split-operator algorithm was used with a time step ∆t = 2 −1 a.u.The grid convergence was obtained by fixing the number of grid points of one of the two dimensions and repeating the calculation with a varying number of grid points for the other dimension.For each increase in the number of points by a factor of k ≈ √ 2 in the varying dimension, the position and momentum ranges as well as densities in that dimension were increased by a factor √ k ≈ 2 1/4 .Then, the error was measured by the distance 2 ) for the grid convergence of the first and second dimension, respectively, with ψ t f (N 1 , N 2 ) denoting the wavefunction at the final time t f , obtained using grids consisting of N 1 and N 2 points in the first and second dimension, respectively.Note also that X and Y are the fixed numbers of grid points of the first and second dimensions, respectively.The comparison of wavefunctions that are represented on different grids was performed by trigonometric interpolation of the wavefunction represented on the sparser grid.In both figures, exponential convergence with respect to the grid density of each dimension is observed.

IV. CONVERGENCE, EFFICIENCY, AND GEOMETRIC PROPERTIES OF THE INTEGRATORS FOR THE OTHER SIMULATIONS
We present here the results obtained when controlling the population (in the absence of nonadiabatic couplings) as well as the energy (in the presence and the absence of nonadiabatic couplings).

V. LOSS OF TIME REVERSIBILITY DUE TO ROUNDOFF ERRORS
, for the population control in the presence of nonadiabatic couplings, the influence of various parameters on the time reversibility of the implicit midpoint method.This figure aims at demonstrating that the loss of time reversibility at longer times is due to the combination of nonlinearity and accumulation of roundoff errors and errors due to the approximate solution of the nonlinear system described by Eq. ( 30) of the main text.Indeed, when using the implicit midpoint method, this nonlinear system is only solved up to a limited accuracy ε (i.e., the nonlinear system is considered solved if f (ψ (k+1) ) ≤ ε).
Panel (a) shows the influence of the time step on the time reversibility of the implicit midpoint method.For this, we propagated the wavefunction to the final time t f = 2048 a.u. with different time steps and kept the same error threshold ε for solving the nonlinear system.We observe that decreasing the time step decreases the time reversibility of the  implicit midpoint method.This is due to a faster accumulation of roundoff errors caused by the greater number of steps performed when using a smaller time step.Panel (b) depicts the influence of the error threshold ε on the time reversibility.Here, we propagated the wavefunction using the implicit midpoint method using the same parameters but only changed the error threshold ε.As expected, a larger error threshold results in a faster loss of time reversibility as the time increases.Finally, panel (c) shows the influence of the floating-point format precision.We compare the single and double precision floating-point formats for the implicit midpoint method and the approximate explicit TVT split-operator algorithm.The results indicate that decreasing the floating-point format from double to single precision, while keeping other parameters fixed, decreases the time reversibility of the implicit midpoint method while it does not affect the approximate explicit TVT split-operator algorithm.Because the implicit midpoint method is by construction time-reversible, decreasing the floating-point format has the same effect as increasing the threshold ε, namely increasing the roundoff error when solving the nonlinear system, hence decreasing its time reversibility.In contrast, the approximate explicit TVT split-operator algorithm is not affected by the change of precision because it is, by construction, time-irreversible and would remain irreversible even if infinite precision were available.S10).Bottom three panels show that no method conserves (c) the inner product, (d) distance between two states (which would imply stability), or (e) total energy E tot (t) := E 0 (t) + VLCT (ψ t ) ψt because even the exact evolution operator does not preserve these properties.State φ 0 is ψ 0 displaced along the reaction coordinate, i.e., a Gaussian wavepacket with parameters q 0 = (0.1, 0), p 0 = (0, 0), and σ 0 = (0.128, 0) a.u.The time step ∆t = 2 −2 a.u. was used for all calculations and line labels are the same as in Fig. S4.Note that only a few points of the Euler methods are visible in some of the plots because the results of the Euler methods leave the range of these plots very rapidly.FIG.S16.Influence of various parameters on the time reversibility of the implicit midpoint method.The three panels show that the time reversibility of the implicit midpoint method deteriorates (a) if the time step ∆t is decreased, leading to a faster accumulation of roundoff errors, (b) if the error ε threshold of the nonlinear solver is increased, and (c) if the precision of the floating-point arithmetic is decreased.Note that decreasing the floating-point precision does not have a visible influence on the approximate explicit split-operator algorithm, which is irreversible by construction.Time reversibility is measured as in Fig. S10.

FIG. 2 .FIG. 3 .
FIG. 2. Convergence of the molecular wavefunction at the final time t f achieved by the local population control in the presence of nonadiabatic couplings.Top: All studied methods, i.e., explicit and implicit Euler methods, approximate explicit TVT split-operator algorithm, trapezoidal rule, implicit midpoint method and its symmetric compositions.Bottom-left: Methods obtained with the Suzuki composition.Bottom-right: Sixth-order methods obtained with different composition schemes.

)
FIG. S1.Error of the wavefunction at the final time t f = 256 a.u. as a function of the number of grid points.The figure indicates that the numerical grids consisting of N 1 = 128 points between θ = ±π/2 a.u. and N 2 = 64 points between q c = ±9.0a.u.are converged.All calculations were performed using the approximate explicit TVT split-operator algorithm with time step ∆t = 2 −1 a.u.
FIG. S2.Error of the wavefunction at the final time t f = 2048 a.u. as a function of the number of grid points.The figure indicates that the numerical grids consisting of N 1 = 256 points between θ = ±3π/2 a.u. and N 2 = 64 points between q c = ±9.0a.u.are converged.All calculations were performed using the approximate explicit TVT split-operator algorithm with time step ∆t = 2 −1 a.u.

Figure
Figure S3 displays the two diabatic potential energy surfaces.
Figures S4-S6 show the convergence with respect to the time step ∆t of the wavefunction at the final time t f for the different simulations and Figs.S7-S9 show the efficiency of various algorithms.The geometric properties are displayed in Figs.S10-S15.
FIG. S4.Convergence of the molecular wavefunction at the final time t f achieved by the local population control in the absence of nonadiabatic couplings.Top: All studied methods, i.e., explicit and implicit Euler methods, approximate explicit TVT split-operator algorithm, trapezoidal rule, implicit midpoint method and its symmetric compositions.Bottom-left: Methods obtained with the Suzuki composition.Bottom-right: Sixth-order methods obtained with different composition schemes.

FIG. S7 .
FIG. S7.Efficiency of the integrators used for the local population control of retinal in the absence of nonadiabatic couplings.Efficiency is measured by plotting the convergence error as a function of the computational (CPU) cost.
FIG. S10.Norm conservation (left) and time reversibility (right) of various integrators at the final time t f as a function of the time step ∆t used for the local population control in the absence of nonadiabatic couplings.Time reversibility is measured by the distance between the initial state ψ ∆t 0 and a "forward-backward" propagated state ψ∆t 0 := Û (0, t) Û (t, 0)ψ ∆t 0 [Eq.(19)].
FIG. S11.Same as Fig S10, but for energy control in the presence of nonadiabatic couplings.
FIG. S13.Geometric properties of various integrators used for the local population control in the absence of nonadiabatic couplings.Panel (a) shows that only the implicit midpoint and approximate explicit split-operator methods conserve the norm, while panel (b) demonstrates that only the implicit midpoint method and the trapezoidal rule are time-reversible.(Reversibility is measured as in Fig.S10).Bottom three panels show that no method conserves (c) the inner product, (d) distance between two states (which would imply stability), or (e) total energy E tot (t) := E 0 (t) + VLCT (ψ t ) ψt because even the exact evolution operator does not preserve these properties.State φ 0 is ψ 0 displaced along the reaction coordinate, i.e., a Gaussian wavepacket with parameters q 0 = (0.1, 0), p 0 = (0, 0), and σ 0 = (0.128, 0) a.u.The time step ∆t = 2 −2 a.u. was used for all calculations and line labels are the same as in Fig.S4.Note that only a few points of the Euler methods are visible in some of the plots because the results of the Euler methods leave the range of these plots very rapidly.
. A symmetric operator is also time-reversible because propagating a molecular state ψ t0 forward to time t and then backward to time t 0 yields ψ t0 again, i.e., An operator Û (t, t 0 ; ψ) is called symmetric if it is equal to its adjoint, i.e., if Û (t, t 0 ; ψ) = Û (t, t 0 ; ψ) * Û (t 0 , t; ψ) Û (t, t 0 The quantum state |ψ t of a system interacting with an external time-dependent electric field