text:
# Gradient Symplectic Algorithms for Solving the Schrödinger Equation with Time-Dependent Potentials

## I Introduction

Recently, we derived a new class of fourth order algorithms for solving both classical and quantum dynamical problems. These algorithms are based on factorizing the evolution operator $`\mathrm{e}^{ϵ(T+V)}`$ to fourth order with purely positive coefficients and require knowing the gradient of the force in the classical case and the gradient of the potential in the quantum case. The resulting algorithms are symplectic or unitary, respectively. While positive coefficients are absolutely necessary for simulating the diffusion process in Monte Carlo algorithms, or doing imaginary time projections, they are not essential in quantum or classical algorithms. Nevertheless, we have shown that this class of gradient symplectic algorithms is far superior to existing fourth order algorithms with negative coefficients. In this work, using Suzuki’s method of implementing operator time-ordering, we prodouce a class of even more effective algorithms for solving quantum dynamical problems with explicit time-dependent potentials. Despite the vast literature on this subject , we believe our work has initiated a new direction in algorithm development. In the past, one labors assiduously to avoid higher order commutators. Here, we show that their inclusion can yield algorithms of great efficiency. Our algorithm 4A, to be discuss below, is the fastest fourth order algorithm known, needing only four Fast Fourier Transforms (FFT) per iteration. This is only twice the computational effort of the second order split-operator method, but the algorithm can converge at time step sizes an order of magnitude larger. Our optimized algorithms, when compared on an equal effort basis, have fourth order error coefficients that are three orders of magnitude smaller than Forest-Ruth’s algorithm and a factor of 30 smaller than McLachlan’s algorithm; both are fourth order algorithms with negative coefficients.

While these gradient symplectic algorithms are very efficient when the potential gradient is known analytically, they remain equally effective when the gradient is obtained numerically. These algorithms are of particular interest in solving the time-dependent Schrödinger in a large 3D mesh. In 3D, even for a modest grid size of (256)3, the number of mesh points already exceed 16 millions. If the wave function array is double precision and complex, its storage alone would have required 268MB. For such a large number of grid points, any vector-matrix multiplication would be prohibitively expensive and must be avoided. For our algorithms, the costliest computational step is just the use of FFT. The unitary character and the large time step acceptance of these algorithms make them ideal for doing long time quantum simulations.

The key problem in solving the Schrödinger equation with time-dependent potentials is the time ordering of operators. This problem is solved variously in the literature by transforming it into a classical problem, treating time as another “spatial coordinate”, introducing auxiliary variables, etc.. Most end up with some sort of time derivative operator, but none has the simplicity of Suzuki’s method of directly implementing time-ordering via a forward time derivative operator. No time integration is necessary. Since this work is less accessible, but is of special relevance to the operator factorization approach of deriving algorithms, we summarize it in some detail in the next section. In Section III, we apply Suzuki’s method and derive four gradient symplectic algorithms for solving the time-dependent Schrödinger equation. In Section IV, we use these algorithms to solve the Walker-Preston model and compare their convergent properties with existing algorithms. In Section V, we derive one-parameter families of these algorithms and show that they can be further optimized for specific applications. Section VI summarizes our conclusions.

## II Operator Decomposition of Ordered Exponentials

For $`H(t)`$ a time-dependent operator, the evolution equation

$$\frac{}{t}\psi (t)=H(t)\psi (t),$$
(1)
has the operator solution

$$\psi (t+\mathrm{\Delta }t)=T\left(\mathrm{exp}_t^{t+\mathrm{\Delta }t}H(s)𝑑s\right)\psi (t).$$
(2)
The time-ordered exponential not only has the conventional expansion

$$T\left(\mathrm{exp}_t^{t+\mathrm{\Delta }t}H(s)𝑑s\right)=1+_t^{t+\mathrm{\Delta }t}H(s_1)𝑑s_1+_t^{t+\mathrm{\Delta }t}𝑑s_1_t^{s_1}𝑑s_2H(s_1)H(s_2)+\mathrm{},$$
(3)
but also the more intuitive interpretation

$`T\left(\mathrm{exp}{\displaystyle _t^{t+\mathrm{\Delta }t}}H(s)𝑑s\right)=`$ $`\underset{n\mathrm{}}{lim}T\left(\mathrm{e}^{\frac{\mathrm{\Delta }t}{n}_{i=1}^nH(t+i\frac{\mathrm{\Delta }t}{n})}\right),`$ (4)
$`=`$ $`\underset{n\mathrm{}}{lim}\mathrm{e}^{\frac{\mathrm{\Delta }t}{n}H(t+\mathrm{\Delta }t)}\mathrm{}\mathrm{e}^{\frac{\mathrm{\Delta }t}{n}H(t+\frac{2\mathrm{\Delta }t}{n})}\mathrm{e}^{\frac{\mathrm{\Delta }t}{n}H(t+\frac{\mathrm{\Delta }t}{n})}.`$ (5)
There are many ways of solving the time-ordering problem. For this work on operator factorization, we prefer Suzuki’s method, which directly implements time ordering without any additional formalism or auxiliary variables.

Let $`D`$ denotes the forward time derivative operator

$$D=\frac{\stackrel{}{}}{t}$$
(6)
such that for any two time-dependent functions $`F(t)`$ and $`G(t)`$,

$$F(t)\mathrm{e}^{\mathrm{\Delta }tD}G(t)=F(t+\mathrm{\Delta }t)G(t).$$
(7)
Suzuki’s proved that

$$T\left(\mathrm{exp}_t^{t+\mathrm{\Delta }t}H(s)𝑑s\right)=\mathrm{exp}[\mathrm{\Delta }t(H(t)+D)].$$
(8)
Using the more intuitive definition of time-ordering (5), and invoking Trotter’s formula, one proof only requires two lines:

$`\mathrm{exp}[\mathrm{\Delta }t(H(t)+D)]=`$ $`\underset{n\mathrm{}}{lim}(\mathrm{e}^{\frac{\mathrm{\Delta }t}{n}H(t)}\mathrm{e}^{\frac{\mathrm{\Delta }t}{n}D})^n,`$ (9)
$`=`$ $`\underset{n\mathrm{}}{lim}\mathrm{e}^{\frac{\mathrm{\Delta }t}{n}H(t+\mathrm{\Delta }t)}\mathrm{}\mathrm{e}^{\frac{\mathrm{\Delta }t}{n}H(t+\frac{2\mathrm{\Delta }t}{n})}\mathrm{e}^{\frac{\mathrm{\Delta }t}{n}H(t+\frac{\mathrm{\Delta }t}{n})},`$ (10)
where property (7) has been applied repeatedly and accumulatively.

For the widely applicable case of $`H(t)=T+V(t)`$, where only one of the operator is explicitly dependent on time, the short time evolution of (2) can be written using (8) as

$$\psi (t+\mathrm{\Delta }t)=\mathrm{e}^{\mathrm{\Delta }t[\stackrel{~}{T}+V(t)]}\psi (t),$$
(11)
which is just like the time-independent case but with an effective $`\stackrel{~}{T}=T+D`$. This suggests a two-step approach of deriving time-dependent algorithms. First, decompose $`\mathrm{e}^{\mathrm{\Delta }t[\stackrel{~}{T}+V(t)]}`$ in terms of $`\mathrm{e}^{\mathrm{\Delta }t\stackrel{~}{T}}`$ and $`\mathrm{e}^{\mathrm{\Delta }tV(t)}`$ using any factorization scheme applicable in the time-independent case. Next, since $`[T,D]=0`$, factorize exactly

$$\mathrm{e}^{\mathrm{\Delta }t\stackrel{~}{T}}=\mathrm{e}^{\mathrm{\Delta }tD}\mathrm{e}^{\mathrm{\Delta }tT},$$
(12)
and incorporate all time-dependent requirements by applying (7). For example, a second order factorization of (11) gives,

$`𝒯_A^{(2)}`$ $`=\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }t\stackrel{~}{T}}\mathrm{e}^{\mathrm{\Delta }tV(t)}\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }t\stackrel{~}{T}},`$ (15)
$`=\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tD}\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tT}\mathrm{e}^{\mathrm{\Delta }tV(t)}\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tD}\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tT},`$
$`=\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tT}\mathrm{e}^{\mathrm{\Delta }tV(t+\mathrm{\Delta }t/2)}\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tT},`$
which is the well known midpoint algorithm for time-dependent problems. The other second order factorization gives the alternative second order algorithm,

$`𝒯_B^{(2)}`$ $`=\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tV(t)}\mathrm{e}^{\mathrm{\Delta }t\stackrel{~}{T}}\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tV(t)},`$ (18)
$`=\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tV(t)}\mathrm{e}^{\mathrm{\Delta }tD}\mathrm{e}^{\mathrm{\Delta }tT}\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tV(t)},`$
$`=\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tV(t+\mathrm{\Delta }t)}\mathrm{e}^{\mathrm{\Delta }tT}\mathrm{e}^{\frac{1}{2}\mathrm{\Delta }tV(t)}.`$
Thus, for $`H(t)=T+V(t)`$, the effect of time-ordering is to increment the time-dependence of each potential operator $`V(t)`$ by the sum of time steps of all the $`T`$ operators to its right.

For the Schrödinger equation with a time-dependent potential, the wave function is evolved forward in a short time $`\mathrm{\Delta }t`$ by

$$\psi (ϵ)=\mathrm{e}^{ϵ[\stackrel{~}{T}+V(t)]}\psi (0),$$
(19)
where $`ϵ=i\mathrm{\Delta }t`$, $`\stackrel{~}{T}=\stackrel{~}{D}+T`$, $`\stackrel{~}{D}=iD`$,

$$T=\frac{1}{2\mu }\underset{i}{}\frac{^2}{x_i^2}\text{ and }V(t)=V(x_i,t).$$
(20)
We will work in atomic units such that the kinetic energy operator has this standard form. Moreover, to do away with messy notations involving $`i`$, we will use $`ϵ`$ as the time step variable everywhere. When $`ϵ`$ appears as the argument of the wave function or potential, it is to be understood that it denotes only the real time step variable $`\mathrm{\Delta }t`$ without the factor $`i`$. (In this way, algorithms can be directly applied to the classical case with $`ϵ`$ purely real without change in form.) For conciseness, we will always regard the present time step as time zero. Thus, the two second order algorithms for solving the Schrödinger equation with step size $`\mathrm{\Delta }t`$ can be denoted simply as

$`𝒯_A^{(2)}(ϵ)`$ $`=\mathrm{e}^{\frac{1}{2}ϵT}\mathrm{e}^{ϵV(ϵ/2)}\mathrm{e}^{\frac{1}{2}ϵT},`$ (21)
$`𝒯_B^{(2)}(ϵ)`$ $`=\mathrm{e}^{\frac{1}{2}ϵV(ϵ)}\mathrm{e}^{ϵT}\mathrm{e}^{\frac{1}{2}ϵV(0)}.`$ (22)
Since the kinetic energy operator is diagonal in momentum space, the operator $`\mathrm{e}^{ϵT}`$ can be implemented as a vector-vector multiplication in Fourier space. Every occurrence of $`\mathrm{e}^{ϵT}`$ requires two FFTs, one direct to Fourier space for the kinetic energy multiplication, and one inverse back to real space for the potential energy multiplication. To minimize the call for FFTs, one favors algorithms with the fewest occurrence of the kinetic energy operator.

## III Gradient Symplectic Fourth Order Algorithms

Following our two-step approach, we can transcribe any time-independent factorization algorithm into a time-dependent algorithm. For example, the well known Forest-Ruth (FR) algorithm(also discovered independently by Campostrini and Rossi,and Candy and Rozmus) can now be transcribed to solve time-dependent problems as

$$𝒯_{FR}^{(4)}(ϵ)\mathrm{e}^{v_3ϵV(ϵ)}\mathrm{e}^{t_3ϵT}\mathrm{e}^{v_2ϵV(a_2ϵ)}\mathrm{e}^{t_2ϵT}\mathrm{e}^{v_1ϵV(a_1ϵ)}\mathrm{e}^{t_1ϵT}\mathrm{e}^{v_0ϵV(0)},$$
(23)
where, $`s=2^{1/3}`$,

$$v_0=v_3=\frac{1}{2}\frac{1}{2s},v_1=v_2=\frac{1}{2}\frac{s1}{2s},t_1=t_3=\frac{1}{2s},t_2=\frac{s}{2s},$$
(24)
and $`a_1=t_1`$, $`a_2=t_1+t_2`$. For easy identification, we adopt the convention of labelling the time step coefficients of operators $`T`$ and $`V`$ by $`t_i`$ and $`v_i`$ respectively, and denote the intermediate time arguments of $`V`$ by coefficients $`a_i`$. The coefficient of the first operator on the right will be denote by $`v_0`$ or $`t_0`$, followed by paired coefficients of $`v_i`$ and $`t_i`$ for $`i=1,2,\mathrm{}n`$. If one were to decompose $`\mathrm{e}^{ϵ(T+V)}`$ only in terms of $`\mathrm{e}^{ϵT}`$ and $`\mathrm{e}^{ϵV}`$, then Forest-Ruth is the only fourth order algorithm possible with 6 FFTs per iteration. The alternative algorithm with $`T`$ and $`V`$ interchanged, with appropriate modifications of the $`a_i`$ coefficients, is also possible, but would have required 8 FFTs. Notice that some of the coefficients are negative, requiring backward propagation and evaluating the potential at a time prior to the present. This is a consequence of Suzuki’s “no-go” theorem, which proved that beyond second order, $`\mathrm{e}^{ϵ(T+V)}`$ cannot be decomposed into a finite products of $`\mathrm{e}^{t_iϵT}`$ and $`\mathrm{e}^{v_iϵV}`$ with purely positive coefficients $`t_i`$ and $`v_i`$. Thus without exception, all higher order factorization algorithms heretofore proposed in the literature contain negative coefficients.

It is the search for positive coefficient factorizations schemes that led one of us to derive fourth order symplectic algorithms for solving classical and subsequently quantum dynamical problems. To circumvent Suzuki’s “no-go” theorem, these factorization schemes employ an additional operator,

$$[V,[T,V]]=\frac{1}{\mu }\underset{i}{}\left(\frac{V}{x_i}\right)^2,$$
(25)
which is just the square of the gradient of the potential. This can be computed analytically or numerically. For brevity, we will transcribe four previously derived gradient algorithms to their time-dependent form. In the next section, we will describe in more detail the one-parameter family of these algorithms.

Our algorithm 4A (See also Ref.), when applied to the time-dependent case of (19), gives

$$𝒯_A^{(4)}(ϵ)\mathrm{e}^{\frac{1}{6}ϵV(ϵ)}\mathrm{e}^{\frac{1}{2}ϵT}\mathrm{e}^{\frac{2}{3}ϵ\stackrel{~}{V}(ϵ/2)}\mathrm{e}^{\frac{1}{2}ϵT}\mathrm{e}^{\frac{1}{6}ϵV(0)},$$
(26)
with $`\stackrel{~}{V}`$ defined by

$`\stackrel{~}{V}(t)`$ $`=V(t)+{\displaystyle \frac{1}{48}}ϵ^2[V(t),[\stackrel{~}{T},V(t)]],`$ (28)
$`=V(t)+{\displaystyle \frac{1}{48}}ϵ^2[V(t),[T,V(t)]].`$
Note that this is a crucial simplification, because $`[V(t),[D,V(t)]]=0`$ ! Thus the addition of the forward time derivative operator $`D`$ causes no additional complication to the gradient term. ( This is very fortunate, because if one were required to keep the alternative double commutator $`[\stackrel{~}{T},[V(t),\stackrel{~}{T}]]`$, this commutator would have given rise to complicated new terms involving $`V(t)/t`$ and $`^2V(t)/t^2`$! ) With the introduction of the gradient term, algorithm 4A can achieve fourth order acccuracy with only four FFTs. Note also that in order to minimize the evaluation of the gradient term, the double commutator has been placed at the center. This choice seems obvious, but there is intrinsic freedom to redistribute the commutator term among the three potential operator without affecting the fourth order convergence of the algorithm. This is true for all algroithms described below. We will assume that this redistribution can be done when necessary.

Similarly algorithm 4B can be transcribed as

$$𝒯_B^{(4)}(ϵ)\mathrm{e}^{t_2ϵT}\mathrm{e}^{ϵ\frac{1}{2}\overline{V}(a_2ϵ)}\mathrm{e}^{t_1ϵT}\mathrm{e}^{ϵ\frac{1}{2}\overline{V}(a_1ϵ)}\mathrm{e}^{t_0ϵT},$$
(29)
where

$$t_0=t_2=\frac{1}{2}(1\frac{1}{\sqrt{3}}),t_1=\frac{1}{\sqrt{3}},a_1=\frac{1}{2}(1\frac{1}{\sqrt{3}}),a_2=\frac{1}{2}(1+\frac{1}{\sqrt{3}}),$$
(30)
and with $`\overline{V}`$ given by

$$\overline{V}(t)=V(t)+\frac{1}{24}(2\sqrt{3})ϵ^2[V(t),[T,V(t)]].$$
(31)
The time-dependent forms of algorithm 4C and 4D are respectively,

$$𝒯_C^{(4)}(ϵ)\mathrm{e}^{\frac{1}{6}ϵT}\mathrm{e}^{\frac{3}{8}ϵV(5ϵ/6)}\mathrm{e}^{\frac{1}{3}ϵT}\mathrm{e}^{\frac{1}{4}ϵ\stackrel{~}{V}(ϵ/2)}\mathrm{e}^{\frac{1}{3}ϵT}\mathrm{e}^{\frac{3}{8}ϵV(ϵ/6)}\mathrm{e}^{\frac{1}{6}ϵT},$$
(32)
$$𝒯_D^{(4)}(ϵ)\mathrm{e}^{\frac{1}{8}ϵ\stackrel{~}{V}(ϵ)}\mathrm{e}^{\frac{1}{3}ϵT}\mathrm{e}^{\frac{3}{8}ϵV(2ϵ/3)}\mathrm{e}^{\frac{1}{3}ϵT}\mathrm{e}^{\frac{3}{8}ϵV(ϵ/3)}\mathrm{e}^{\frac{1}{3}ϵT}\mathrm{e}^{\frac{1}{8}ϵ\stackrel{~}{V}(0)},$$
(33)
where $`\stackrel{~}{V}(t)`$ is as defined by (28). The number of FFT required by each algorithm is given in Table I.

## IV Solving the Walker-Preston model

To demonstrate the effectiveness of these new algorithms, we use them to solve the Walker and Preston model of a diatomic molecule in a strong laser field. Since this problem has been used by many authors to test their time-dependent algorithms, it is an excellent choice for comparing our algorithms. The model is defined by the one dimensional Hamiltonian (in atomic units),

$$H=\frac{1}{2\mu }\frac{^2}{x^2}+V(x,t),$$
(34)
with

$$V(x,t)=V_0(1\mathrm{e}^{\alpha x})^2+Ax\mathrm{cos}(\omega t),$$
(35)
and where $`\mu =1745`$, $`V_0=0.2251`$, $`\alpha =1.1741`$, $`A=0.011025`$, and $`\omega =0.01787`$. The wave function is initially chosen to be the Morse oscillator ground state, which has the corresponding ground state energy $`E_0=(\omega _0/2)[1\omega _0/8V_0]`$ with $`\omega _0=\alpha \sqrt{2V_0/\mu }`$. In conformity with the above authors, we discretize the wave function $`\psi (x,t)`$ using 64 uniform points at $`x_k=0.8+k\mathrm{\Delta }x`$, with $`\mathrm{\Delta }x=0.08`$. The natural time scale defined by the laser oscillation frequence is $`\tau =2\pi /\omega =351.6`$.

In Fig.1, the percentage errors of the total energy, $`100|E(1000\tau )E_{con}|/E_{con}`$, is plotted as a function of the step size $`\mathrm{\Delta }t`$ used in the calculation. $`E_{con}=5.029155E_0`$ is the converged value at very small time steps. The plotting symbols indicate calculated results. The lines are monomial fits in either $`(\mathrm{\Delta }t)^2`$ or $`(\mathrm{\Delta }t)^4`$. SO is our calculation using the second order split-operator algorithm (22). RS3 is Gray and Verosky’s best convergent results. They use a third order algorithm, but the error is degraded by the Magnus approximation to second order. Both can be well-fitted by a single quadratic $`b_i(\mathrm{\Delta }t)^2`$, with coefficients $`b_{SO}=1.1`$ and $`b_{RS3}=0.67`$. These are plotted as lines running through the data. The convergent range for both is below $`\mathrm{\Delta }t<\tau /100=3.52`$. Algorithm RS3 requires three times as much effort as SO. For the same amount of computational effort, one can run SO three times at $`\mathrm{\Delta }t/3`$ and gain an reduction of $`(1/3)^2`$ in its error coefficient. This projected convergence curve $`(b_{SO}/9)(\mathrm{\Delta }t)^2`$ is plotted as a broken line and labelled as SO. For the same effort (6 FFTs) one can also run the Forest-Ruth algorithm and obtain results shown as solid triangles, with a fitted line $`9.7\times 10^4(\mathrm{\Delta }t)^4`$. This equal effort comparison clearly demonstrates the greater efficiency of fourth order algorithms. For errors in the range of 0.1 to 0.01 percent, FR’s time step size can be 4 to 6 times as large as SO’s. This ratio would increase if greater accuracy is required. Algorithm 4A, which requires only 2/3 the effort of FR, can be fitted by $`4.8\times 10^6(\mathrm{\Delta }t)^4`$. This fit is plotted as a dotted line barely visible above the zero error line over the range of the plot.

To discuss the convergence of our gradient symplectic algorithms, we greatly expanded the plotting range in Fig.2. Here, we plot directly the convergence of $`E(1000\tau )/E_0`$. The plotting symbols indicate calculated results. We retained SO and FR for comparison. The SO result is well-fitted by $`E_{con}/E_0+0.054(\mathrm{\Delta }t)^2`$. Algorithm 4C and 4D yielded indistinguishable results, despite the fact that 4D uses only 6 FFTs, two less than 4C. Obviously one should not use algorithm 4C in the present case; we only included it here for completeness. Since the FR algorithm is known to have rather large errors, we have also implemented McLachlan’s fourth order algorithm and obtained results labelled as M. This algorithm requires 8 FFTs per iteration and is the best algorithm tested by Sanz-Serna and Portillo. Both FR and M are examples of fourth order algorithms with negative coefficients. The step size convergence of all fourth order algorithms can be very well-fitted by $`E_{con}/E_0+d_i(\mathrm{\Delta }t)^4`$. These are plotted as lines going through data points. The coefficient $`|d_i|`$ for each algorithm is listed in Table I. Since the computational effort per time step is proportional to the number of FFTs, $`N_i`$, the error per unit effort, taking into account the variation of step size with effort, would be proportinal to $`|d_i|N_i^4`$, i.e., this is a measure of error of equal computational effort for all fourth order algorithms. Again, for example, algorithm 4A only requires half the number of of FFTs as McLachlan’s algorithm. At a given $`\mathrm{\Delta }t`$ in running McLachlan’s algorithm, one can execute algorithm 4A twice at $`\mathrm{\Delta }t/2`$ and reduce its error by a factor of $`(1/2)^4`$. Thus despite the appearance in Fig.2, algorithm 4A actually has a much smaller error per unit effort than McLachlan’s algorithm. We normalize this equal effort error to the value of FR and define the normalized, equal-effort error coefficient as $`\delta _{eq}=|d_i/d_{FR}|(N_i/N_{FR})^4`$. This value for each algorithm is listed in Table I. Thus, excluding 4C, our gradient algorithms are roughly a factor of $`10^3`$ smaller in error than FR and a factor of 3 to 8 smaller than McLachlan’s algorithm. A even more useful measure is to consider $`\tau _{eff}=\delta _{eq}^{1/4}`$, which would give the time step size relative to RF for achieving the same error with equal effort. This is also given in Table I. From examining Fig.1 and 2, it certainly seems reasonable that our algorithms can converge at time steps nearly 5 times as large as RF’s and on the order of 30 times as large as SO’s. We will discuss further optimized algorithms 4ACB and 4BDA in the next section.

## V One-Parameter Family of Algorithms

Algorithms 4A and 4C are special cases of the more general algorithm

$$𝒯_{ACB}^{(4)}(ϵ)\mathrm{e}^{t_3ϵT}\mathrm{e}^{v_3ϵV(a_3ϵ)}\mathrm{e}^{t_2ϵT}\mathrm{e}^{v_2ϵ\stackrel{~}{V}(a_2ϵ)}\mathrm{e}^{t_1ϵT}\mathrm{e}^{v_1ϵV(a_1ϵ)}\mathrm{e}^{t_0ϵT}.$$
(36)
By use of a Mathematica program to symbolically combine exponentials of operators, it is easy to determine these positive coefficients to yield a fourth order algorithm. For the above form of the algorithm, there is one free parameter, which we will take to be $`t_0`$. Given $`t_0`$, the rest of the coefficients are:

$$t_1=t_2=\frac{1}{2}t_0,t_3=t_0,v_1=v_3=\frac{1}{6}\frac{1}{(12t_0)^2},v_2=1(v_1+v_3),$$
(37)
where $`\stackrel{~}{V}`$ here is given by

$$\stackrel{~}{V}(t)=V(t)+\frac{u_0}{v_2}ϵ^2[V(t),[T,V(t)]],$$
(38)
with $`a_1=t_0`$, $`a_2=1/2`$, $`a_3=1t_0`$, and

$$u_0=\frac{1}{12}\left[1\frac{1}{12t_0}+\frac{1}{6(12t_0)^3}\right].$$
(39)
For $`t_0=0`$, one has algorithm 4A. For $`t_0=1/60.167`$, one recovers algorithm 4C. Thus one can change continuously from algorithm 4A to 4C, and beyond, by “dialing” the parameter $`0t_0\frac{1}{2}(1\frac{1}{\sqrt{3}})0.21`$. At the upper limit of $`t_0=\frac{1}{2}(1\frac{1}{\sqrt{3}})`$, $`v_2=0`$, and the algorithm becomes a variant of algoirthm 4B, where the commutator term remained at the center rather than equally distributed between the two potential operators at positions $`v_1`$ and $`v_3`$. We shall refer to this algorithm as 4B. Algorithm 4B can be continuously transformed to 4B by redistributing the commutator term from the center to both sides. For example, one can multiply the central commutator term by a factor $`1\alpha `$ and add $`\alpha /2`$ times the commutator term to each potential operator on the side. As $`\alpha `$ ranges from 0 to 1, algorithm 4B is continuously transformed to 4B. When $`\alpha `$ reaches 1, the central commutator disappears and the number of FFTs collapses from 8 to 6. Thus both algorithms 4A and 4B are singular end points of this general algorithm with discontinuous changes in the numbers of FFT.

Algorithm 4B and 4D are special cases of the general algorithm

$$𝒯_{BDA}^{(4)}(ϵ)\mathrm{e}^{v_3ϵV(ϵ)}\mathrm{e}^{t_3ϵT}\mathrm{e}^{v_2ϵ\stackrel{~}{V}(a_2ϵ)}\mathrm{e}^{t_2ϵT}\mathrm{e}^{v_1ϵ\stackrel{~}{V}(a_1ϵ)}\mathrm{e}^{t_1ϵT}\mathrm{e}^{v_0ϵV(0)},$$
(40)
all requiring 6 FFTs. Here, we choose $`t_1`$ as the free parameter. The coefficients are then

$$t_3=t_1,t_2=12t_1,v_3=v_0=\frac{6t_1(t_11)+1}{12(t_11)t_1},v_2=v_1=\frac{1}{2}v_0.$$
(41)
Here $`\stackrel{~}{V}`$ is given by

$$\stackrel{~}{V}(t)=V(t)+\frac{u_0}{v_2}ϵ^2[V(t),[T,V(t)]],$$
(42)
with $`a_1=t_1`$, $`a_2=1t_1`$, and

$$u_0=\frac{1}{48}\left[\frac{1}{6t_1(1t_1)^2}1\right].$$
(43)
Algorithm 4B corresponds to setting $`v_0=0`$ and selecting the smaller of the two quadratic solution, $`t_1=\frac{1}{2}(1\frac{1}{\sqrt{3}})`$. When $`t_1=1/3`$, one obtains a variant of algorithm 4D, which we will denote as 4D. Again, one can transform 4D to 4D continuously by distributing the commutator term at positions $`v_1`$ and $`v_2`$ to $`v_0`$ and $`v_3`$. However, in order to keep the number of gradient terms at a minimum, we will not bother with this refinement here. For positive coefficients, we must have $`\frac{1}{2}(1\frac{1}{\sqrt{3}})t_1\frac{1}{2}`$. At the upper limit of $`t_1=\frac{1}{2}`$, $`t_2=0`$, and the algorithm collapses back to algorithm 4A. Thus, there are two continuous families of algorithms, the C-type requiring eight FFTs and the D-type requiring six. They are joined at both ends by algorithms 4A and 4B with discontinuous changes in the numbers of FFTs. (Recently, Omelyan, Mryglod and Folk have considered the one-parameter factorization schemes from 4A to 4C and 4B to 4D separately, without realizing that they can be further extended and joined into one.)

Fig.2 shows that the fourth order convergence error is negative for algorithm 4A and 4B, and positive for 4C and 4D. Since algorithm 4A can be continuously changed to 4C, and algorithm 4B be transformed to 4D, there must be parameter values such that this fourth order error can be made to vanish. This immediately suggests a simple strategy for further optimization: for each application, one can first minimized the convergence error for a short time run with respect to $`t_0`$ in $`𝒯_{ACB}^{(4)}`$ or $`t_1`$ in $`𝒯_{BDA}^{(4)}`$. One can then use the algorithm at that value and at a large time step, to do long time simulations. Since the error term is generally time dependent, there is no guarantee that when it is small for one time, it will remain small for all time. For the important case of periodic time-dependence, one can hope to reduce the bound within which the error can fluctuate. This seems to work for the present case.

From Fig. 2, as the algorithm is changed continuously from 4A to 4C and back to 4B, one should observe two zero error crossings. In varying $`t_0`$ in $`𝒯_{ACB}^{(4)}`$, we only see one. This suggests that the second crossing is when 4B$`{}_{}{}^{}`$4B, which we did not consider in this work. The one zero crossing we observed is at $`t_0=0.144`$, precisely between 4A($`t_0`$=0) and 4C ($`t_0`$=1/6). A short run ($`t=100\tau `$) to determine the optimal value for $`t_0`$ is shown as solid symbols with solid fitting lines in Fig.3. For $`t_0=0.142`$ and $`t_0=0.146`$, the convergence errors can be well fitted by fourth order lines as shown. At the value of $`t_0=0.144`$, the fourth order coefficient is an order of magnitude smaller. We deem it sufficient to determine the zero-crossing parameter value to three-digit accuracy. The resulting long run is shown as 4ACB in Fig.2. As shown in Table I, $`\delta _{eq}`$ in this case is reduce by an order of magnitude and $`\tau _{eff}`$ lengthened to 8.5 . The convergence error curve for 4ACB remained very flat in going from $`t=100\tau `$ to $`t=1000\tau `$.

In the case of $`𝒯_{BDA}^{(4)}`$, algorithm 4D is sufficiently dissimilar to 4D such that the convergence error remained negative as 4B is morphed to 4D at $`t_1=1/3`$. However, as $`t_1`$ increases above $`1/3`$, the convergence error moves up and turns positive. Since $`𝒯_{BDA}^{(4)}`$ ends up at 4A, the error must cross zero again on its way down. The first crossing is at $`t_1=0.35`$. This is shown in Fig.3 with hollow symbol connected by fitted, broken lines. The second crossing is near $`t_1=0.46`$, but we do not bother to show this similar result. The long-run results using $`t_1=0.35`$ in $`𝒯_{BDA}^{(4)}`$ is shown on Fig.2 as 4BDA. As listed in Table I, its equal effort error and effective time step factor are virtually identical to those of 4ACB. In both cases, without any additional computational effort, we have reduced the error by a factor of five and lengthen the time step size by a factor 1.5 . Both algorithms can achieve 10-4 accuracy in total energy at a time step size 50 times as large as SO.

A key advantage of these factorization algorithms is that one can obtain analytically their energy error terms by use of Mathematica. For fourth order gradient algorithms, there are four non-vanishing error operators in the Hamiltonian,

$$\mathrm{\Delta }E^{(4)}=ϵ^4(e_1[VTVTV]+e_2[TTVTV]+e_3[VTTTV]+e_4[TTTTV]),$$
(44)
where $`e_i`$ are coefficient functions depending on $`t_0`$ or $`t_1`$, and $`[ABC][A,[B,[C]]]`$, etc.. The zero-crossing values of $`t_0`$ in $`𝒯_{ACB}^{(4)}`$ and $`t_1`$ in $`𝒯_{BDA}^{(4)}`$, are closely matched to the crossing points of the first two error coefficients, $`e_1=e_2`$. For 4ACB and 4BDA, the predicted zero-crossing values from solving $`e_1e_2=0`$ numerically are $`t_0=0.14233`$ and $`t_1=0.35023`$ respectively, in excellent agreement with empirical values. It seems that out of the four error operators, only two are dominant and they are opposite in sign. It should be noted that these crossing points have nothing to do with the minimum of $`_{i=1}^4e_i^2`$. These algorithms are not optimized by simply minimizing the sum of squares of error coefficients. Further discussions on the above one-parameter family of gradient algorithms will be presented in a separate publication.

## VI Conclusions

In this work, we have derived a new class of fourth order algorithms for solving the time-dependent Schrödinger equation with explicit time-dependent potentials. This class of algorithms is characterized by having only positive factorization coefficients and can achieve great efficiency by knowing the gradient of the potential. For the Walker-Preston model, on an equal effort basis, the convergence errors of these algorithms can be 5000 and 30 times smaller than negative-coefficient algorithms such as the Forest-Ruth and the McLachlan algorithm, respectively. These gradient algorithms should be further tested on more realistic problems where the gradient of potential may have to be computed numerically.

Our discovery of the one parameter family of gradient algorithms illustrates the power of the operator factorization approach of solving any evolution equation. The ability to tailor a specific algorithm for a particular application reflects great versatility of our method. Our work encourages further systematic study of these and other parametrized families of algorithm and their optimization. In particular, one should explore the freedom in redistributing the commutator term in going from 4B to 4B and from 4D to 4D. We have not done so here in order to concentrate on algorithms with the minimum number of gradient evaluations. With the development of these powerful gradient algorithms for solving the Schrödinger equation directly, we can see no practical advantage in solving the corresponding classical problem using classical methods. This is in agreement with Sanz-Serna and Portillo’s earlier assessment.

This work spurs further interest in finding six and higher order factorization algorithms with purely positive coefficients. Despite intense effort, we have yet to find a sixth order factorization scheme with positive coefficients. It is most likely that one doesn’t exist, even with the inclusion of the double commutator term. A recent work reached the same conclusion, but offered no proof either. Such a proof, if exists, would make these fourth order gradient algorithms unique. There would be no higher order generalizations. This would raise many other interesting questions, such as what other operators are necessary for a higher order positive time step factorization, what would be the optimal sixth order algorithm with minimal of negative coefficients, etc.. All such investigations would deepen our understanding of the operator factorization method of solving evolution equations.

###### Acknowledgements.

This work was supported, in part, by the National Science Foundation grants No. PHY-0100839 to SAC.