An implicit split-operator algorithm for the nonlinear time-dependent Schr\"{o}dinger equation

The explicit split-operator algorithm is often used for solving the linear and nonlinear time-dependent Schr\"{o}dinger equations. However, when applied to certain nonlinear time-dependent Schr\"{o}dinger equations, this algorithm loses time reversibility and second-order accuracy, which makes it very inefficient. Here, we propose to overcome the limitations of the explicit split-operator algorithm by abandoning its explicit nature. We describe a family of high-order implicit split-operator algorithms that are norm-conserving, time-reversible, and very efficient. The geometric properties of the integrators are proven analytically and demonstrated numerically on the local control of a two-dimensional model of retinal. Although they are only applicable to separable Hamiltonians, the implicit split-operator algorithms are, in this setting, more efficient than the recently proposed integrators based on the implicit midpoint method.

The best known NL-TDSE is the Gross-Pitaevskii equation, [22][23][24][25][26] which models the dynamics of Bose-Einstein condensates. 27,28To solve this NL-TDSE with cubic nonlinearity, the explicit second-order splitoperator algorithm [29][30][31][32] is frequently used 33 because it is efficient and geometric. 34However, we recently showed 35 that the success of the explicit split-operator algorithm in solving this NL-TDSE is due to its simple nonlinearity and, therefore, is rather an exception than a rule.Indeed, for other nonlinearities, the algorithm becomes time-irreversible and inefficient because its accuracy decreases to the first order in the time step. 35In many applications, this is not an issue, but it can become a problem if accurate wavefunctions are needed.
As we have recently demonstrated, 35 an example of a NL-TDSE, on which the split-operator algorithm loses the second-order accuracy and time reversibility, is proa) Electronic mail: julien.roulet@epfl.chb) Electronic mail: jiri.vanicek@epfl.ch7][38][39][40][41][42][43][44][45][46] LCT is a technique that aims at controlling the expectation value of a specified operator by computing, from the state, an electric field that will either increase or decrease the chosen expectation value.Because the electric field is statedependent, the interaction between this electric field and the system is nonlinear.
To overcome the limitations of the explicit splitoperator algorithm applied to general NL-TDSEs, in our previous work 35 we developed high-order integrators by symmetrically composing the implicit midpoint method.These integrators are applicable to the general nonlinear Schrödinger equation with both separable and nonseparable Hamiltonians and, in contrast to the explicit split-operator algorithm, are efficient, while preserving the geometric properties of the exact solution.
Here, we show that it is not necessary to abandon the split-operator algorithm altogether, but only its explicit nature.In the linear case, the second-order split-operator algorithms are obtained by composing two adjoint firstorder split-operator methods, which are both explicit.We show that to achieve a second-order accuracy in the nonlinear case, one of the two adjoint algorithms must be implicit.Although implicit generalizations of the Verlet algorithm exist [47][48][49][50] for classical systems with nonseparable Hamiltonians, to the best of our knowledge no implicit splitting methods were developed for quantum systems with separable but nonlinear Hamiltonians.Therefore, we present an implicit generalization of the second-order split-operator algorithm, which is geometric, applicable to the general NL-TDSE, and can be composed with various composition methods [51][52][53][54] to further increase its order of convergence and efficiency.
The remainder of this paper is organized as follows: In Sec.II, we present the NL-TDSE, discuss the geometric properties of its evolution operator and describe how LCT generates a NL-TDSE.In Sec.III, we present the algorithms, their geometric properties and the procedure employed to perform the implicit propagation required in the implicit split-operator algorithms.In Sec.IV, we verify the convergence and the geometric properties of the proposed integrators by performing LCT on a twodimensional model of retinal.Section V concludes this work.

II. NONLINEAR SCHR ÖDINGER EQUATION
The nonlinear time-dependent Schrödinger equation describes the time evolution of the molecular state ψ t under the influence of a state-dependent Hamiltonian operator Ĥ(ψ t ).We will assume that the Hamiltonian is separable into a sum of a momentum-dependent kinetic energy operator T (p) and a position-dependent nonlinear potential energy operator V tot (q, ψ), with p and q denoting the momentum and position operators, respectively.

A. Geometric properties of the exact evolution operator
The formal solution of Eq. ( 1) with initial condition ψ t0 can be expressed as where the exact evolution operator Û is given by the time-ordered exponential with T denoting the time-ordering operator.Because the Hamiltonian in Eq. ( 1) is nonlinear, the exact evolution operator (4) is also nonlinear; we emphasize this by including ψ as an explicit argument of Û .When solving Eq. ( 1), some geometric properties conserved by the linear time-dependent Schrödinger equation are not conserved in the nonlinear case.The exact nonlinear evolution operator (4) does not conserve the inner product and, as a result, is not symplectic. 35Furthermore, the energy is not conserved because the state dependence of the Hamiltonian makes it implicitly time-dependent.However, we demonstrated that the exact nonlinear evolution operator conserves the norm and is time-reversible. 35

B. Nonlinear Hamiltonian of local control theory
Local control theory 36,37 aims at controlling the expectation value Ô ψt := ψ t | Ô|ψ t of a chosen operator Ô in the state ψ t .To this end, a state-dependent electric field E LCT (ψ t ) is computed on the fly so that the expectation value Ô ψt increases or decreases monotonously when the electric field interacts with the system.Within the electric-dipole approximation, 55 this interaction is described by the interaction potential where ˆ µ is the electric dipole moment operator.Due to the state-dependent control field, both the operator (5)  and the total potential energy operator Vtot (ψ) := V0 + VLCT (ψ) [see Eq. ( 2)], where V0 denotes the molecular potential energy operator, are nonlinear.To control the expectation value Ô ψt , the control field used in Eq. ( 5) is computed as where λ > 0 is a parameter that scales the amplitude of the control field and the sign is chosen according to whether one wants to increase or decrease Ô ψt .The control field (6) ensures [35][36][37] that the time derivative d Ô ψt /dt remains positive [or negative, depending on the sign in Eq. ( 6)], indicating a monotonic evolution of Ô ψt .However, this monotonic behavior is only guaranteed if the chosen operator Ô commutes with the unperturbed molecular Hamiltonian Ĥ0 := T + V0 , i.e., if [ Ô, Ĥ0 ] = 0. 35,38,44

III. GEOMETRIC INTEGRATORS FOR THE NONLINEAR TIME-DEPENDENT SCHR ÖDINGER EQUATION
To solve the NL-TDSE (1), numerical propagation methods obtain the state ψ t+∆t at time t + ∆t from the state ψ t at time t using the relation where Ûappr (t + ∆t, t; ψ) denotes an approximate evolution operator which depends on ψ and where ∆t is the numerical time step.While all reasonable numerical methods give the exact solution in the limit ∆t → 0, some geometric properties of the exact evolution operator may not be preserved by the numerical methods using a finite ∆t.In this section, we present the different numerical methods and discuss their geometric properties.Detailed proofs of the geometric properties of the presented numerical methods are shown in Appendix A.

A. Loss of geometric properties by the first-order split-operator algorithms
For separable Hamiltonians, the simplest split-step methods are the explicit TV and implicit VT split-operator algorithms, which approximate the exact evolution operator, respectively, as ÛVT (t + ∆t, t; ψ t+∆t ) : = Û Vtot(ψt+∆t) (∆t) Û T (∆t), (9)   where Û Â(∆t) := e −i Â∆t/ denotes an evolution operator associated with a time-independent Hermitian operator Â and time step ∆t.Both integrators ( 8) and ( 9) are norm-conserving.However, both lose the symmetry and time reversibility of the exact evolution operator.Moreover, both integrators are only first-order accurate in the time step, and therefore, very inefficient.Note also that the TV split-operator algorithm is explicit because it depends on the state ψ t while the VT split-operator algorithm is, due to its dependence on the state ψ t+∆t , implicit and, therefore, requires solving a nonlinear system of equations.
The implicit midpoint method is obtained by composing the first-order accurate explicit Ûexpl (t + ∆t, t; ψ t ) := 1 − i Ĥ(ψ t )∆t/ and implicit Ûimpl (t + ∆t, t; ψ t+∆t ) := [1 + i Ĥ(ψ t+∆t )∆t/ ] −1 Euler methods, which are adjoints of each other.For a detailed description of the implicit midpoint and Euler methods in the context of NL-TDSEs, we refer the reader to Ref. 35.
The second-order methods ( 10)-( 12) are all symmetric and time-reversible regardless of the size of the time step.Therefore, they can be further composed using symmetric composition methods 34,35,[49][50][51][52][53][54]56,57 in order to obtain integrators of arbitrary even orders of convergence. To ths end, starting from an integrator Ûp of even order p, an integrator Ûp+2 of order p + 2 is generated using the symmetric composition where ξ n := n j=1 γ j is the sum of the first n real composition coefficients γ j and M denotes the total number of composition steps.Composition coefficients γ 1 , . . ., γ M satisfy the relations M n=1 γ n = 1 (consistency), γ M +1−n = γ n (symmetry), and M j=1 γ p+1 j = 0 (order increase). 50In this work, we will use the triplejump 51 (M = 3) and Suzuki's fractal 52 (M = 5) composition methods, which can both generate integrators of arbitrary even orders of convergence.However, the number of composition steps increases exponentially with the order of convergence, increasing drastically the cost of performing a single time step.To circumvent this, we will also use nonrecursive methods 53,54 for obtaining sixth-eight-and tenth-order integrators.These composition methods, which will be referred to as "optimal", were designed so that they minimize either the sum M n=1 |γ n | or the maximum max n |γ n | of the magnitudes of the composition steps and, therefore, are more efficient than both the triple-jump and Suzuki's fractal.For more details on these composition methods, see Ref. 57.Note that by "order" we mean the formal order because, as shown by Lubich 58 and Thalhammer, 59 who performed rigorous convergence analysis of splitting methods applied to the NL-TDSE, the actual order depends on the regularity of the initial state.Because we do not perform this analysis here, we will verify the predicted (formal) order numerically in Sec.IV.

C. Approximate application of the explicit split-operator algorithm
Because implicit algorithms require more expensive iterative solvers, it is tempting to ignore the implicit character of the above-described VTV and TVT algorithms, and instead employ their explicit versions, which consist of using the state ψ that is available for computing the evolution operator for Vtot .For example, instead of using the state ψ t+∆t in Eq. ( 9) (i.e., performing the implicit propagation exactly for the VT algorithm), the state ψ T ,∆t/2 := Û T (∆t/2)ψ t obtained after the kinetic propagation is often used.After composition with the TV algorithm, this yields the approximate explicit TVT 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 ).5][46] However, despite conserving the norm, the integrator is only first-order accurate and neither symmetric nor time-reversible, as shown in Ref. 35.Indeed, any explicit version of the integrators ( 10) and ( 11) will be first-order accurate and time-irreversible due to ignoring the implicit character of the VT algorithm.

D. Solving the implicit propagation
Both TVT and VTV implicit split-operator algorithms rely on the implicit VT method.Using the evolution operator ÛV T from Eq. ( 9), the implicit VT propagation of a state |ψ t is translated into solving the nonlinear system ÛVT (t + ∆t, t; ψ t+∆t ) −1 |ψ t+∆t = |ψ t . ( This nonlinear system can be written as f (ψ t+∆t ) = 0 with the nonlinear functional where we have, for convenience, included a nonzero factor of Û T (∆t) into the definition of f (ψ).Following Ref. 35, we employed the Newton-Raphson method to solve this nonlinear system.This method computes, until convergence, iterative solutions of the nonlinear system using the iterative map where ψ (k) denotes the solution obtained at the kth iteration and δψ (k) is the state obtained by solving the linear system with Ĵ := δf (ψ)/δψ denoting the Jacobian of the nonlinear mapping f (ψ).The linear system (17) is solved using the generalized minimal residual method, [60][61][62] an iterative method based on the Arnoldi process 63,64 (see the supplementary material of Ref. 35 for a detailed presentation of this algorithm).Similarly to Ref. 35, we employ the solution from the explicit propagation, i.e., the solution obtained using Eq. ( 8), as the initial guess ψ (0) .However, if the initial guess is too far from the implicit solution, which happens at large time steps, the algorithm fails to converge.The procedure described above differs from that presented in Ref. 35 only by the nonlinear system one needs to solve.Fortunately, the use of approximations for estimating the Jacobian is, as in Ref. 35, avoided because the Jacobian Ĵ(ψ) of the nonlinear function ( 15) can be obtained analytically: where we employed, in the third line, the generalized complex derivative 65 of the nonlinear potential, which is given by the bra vector

IV. NUMERICAL EXAMPLES
Integrators presented in Sec.III were tested on a local control simulation in a two-dimensional model describing the cis-trans photo-isomerization of retinal.This ultrafast reaction, which is mediated by a conical intersection, is the first event occurring in the biological process of vision.The model, which uses the reaction coordinate θ, an angle describing the torsional motion of the retinal molecule, and a vibronically active coupling mode q c , was taken from Ref. 66 and we used it as described in Ref. 35 (see Fig. S3 of the supplementary material of Ref. 35 for the two diabatic potential energy surfaces of the model).We used the same grid [a regular direct-product grid consisting of 128 points between θ = ±π/2 a.u. and 64 points between q c = ±9 a.u.], the same initial state [a two-dimensional Gaussian wavepacket ψ 0 (x) = 2 j=1 (σ 2 0,j π) −1/4 exp[ip 0,j (x j − q 0,j )/ − (x j − q 0,j ) 2 /2σ 2 0,j ], with x := (θ, q c ), initial positions and momentum q 0 = p 0 = (0, 0) a.u. and initial width σ 0 = (0.128, 1) a.u., corresponding to the ground vibrational state of the harmonic fit of the ground electronic potential energy surface], and the same initial populations P 1 (0) = 0.999 and P 2 = 0.001 of the ground and excited electronic states, respectively.The kinetic and potential propagations were performed using the dynamic Fourier method [29][30][31][32] with the Fastest Fourier Transform in the West 3 (FFTW3) library 67 to change between position and momentum representations.Following Ref. 35, we assumed that the electric dipole moment operator ˆ µ was coordinate independent (Condon approximation) and aligned with the electric field E LCT .Consequently, we could drop the vector symbols in Eq. ( 5), i.e., replace the vectors ˆ µ and E LCT with the scalars μ and E LCT .In all simulations, the electric dipole moment operator had unit transition (offdiagonal) elements (μ 12 = μ21 = 1 a.u.) and zero diagonal elements (μ 11 = μ22 = 0 a.u.); which allowed the control field to couple the two electronic states and simultaneously avoid coupling between the vibrational states.Throughout this section, the bold font denotes electronic operators expressed as S × S matrices in the basis of S electronic states and that the hat ˆdenotes nuclear operators acting on the Hilbert space of nuclear wavefunctions, i.e., square-integrable functions of D continuous degrees of freedom.
In all simulations, we used LCT for increasing the molecular energy E 0 (t) := Ĥ0 ψt of the system, which required employing the molecular energy operator Ĥ0 as the target observable.First, we performed the local control by solving the NL-TDSE (1) up to the final time t f = 256 a.u.using the implicit TVT split-operator algorithm and intensity parameter λ = 1.534×10 −1 , which was chosen arbitrarily so that the amplitude of the obtained control field was not too high and, at the same time, strong enough to induce a significant increase of molecular energy (see Fig. 1).The results indicate a successful increase in molecular energy [panel (a)] and, as predicted in Sec.II, this increase is monotonic because the molecular energy operator commutes with itself.Indeed, by construction, there are no nonzero diagonal elements in the electric dipole moment operator and, as a result, vibrational energy cannot be added by the control pulse ( [ μ, T] ψt = 0).Instead, due to the presence of nonzero transition (offdiagonal) elements in the electric dipole moment, the control pulse increases the electronic energy of our system by increasing the excited state population P 2 (t) [panel (b)]; this is confirmed by the carrier frequency of the control pulse [panel (c)], which corresponds to an electronic transition between the two states.To verify the order of convergence of the integrators presented in Sec.III, the same simulation was repeated for each integrator with different time steps, and the errors in the obtained wavefunctions were compared at the final time t f = 256 a.u.To measure the convergence er-ror, we used the L 2 -norm ψ t f (∆t) − ψ t f (∆t/2) where ψ t (∆t) denotes the wavefunction at time t obtained using a time step ∆t. Figure 2 shows the convergence behavior of various integrators, including higher-order integrators obtained by composing the implicit TVT method with the triple-jump, Suzuki's fractal, and optimal composition schemes.The results in panel (a) indicate that the implicit TVT method has the expected order of convergence and that it is, for a given time step, more accurate than all first-order methods, including the approximate explicit TVT algorithm.Comparison between different orders of the optimal composition of the implicit TVT algorithm [panel (b)] shows that, for a given time step, a higher order of composition yields more accurate integrators.Similarly, comparing sixth-order methods obtained with different composition schemes [panel (c)] indicates that, for a given time step, Suzuki's fractal composition is more accurate than both the optimal and triple-jump compositions.Note that after reaching a machine precision plateau, the higher-order integrators show a slight increase in the error with a decreasing time step, which is due to the accumulation of roundoff errors since the number of steps increases for a fixed total time of simulation.Moreover, some results for high-order integrators could not be obtained because they did not converge at large time steps (the difference between the initial guess and the implicit solution was too large for the Newton-Raphson method to converge) and became computationally unaffordable at smaller time steps, when the Newton-Raphson method was converging.The higher-order integrators, obtained by composition, require performing many substeps at each time step, which increases their cost.To check that greater accuracy for a given time step is not detrimental to efficiency, in Fig. 3 we plot the dependence of errors of the wavefunctions obtained by various methods on the computational cost, measured by the central processing unit (CPU) time (see also Fig. S2 of the supplementary material, which displays the efficiency results for all studied methods).
Figure 3 demonstrates that, if high accuracy is desired, the higher-order integrators are more efficient even though they require performing many substeps at each time step.For example, below an error of 3 × 10 −4 , the second-order implicit TVT split-operator algorithm is already more efficient than the approximate explicit TVT split-operator algorithm.Figure 3 also shows that the implicit TVT split-operator algorithm is more efficient than the implicit midpoint method, indicating that the implicit split-operator algorithm is the method of choice for separable Hamiltonians and that the implicit midpoint method should only be used when the Hamiltonian is not separable.Indeed, for errors below 7 × 10 −4 , the TVT split-operator algorithm is more efficient than the implicit midpoint method.
In Fig. 4, we checked the preservation of geometric properties by the implicit and approximate explicit TVT methods (see Fig. S3 of the supplementary material for a version of this figure which displays the results for all the elementary methods).Since geometric integrators preserve geometric properties exactly regardless of the size of the time step, we intentionally used a rather large time step ∆t = 2 −2 a.u.We also extended the final time of the simulation to t f = 2048 a.u. in order to induce more dynamics.Following Ref. 35, we modified the grid to 256 points between θ = ±3π/2 a.u. and 64 points between q c = ±9 a.u. to ensure that the grid representation of the wavefunction at the new final time t f remains converged.The results show that while both the implicit and approximate explicit TVT integrators conserve the norm [panel (a)], only the implicit TVT method is timereversible [panel (b)].However, due to the nonlinearity of the time-dependent Schrödinger equation and the accumulation of roundoff errors, one observes a gradual loss of time reversibility as the time increases.(See Sec.V of the supplementary material of Ref. 35 for a detailed analysis of this loss of time-reversibility.)The bottom three panels of Fig. 4 (and Fig. S3) demonstrate that none of the methods conserves the inner product [panel (c)], distance between two states [panel (d)], or total energy [panel (e)], because these properties are not conserved even by the exact nonlinear evolution operator (4).

V. CONCLUSION
We presented high-order integrators for solving the NL-TDSE with separable Hamiltonians.In contrast to their first-order explicit versions, the proposed methods, obtained by composing an implicit split-operator algo- rithm, preserve all geometric properties of the exact solution: they are symmetric, time-reversible, and normconserving.Moreover, the proposed integrators are more efficient than both the explicit split-operator algorithm and the recently proposed 35 compositions of the implicit midpoint method.

SUPPLEMENTARY MATERIAL
The supplementary material contains analogues of Figs.2-5 of the main text that display the numerical results for all studied methods.will always yield a symmetric method. 50Therefore, the proposed high-order integrators are also symmetric and time-reversible.

FIG. 1 .
FIG. 1. Local control simulation whose goal is increasing the molecular energy E0(t).(a) Molecular energy.(b) Excited state population.(c) Pulse obtained by LCT.

FIG. 2 .
FIG. 2. Convergence of the molecular wavefunction at the final time t f = 256 a.u.achieved by the local energy control.(a) First-order and implicit TVT methods.(b) Methods obtained with the optimal composition (Suzuki's fractal is the optimal fourth-order composition scheme 57 ).(c) Sixth-order methods obtained with different composition schemes.

FIG. 3 .
FIG. 3. Efficiency of various integrators used for simulating the local energy control of retinal up to the final time t f = 256 a.u.

FIG. 5 .
FIG. 5. Norm conservation (a) and time reversibility (b) of various integrators at the final time t f = 256 a.u. as a function of the time step ∆t used for the local energy control of retinal.Reversibility is measured as in Fig. 4 and line labels are the same as in Fig. 2.
FIG. S2.Efficiency of the integrators used for simulating the local energy control of retinal up to the final time t f = 256 a.u.Line labels are the same as in Fig. S1.