# Applying splitting methods with complex coefficients to the numerical
integration of unitary problems

S. Blanes1, F. Casas2, A. Escorihuela-Tomàs3
1 Universitat Politècnica de València, Instituto de Matemática
Multidisciplinar, 46022-Valencia, Spain
email<EMAIL_ADDRESS>
2 Departament de Matemàtiques and IMAC, Universitat Jaume I, 12071-Castellón,
Spain
email<EMAIL_ADDRESS>
3 Departament de Matemàtiques, Universitat Jaume I, 12071-Castellón, Spain
email<EMAIL_ADDRESS>

###### Abstract

We explore the applicability of splitting methods involving complex
coefficients to solve numerically the time-dependent Schrödinger equation. We
prove that a particular class of integrators are conjugate to unitary methods
for sufficiently small step sizes when applied to problems defined in the
group $\mathrm{SU}(2)$. In the general case, the error in both the energy and
the norm of the numerical approximation provided by these methods does not
possess a secular component over long time intervals, when combined with
pseudo-spectral discretization techniques in space.

##  1 Introduction

Splitting methods constitute a natural choice for the numerical time
integration of differential equations of the form

$\frac{du}{dt}=A(u)+B(u),\qquad u(0)=u_{0},$ (1.1)

when each subproblem

$\frac{du}{dt}=A(u),\qquad\frac{du}{dt}=B(u)$

with $u(0)=u_{0}$ can be solved explicitly [3, 15, 19]. Then, by composing the
solution of each part with appropriately chosen coefficients, it is possible
to construct an integrator of a given order $r\geq 1$ for (1.1). In the
particular case of a linear problem,

$\frac{du}{dt}=Au+Bu,$ (1.2)

a splitting method is a composition of the form

$\Psi^{[r]}(h)=\mathrm{e}^{b_{s+1}hB}\mathrm{e}^{a_{s}hA}\mathrm{e}^{b_{s}hB}\cdots\mathrm{e}^{a_{1}hA}\,\mathrm{e}^{b_{1}hB},$
(1.3)

where $h:=\Delta t$ is the time step and the coefficients $a_{j}$, $b_{j}$ are
chosen as solutions of the order conditions, a set of polynomial equations
that must be satisfied to achieve an order of accuracy $r$, i.e., so that
$\exp(h(A+B))u_{0}-\Psi^{[r]}(h)u_{0}=\mathcal{O}(h^{r+1})$.

The simplest example within this class is the Lie–Trotter splitting,

$\mathrm{e}^{hA}\,\mathrm{e}^{hB}\qquad\mbox{ or
}\qquad\mathrm{e}^{hB}\,\mathrm{e}^{hA},$ (1.4)

providing a first order approximation ($r=1$), whereas the palindromic
versions

$\mathcal{S}(h)=\mathrm{e}^{h/2A}\,\mathrm{e}^{hB}\,\mathrm{e}^{h/2A}\qquad\mbox{
or
}\qquad\mathcal{S}(h)=\mathrm{e}^{h/2B}\,\mathrm{e}^{hA}\,\mathrm{e}^{h/2B},$
(1.5)

known as Strang splittings, are methods of order $r=2$.

Although very efficient high order splitting methods can be found in the
literature for the numerical integration of Eq. (1.1), it is important to
remark that if the order $r\geq 3$, then necessarily some of the coefficient
$a_{j}$ and $b_{j}$ have to be negative [2, 21, 22]. This, while does not
constitute a particular problem when the differential equation is reversible,
makes unfeasible their application in parabolic differential equations of
evolutionary type, when the operators $A$ and $B$ are only assumed to generate
$C^{0}$ semi-groups (and not groups): in that case the flows $\mathrm{e}^{tA}$
and/or $\mathrm{e}^{tB}$ may not be defined for $t<0$ [10, 16, 17]. Notice
that this is the case, in particular, if $A$ is the Laplacian operator.

Moreover, even in problems where splitting methods of order $r\geq 3$ can be
safely applied, the presence of negative coefficients usually leads to large
truncation errors, so that more stages than strictly necessary to achieve a
given order have to be included in the composition to reduce these errors and
improve the overall efficiency [6].

It is with the aim of circumventing these drawbacks that splitting methods
with complex coefficients (with positive real part) have entered into the
literature, mainly in the context of the integration of parabolic differential
equations [10, 17, 5], but also for ordinary differential equations (ODEs)
when structure-preserving (symplecticity, energy conservation, reversibility)
is at stake [11].

Splitting and composition methods with complex coefficients, although
computationally between 2 and 4 times more costly than their real counterparts
when applied to ODEs involving real vector fields, possess however some
remarkable properties: their truncation errors with the minimum number of
stages are typically very small, and their stability threshold is
comparatively large. Moreover, when the numerical solution is projected at
each time step, they lead to approximations that still preserve important
qualitative features (such as symplecticity and time-symmetry) up to an order
much higher than the order of the method itself [7, 9, 4].

To better illustrate these points, let us consider a time-symmetric second
order method $\mathcal{S}(h)$ (such as one of the compositions (1.5). Then, a
fourth-order method can be obtained by composition. More specifically, since
the coefficients of such a scheme have to satisfy three order conditions, it
makes sense to take three maps,

$\mathcal{S}(\gamma_{3}h)\,\mathcal{S}(\gamma_{2}h)\,\mathcal{S}(\gamma_{1}h).$

In that case, the order conditions read [6, 15]

$\sum_{j=1}^{3}\gamma_{j}=1,\qquad\sum_{j=1}^{3}\gamma_{j}^{3}=0,\qquad\sum_{j=1}^{2}\left(\gamma_{j}^{3}\left(\sum_{k=j+1}^{3}\gamma_{k}\right)-\gamma_{j}\left(\sum_{k=j+1}^{3}\gamma_{k}^{3}\right)\right)=0.$
(1.6)

and admit only one real solution, namely

$\gamma_{1}=\gamma_{3}=\frac{1}{2-2^{1/3}},\qquad\gamma_{2}=1-2\gamma_{1}$

leading to a time-symmetric composition scheme, usually referred to as
Yoshida’s method, here denoted as $\mathcal{S}^{[4]}(h)$. Notice, however,
that there are four more complex solutions. The first pair,

$\gamma_{1}=\gamma_{3}\equiv\gamma=\frac{1}{2-2^{1/3}\mathrm{e}^{2ik\pi/3}},\qquad\gamma_{2}=1-2\gamma,\qquad
k=1,2$ (1.7)

leads again to two time-symmetric methods, denoted as $\Psi_{P,c}^{[4]}(h)$,
whereas the second one, denoted as $\Psi_{SC,c}^{[4]}(h)$,

$\gamma_{1}=\bar{\gamma}_{3}=\frac{1}{4}\pm
i\,\frac{1}{4}\sqrt{\frac{5}{3}},\qquad\gamma_{2}=\frac{1}{2}$ (1.8)

(here the bar indicates the complex conjugate), corresponds to a so–called
_symmetric-conjugate_ composition method [4]: it is symmetric in the real part
of the coefficients and skew-symmetric in the imaginary part. Here and in the
sequel, the first sub-index in a method (either $P$ or $SC$) refers to its
type (either palindromic or symmetric-conjugate, respectively), whereas the
second sub-index (either $r$ or $c$) indicates that the $a_{i}$ coefficients
in the splitting are real or complex, respectively.

At order five there are two additional order conditions. One of them,
$\omega_{5,1}=\sum_{j=1}^{3}\gamma_{j}^{5}$, has been typically used to
measure the relative error of methods of the same class. If one defines the
error as ${\cal E}=|\omega_{5,1}|$, then one has for the previous methods the
following values of ${\cal E}$:

$\begin{array}[]{ccc}\mathcal{S}^{[4]}(h)&&{\cal E}=5.29\ldots,\\\
\Psi_{P,c}^{[4]}(h)&&{\cal E}=0.024\ldots,\\\ \Psi_{SC,c}^{[4]}(h)&&{\cal
E}=0.027\ldots\end{array}$

Notice that the error of methods with complex coefficients is about 200
smaller than in the real case.

In the particular case in which $\mathcal{S}(h)$ is given by (1.5), the
previous methods can also be written as

$\mathrm{e}^{b_{4}hB}\,\mathrm{e}^{a_{3}hA}\,\mathrm{e}^{b_{3}hB}\,\mathrm{e}^{a_{2}hA}\,\mathrm{e}^{b_{2}hB}\,\mathrm{e}^{a_{1}hA}\,\mathrm{e}^{b_{1}hB},$
(1.9)

with

$b_{1}=\frac{1}{2}\gamma_{1},\ \ \ a_{1}=\gamma_{1},\ \ \
b_{2}=\frac{1}{2}(\gamma_{1}+\gamma_{2}),\ \ \ a_{2}=\gamma_{2},\ \ \
b_{3}=\frac{1}{2}(\gamma_{2}+\gamma_{3}),\ \ \ a_{3}=\gamma_{3},\ \ \
b_{4}=\frac{1}{2}\gamma_{3}.$

As a matter of fact, the simplest symmetric-conjugate composition corresponds
to the third order scheme

$\Psi_{SC,c}^{[3]}(h)=\mathcal{S}_{\alpha_{2}h}^{[2]}\circ\mathcal{S}_{\alpha_{1}h}^{[2]},$
(1.10)

with

$\alpha_{1}=\bar{\alpha}_{2}\equiv\alpha=\frac{1}{2}+i\frac{\sqrt{3}}{6}.$

For equation (1.2), method (1.10) can be written as

$\Psi_{SC,c}^{[3]}(h)=\mathrm{e}^{\bar{b}_{1}hB}\,\mathrm{e}^{\bar{a}_{1}hA}\,\mathrm{e}^{b_{2}hB}\,\mathrm{e}^{a_{1}hA}\,\mathrm{e}^{b_{1}hB},$
(1.11)

with $a_{1}=\alpha$, $b_{1}=\alpha/2$, $b_{2}=1/2$.

Although (1.11) is of order 3, if $A$ and $B$ are real, then it renders a
scheme of order 4 when it is projected on the real axis after each time step.
In addition, it verifies
$\Psi_{SC,c}^{[3]}(-h)\circ\Psi_{SC,c}^{[3]}(h)=I+\mathcal{O}(h^{8})$. It is
said that the scheme is pseudo-symmetric of order 7, since it preserves the
time-symmetry property up to terms of order $h^{7}$ [9].

Schemes with complex coefficients have been proposed before for the treatment
of quantum problems, mainly in the context of imaginary time propagation, with
the purpose of computing ground state energies [1] and in quantum Monte Carlo
simulations [23, 14], but also in the decomposition of unitary operators [20].
In the later case it is shown, both for unitary $2\times 2$ matrices and
empirically for exponentials of Gaussian random Hermitian matrices, that a
splitting method does indeed possess a maximal time step for which the scheme
is numerically stable. We generalize the treatment to differential equations
defined in $\mathrm{SU}(2)$ for methods possessing a particular symmetry and
eventually examine their behavior when they are applied to the time dependent
Schrödinger equation.

##  2 Splitting methods in $\mathrm{SU}(2)$

In the study of the evolution of two-level quantum systems one has to deal
with the Schrödinger equation, which in this context reads ($\hbar=1$)

$i\frac{dU}{dt}=H\,U,\qquad U(0)=I,$ (2.1)

where $U(t)$ is a $2\times 2$ unitary matrix with unit determinant and the
skew-Hermitian Hamiltonian $H$ can be expressed as a linear combination of
Pauli matrices,

$\sigma_{1}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right),\qquad\sigma_{2}=\left(\begin{array}[]{cr}0&-i\\\
i&0\end{array}\right),\qquad\sigma_{3}=\left(\begin{array}[]{cr}1&0\\\
0&-1\end{array}\right).$ (2.2)

Since our purpose is to analyze splitting methods in this context, we assume
that $H$ can be written as

$H=\mathbf{a}\cdot\mathbf{\sigma}+\mathbf{b}\cdot\mathbf{\sigma}$ (2.3)

for given vectors $\mathbf{a},\mathbf{b}\in\mathbb{R}^{3}$ and
$\mathbf{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})$, so that, by comparing
with (1.2), one has $A\equiv-i\,\mathbf{a}\cdot\mathbf{\sigma}$ and
$B\equiv-i\,\mathbf{b}\cdot\mathbf{\sigma}$. The exact solution of Eq. (2.1)
after one time step $h$ is

$U_{\mathrm{ex}}(h)=\mathrm{e}^{-ihH}=\mathrm{e}^{-ih(\mathbf{a}\cdot\mathbf{\sigma}+\mathbf{b}\cdot\mathbf{\sigma})}.$

On the other hand, if a splitting method of the form (1.3) of order $r$ with
_real_ coefficients is applied to solve this very simple problem, it is clear
that the corresponding approximation can be written as

$U_{\mathrm{app}}(h)=\mathrm{e}^{-i\,h\,\mathbf{d}(h)\cdot\mathbf{\sigma}},\qquad\mbox{
where
}\qquad\mathbf{d}(h)=\mathbf{a}+\mathbf{b}+\mathcal{O}(h^{r})\in\mathbb{R}^{3},$

and thus the method still renders an approximation in $\mathrm{SU}(2)$. The
situation is different, however, when the splitting method (1.3) involves
_complex_ coefficients, since in that case $\mathbf{d}(h)\in\mathbb{C}^{3}$
and the approximation is no longer unitary. In general, the scheme will be
unstable and the errors will grow exponentially along the integration.

#### Example.

At this point it is worth testing the previous third- and fourth-order schemes
obtained by composing the Strang splitting and involving complex coefficients,
namely $\Psi_{SC,c}^{[3]}(h)$, $\Psi_{SC,c}^{[4]}(h)$, and
$\Psi_{P,c}^{[4]}(h)$. To do that, we consider the following simple
Hamiltonian in $\mathrm{SU}(2)$: $H=\sigma_{1}+\sigma_{2}$, or alternatively,
$\mathbf{a}=(1,0,0)$, $\mathbf{b}=(0,1,0)$ in (3.5).

In our experiment, we take as initial condition $U(t_{0}=0)=I$, integrate Eq.
(2.1) with different values of the time step $h$ and compute the error of the
approximation (in the 2-norm) at the final time $t_{f}=10$ as a function of
the computational cost (estimated as the number of exponentials involved in
the whole integration). The results obtained with each method are displayed in
Figure 1 (left panel). We notice that all schemes involving complex
coefficients provide considerably more accurate results than
$\mathcal{S}^{[4]}(h)$ (black solid line), the fourth-order methods being also
more efficient than $\Psi_{SC,c}^{[3]}(h)$ for high accuracy.

In order to check how each scheme with complex coefficients behaves with
respect to unitarity, we take as a final time $t_{f}=1000$, and adjust $h$
(and therefore the number of iterations $N$) so that they require the same
computational cost. Specifically, $N=6000$ ($h=1/6$) for scheme
$\Psi_{SC,c}^{[3]}(h)$ and $N=4000$ ($h=1/4$) for all methods of order 4.
Finally we compute $|\|U_{\mathrm{app}}(nh)\|-1|,\
n=1,2,\ldots,N=(t_{f}-t_{0})/h$, where $U_{\mathrm{app}}(nh)$ denotes the
approximate solution after $n$ steps. The outcome is depicted in Figure 1
(right panel). Notice how the error in unitarity grows for
$\Psi_{P,c}^{[4]}(h)$, whereas it is bounded, even for large intervals, for
the symmetric-conjugate methods $\Psi_{SC,c}^{[3]}(h)$ and
$\Psi_{SC,c}^{[4]}(h)$. Among them, the later clearly provides more accurate
results.

Figure 1:  Left: 2-norm error vs. computational cost (number of exponentials)
for $\Psi_{SC,c}^{[3]}$ (dotted line), $\mathcal{S}^{[4]}$ (real coefficients,
solid line), $\Psi_{P,c}^{[4]}$ (complex coefficients, dash-dotted line) and
$\Psi_{SC,c}^{[4]}$ (dashed line) . Right: Error in unitarity for
$\Psi_{P,c}^{[4]}$ (solid line) and the symmetric-conjugate methods
$\Psi_{SC,c}^{[3]}$ (dotted line) and $\Psi_{SC,c}^{[4]}$ (dashed line).

This marked difference of both types of integrators can also be illustrated by
computing the eigenvalues $\lambda_{1}$, $\lambda_{2}$ of the approximate
solution after one step size, i.e., of the corresponding matrix
$U_{\mathrm{app}}(h)$. In the exact case, of course, both evolve on the unit
circle in the complex plane, whereas here one has still
$\lambda_{1}\lambda_{2}=1$, since the determinant is one. In Figure 2 we
depict $|\lambda_{j}|$, $j=1,2$, as a function of $h$ for the palindromic
scheme $\Psi_{P,c}^{[4]}$ with $k=1$ (black, dashed lines) and the symmetric-
conjugate splittings $\Psi_{SC,c}^{[3]}$ (blue, dotted line) and
$\Psi_{SC,c}^{[4]}$ (red, solid line) in the range $1\leq h\leq 3$. It is
worth remarking that for the symmetric-conjugate methods both $|\lambda_{j}|$
are exactly 1 for $0\leq h\leq h^{*}$, with $h^{*}=1.7570473$ for
$\Psi_{SC,c}^{[3]}$ and $h^{*}=2.9139468357$ for $\Psi_{SC,c}^{[4]}$. In other
words, they behave as unitary maps when $h\leq h^{*}$. On the other hand, it
can be checked that $|\lambda_{1}|>1$ for any $h>0$ for $\Psi_{P,c}^{[4]}$.
$\Box$

Figure 2:  Absolute value of the eigenvalues of the approximate solution
matrix obtained with $\Psi_{P,c}^{[4]}$ with complex coefficients ($k=1$,
black dashed line), $\Psi_{SC,c}^{[3]}$ (blue dotted line) and
$\Psi_{SC,c}^{[4]}$ (red, solid line).

The previous example illustrates in fact a general pattern exhibited by
symmetric-conjugate methods for this problem, as we next prove.

###### Proposition 2.1

Suppose a symmetric-conjugate splitting method of the form (1.3), with
$a_{s+1-j}=\bar{a}_{j}$, $b_{s+2-j}=\bar{b}_{j}$, is applied to the numerical
integration of the Schrödinger equation (2.1) with the Hamiltonian given by
(3.5). In that case, the following statements hold:

  * (a)

The eigenvalues of the matrix approximating the solution after on time step
$h$ lie on the unit circle in the complex plane for sufficiently small $h$.

  * (b)

The symmetric-conjugate splitting method is itself conjugate to a unitary
method for sufficiently small $h$.

Proof.  When a splitting method of the form (1.3) is applied to solve Eq.
(2.1), the corresponding approximation after one step can be written as
$U_{\mathrm{app}}(h)=\exp(V(h))$, where $V(h)$ is a linear combination of $A$,
$B$ and all their nested commutators,

$V(h)=h(w_{1,1}A+w_{1,2}B)+h^{2}w_{2,1}[A,B]+h^{3}(w_{3,1}[A,[A,B]]+w_{3,2}[B,[A,B]])+\mathcal{O}(h^{4}),$
(2.4)

and $w_{n,k}$ are polynomials in the coefficients $a_{j}$, $b_{j}$. Method
(1.3) is of order $r$ if $w_{1,1}=w_{1,2}=1$ and the polynomials $w_{n,k}$
vanish for $1<n<r$. In our case, since

$[-i\,\mathbf{a}\cdot\mathbf{\sigma},-i\,\mathbf{b}\cdot\mathbf{\sigma}]=-i\,2\,(\mathbf{a}\times\mathbf{b})\cdot\mathbf{\sigma},$

it is clear that the vector fields associated with all commutators in (2.4)
containing an even number of operators are perpendicular to the plane
generated by the vectors $\mathbf{a}$ and $\mathbf{b}$, whereas those
containing an odd number of operators $A$ and $B$ are in such a plane. If in
addition the method is symmetric-conjugate, then a straightforward computation
shows that $V^{{\dagger}}(h)=V(-h)$, and so all polynomials $w_{2j+1,k}$ are
real whereas all polynomials $w_{2j,k}$ are pure imaginary. Therefore, $V(h)$
can be written as

$V(h)=-i\,h\tilde{H}(h),\qquad\mbox{ with
}\qquad\tilde{H}(h)=\mathbf{d}(h)\cdot\mathbf{\sigma}+i\,\mathbf{c}(h)\cdot\mathbf{\sigma},$
(2.5)

for two vectors $\mathbf{c},\mathbf{d}\in\mathbb{R}^{3}$ verifying
$\mathbf{c}\cdot\mathbf{d}=0$ and

$\mathbf{d}(h)=\mathbf{a}+\mathbf{b}+\mathcal{O}(h^{r+1}),\qquad\qquad\mathbf{c}(h)=\mathcal{O}(h^{r}).$
(2.6)

This special structure of $\tilde{H}(h)$ allows one to obtain statements (a)
and (b) above. First, if we write

$\mathrm{e}^{-i\,h\tilde{H}(h)}=\mathrm{e}^{h\,\mathbf{u}\cdot\mathbf{\sigma}},\qquad\mbox{
with }\qquad\mathbf{u}=\mathbf{c}-i\,\mathbf{d},$

then

$U_{\mathrm{app}}(h)=\mathrm{e}^{h\,\mathbf{u}\cdot\mathbf{\sigma}}=\cosh(hu)I+\frac{\sinh(hu)}{u}\,\mathbf{u}\cdot\mathbf{\sigma},$

with
$u=\sqrt{\mathbf{u}\cdot\mathbf{u}}=(\|\mathbf{c}\|^{2}-\|\mathbf{d}\|^{2})^{1/2}$.
Of course, if $\|\mathbf{c}\|<\|\mathbf{d}\|$, then $\cosh(u)=\cos\alpha$,
$\sinh(u)=i\sin\alpha$, with
$\alpha=(\|\mathbf{d}\|^{2}-\|\mathbf{c}\|^{2})^{1/2}$, and the eigenvalues of
$U_{\mathrm{app}}(h)$ are $\lambda_{1,2}=\exp(\pm ih\alpha(h))$. But, in
virtue of (2.6), this always holds for sufficiently small values of $h$.

Statement (b) can demonstrated as follows. Let us introduce the unitary vector

$\mathbf{C}=\frac{\mathbf{d}\times\mathbf{c}}{\|\mathbf{d}\times\mathbf{c}\|}.$

A trivial computation shows that

$\mathbf{d}\times\mathbf{C}=-\frac{\|\mathbf{d}\|}{\|\mathbf{c}\|}\,\mathbf{c},\qquad\qquad\mathbf{c}\times\mathbf{C}=\frac{\|\mathbf{c}\|}{\|\mathbf{d}\|}\,\mathbf{d},$

and furthermore, for a given parameter $s\in\mathbb{R}$,

$\mathrm{e}^{s\,\mathbf{C}\cdot\mathbf{\sigma}}\,\mathrm{e}^{h\mathbf{u}\cdot\mathbf{\sigma}}\,\mathrm{e}^{-s\,\mathbf{C}\cdot\mathbf{\sigma}}=\exp\left(\mathrm{e}^{s\,\mathbf{C}\cdot\mathbf{\sigma}}\,(h\,\mathbf{u}\cdot\mathbf{\sigma})\,\mathrm{e}^{-s\,\mathbf{C}\cdot\mathbf{\sigma}}\right).$

From the definition of $\mathbf{C}$ and the properties of the Pauli matrices
[13], one has

$\displaystyle\mathrm{e}^{s\,\mathbf{C}\cdot\mathbf{\sigma}}\,(\mathbf{c}\cdot\mathbf{\sigma})\,\mathrm{e}^{-s\,\mathbf{C}\cdot\mathbf{\sigma}}$
$\displaystyle=$
$\displaystyle\cosh(2s)(\mathbf{c}\cdot\mathbf{\sigma})-i\,\sinh(2s)\frac{\|\mathbf{c}\|}{\|\mathbf{d}\|}\,\mathbf{d}\cdot\mathbf{\sigma}$
$\displaystyle\mathrm{e}^{s\,\mathbf{C}\cdot\mathbf{\sigma}}\,(\mathbf{d}\cdot\mathbf{\sigma})\,\mathrm{e}^{-s\,\mathbf{C}\cdot\mathbf{\sigma}}$
$\displaystyle=$
$\displaystyle\cosh(2s)(\mathbf{d}\cdot\mathbf{\sigma})+i\,\sinh(2s)\frac{\|\mathbf{d}\|}{\|\mathbf{c}\|}\,\mathbf{c}\cdot\mathbf{\sigma},$

and thus

$\mathrm{e}^{s\,\mathbf{C}\cdot\mathbf{\sigma}}\,(\mathbf{u}\cdot\mathbf{\sigma})\,\mathrm{e}^{-s\,\mathbf{C}\cdot\mathbf{\sigma}}=\left(\sinh(2s)\frac{\|\mathbf{d}\|}{\|\mathbf{c}\|}+\cosh(2s)\right)\mathbf{c}\cdot\mathbf{\sigma}-i\left(\sinh(2s)\frac{\|\mathbf{c}\|}{\|\mathbf{d}\|}+\cosh(2s)\right)\mathbf{d}\cdot\mathbf{\sigma}.$

If we now take $s$ such that

$\sinh(2s)\frac{\|\mathbf{d}\|}{\|\mathbf{c}\|}+\cosh(2s)=0,\qquad\mbox{ i.e.,
}\qquad\tanh(2s)=-\frac{\|\mathbf{c}\|}{\|\mathbf{d}\|},$ (2.7)

then, clearly

$\mathrm{e}^{s\,\mathbf{C}\cdot\mathbf{\sigma}}\,(-ih\tilde{H})\,\mathrm{e}^{-s\,\mathbf{C}\cdot\mathbf{\sigma}}=-i\,h\,\mathbf{D}(h)\cdot\mathbf{\sigma},$

with

$\mathbf{D}(h)=\sinh(2s)\left(\frac{\|\mathbf{c}\|}{\|\mathbf{d}\|}-\frac{\|\mathbf{d}\|}{\|\mathbf{c}\|}\right)\mathbf{d}\in\mathbb{R}^{3}.$

In consequence,

$\mathrm{e}^{s\,\mathbf{C}\cdot\mathbf{\sigma}}\,\mathrm{e}^{-ih\tilde{H}}\,\mathrm{e}^{-s\,\mathbf{C}\cdot\mathbf{\sigma}}=\mathrm{e}^{-ih\mathbf{D}(h)\cdot\mathbf{\sigma}}.$
(2.8)

In other words, if $s$ is such that Eq. (2.7) holds, then the map obtained by
applying a symmetric-conjugate splitting method is conjugate to a unitary
matrix. Notice that if $\|\mathbf{c}\|<\|\mathbf{d}\|$ this is always
possible, in agreement with statement (a) for the eigenvalues of the
approximate solution matrix.

Proposition 2.1 thus provides a rigorous justification of the results shown in
Figures 1 and 2: since a symmetric-conjugate splitting method is ultimately
conjugate to a unitary map in the sense of eq. (2.8) for sufficiently small
values of $h$, then the error in the unitarity of the numerical solution is
bounded, whereas the eigenvalues remain on the unit circle in the complex
plane.

Methods of the form (2.8) are called _processed_ or _corrected_ in the
literature (see, e.g. [3, 6, 15, 19]). In that context, method
$\mathrm{e}^{-ih\tilde{H}}$ is called the kernel, and
$\mathrm{e}^{s\,\mathbf{C}\cdot\mathbf{\sigma}}$ the processor. For
integrators of this class, only the error terms in the kernel that cannot be
removed by a processor are relevant in the long run. In the case of unitary
problems in $\mathrm{SU}(2)$ we have shown that any symmetric-conjugate
splitting method is indeed the kernel of a processed unitary scheme.

##  3 Application to the time-dependent Schrödinger equation

In view of the previous results in $\mathrm{SU}(2)$, it is natural to examine
the situation when splitting methods with complex coefficients, and in
particular symmetric-conjugate schemes, are applied in a more general setting.
To this end, we next consider the numerical integration of the general time
dependent Schrödinger equation

$i\frac{\partial}{\partial
t}\psi(x,t)=-\frac{1}{2\mu}\Delta\psi(x,t)+V(x)\psi(x,t),$ (3.1)

where now $\psi:\mathbb{R}^{d}\times\mathbb{R}\longrightarrow\mathbb{C}$ is
the wave function representing the state of the system and the initial state
is $\psi(x,0)=\psi_{0}(x)$. We take again $\hbar=1$ and a reduced mass $\mu$.
Equation (3.1) can be written as

$i\frac{\partial}{\partial t}\psi=(\hat{T}(P)+\hat{V}(X))\psi,$ (3.2)

with $\hat{T}(P)=P^{2}/(2\mu)$, and the operators $X$ and $P$ are defined by
their actions on $\psi(x,t)$ as

$X\psi(x,t)=x\,\psi(x,t),\qquad\quad P\ \psi(x,t)=-i\,\nabla\psi(x,t).$ (3.3)

The usual procedure for applying splitting methods in this setting consists
first in discretizing the space variables $x$, so as to get a system of
ordinary differential equations (ODEs) which is subsequently integrated in
time by the splitting scheme. If, for simplicity, we consider the one-
dimensional problem, $d=1$, and suppose that it is defined in
$x\in[x_{0},x_{N}]$, first this interval is partitioned into $N$ parts of
length $\Delta x=(x_{N}-x_{0})/N$ and the vector
$u=(u_{0},\ldots,u_{N-1})^{T}\in\mathbb{C}^{N}$ is formed, with
$u_{n}=\psi(x_{n},t)$ and $x_{n}=x_{0}+n\Delta x$, $n=0,1,\ldots,N-1$. The
partial differential equation (3.1) is then replaced by the $N$-dimensional
linear ODE

$i\frac{d}{dt}u(t)=H\,u(t),\qquad u(0)=u_{0}\in\mathbb{C}^{N},$ (3.4)

where now $H$ represents the (real symmetric) $N\times N$ matrix associated
with the Hamiltonian.

When a Fourier spectral collocation method is used, then the matrix $H$ in
(3.4) is

$H=T+V,$ (3.5)

where $V$ is a diagonal matrix associated with the potential $\hat{V}$ and $T$
is a (full) differentiation matrix related with the kinetic energy $\hat{T}$.
Their action on the wave function vector $u$ is trivial: on the one hand,
$(Vu)_{n}=V(x_{n})u_{n}$ and thus the product $Vu$ requires to compute $N$
complex multiplications. On the other hand,
$Tu=\mathcal{F}^{-1}D_{T}\mathcal{F}u$, where $\mathcal{F}$ and
$\mathcal{F}^{-1}$ are the forward and backward discrete Fourier transform,
and $D_{T}$ is again diagonal. The transformation $\mathcal{F}$ from the
discrete coordinate representation to the discrete momentum representation
(and back) is done via the fast Fourier transform (FFT) algorithm, requiring
$\mathcal{O}(N\log N)$ operations.

Notice that, since

$\left(\mathrm{e}^{\tau V}u\right)_{i}=\mathrm{e}^{\tau
V(x_{i})}u_{i},\qquad\quad\mathrm{e}^{\tau
T}u=\mathcal{F}^{-1}\mathrm{e}^{\tau D_{T}}\mathcal{F}\,u,$

splitting methods constitute a valid alternative to approximate the solution
$u(t)=\mathrm{e}^{\tau(T+V)}u_{0}$ for a time step $\Delta t$, with
$\tau=-i\Delta t$. Thus, with the Lie–Trotter scheme (1.4) one has

$\mathrm{e}^{\tau(T+V)}=\mathrm{e}^{\tau T}\,\mathrm{e}^{\tau
V}+\mathcal{O}(\tau^{2}),$

whereas the 2nd-order Strang splitting (1.5) constructs the numerical
approximation $u_{n+1}$ at time $t_{n+1}=t_{n}+\Delta t$ by

$u_{n+1}=\mathrm{e}^{\tau/2V}\,\mathrm{e}^{\tau
T}\,\mathrm{e}^{\tau/2V}\,u_{n}\equiv\mathcal{S}(\tau)u_{n}.$

The resulting scheme is called split-step Fourier method in the chemical
literature, and has some remarkable properties. In particular, it is both
unitary and symplectic [18], as well as time-reversible. In addition, for
suitable regularity assumptions on the potential and on the norm of the
commutators $[\hat{T},\hat{V}]=\hat{T}\hat{V}-\hat{V}\hat{T}$ and
$[\hat{T},[\hat{T},\hat{V}]]$, the error at $t_{n}$ is bounded by

$\|u_{n}-u(t)\|\leq C\,\Delta t^{2}\,t\,\max_{0\leq s\leq t}\|u(s)\|_{2}.$

Higher order methods can be obtained of course by considering compositions
(1.3), which in this setting read

$\Psi^{[r]}(\tau)=\mathrm{e}^{b_{s+1}\tau V}\mathrm{e}^{a_{s}\tau
T}\mathrm{e}^{b_{s}\tau V}\cdots\mathrm{e}^{a_{1}\tau
T}\,\mathrm{e}^{b_{1}\tau V},$ (3.6)

and in fact, a large collection of practical schemes of different orders exist
for carrying out the numerical integration (see e.g. [6, 15, 19] and
references therein). In addition, from (3.3), it is clear that
$[\hat{X},\hat{P}]\ \psi(x,t)=i\,\psi(x,t),$ and so

$[\hat{V},[\hat{V},[\hat{V},\hat{T}]]]\ \psi(x,t)=0.$ (3.7)

This property leads to a reduction in the number of order conditions necessary
to achieve a given order $r$ and allows one to construct highly efficient
schemes.

##  4 Splitting methods with complex coefficients

When exploring the applicability of splitting methods with complex
coefficients to the general time-dependent Schrödinger equation, several
aspects must be addressed. First, since the computational cost of method (3.6)
is dominated by the number of FFTs per step, the presence of complex $a_{j}$,
$b_{j}$ does not contribute significantly to increase this cost. In addition,
it has been shown in other problems that splitting methods with complex
coefficients involving the minimum number of flows to achieve a given order
already provide good efficiency, in contrast with their real counterparts. On
the other hand, however, since $\sum_{j}a_{j}=1$ for a consistent method, if
$a_{j}\in\mathbb{C}$, then imaginary parts positive _and_ negative enter into
the game, with the result that severe instabilities may arise in practice due
to the unboundedness of the Laplace operator. With respect to the potential,
since in regions where it takes large values the wave function typically is
close to zero, we can introduce an artificial cut-off bound in the computation
if necessary, so that complex $b_{j}$ can in principle be used, at least for a
sufficiently small $\Delta t$. It makes sense therefore to construct and
examine in detail methods with real $a_{j}$ and complex $b_{j}$ coefficients.

In the following, and for simplicity, we restrict ourselves to splitting
methods (3.6) of order $r\leq 4$, which we denote by their sequence of
coefficients as

$(b_{s+1},a_{s},b_{s},\ldots,a_{2},b_{2},a_{1},b_{1}).$

The order conditions are then

Order 1:
$\displaystyle\sum_{i=1}^{s}a_{i}=1,\qquad\qquad\qquad\qquad\sum_{i=1}^{s}b_{i}=1,$
Order 2:
$\displaystyle\sum_{i=1}^{s}b_{i}\left(\sum_{j=1}^{i}a_{j}\right)=\frac{1}{2},$
Order 3:
$\displaystyle\sum_{i=1}^{s}b_{i}\left(\sum_{j=1}^{i}a_{j}\right)^{2}=\frac{1}{3},\qquad\quad\sum_{i=1}^{s}a_{i}\left(\sum_{j=i}^{s}b_{j}\right)^{2}=\frac{1}{3},$
(4.1) Order 4:
$\displaystyle\sum_{i=1}^{s}b_{i}\left(\sum_{j=1}^{i}a_{j}\right)^{3}=\frac{1}{4},\qquad\quad\sum_{i=1}^{s}a_{i}\left(\sum_{j=i}^{s}b_{j}\right)^{3}=\frac{1}{4},$
$\displaystyle\sum_{i=1}^{s}a_{i}^{2}\left(\sum_{j=i}^{s}b_{j}\right)^{2}+2\,\sum_{i=2}^{s}a_{i}\left(\sum_{j=1}^{i-1}a_{j}\right)\left(\sum_{k=i}^{s}b_{k}\right)^{2}=\frac{1}{6}.$

In typical applications of splitting methods with real coefficients, only
palindromic sequences of coefficients, i.e., methods (3.6) with
$b_{s+2-j}=b_{j}$, $a_{s+1-j}=a_{j}$ for all $j$ are used. In that case, all
the conditions at even order are automatically satisfied and the resulting
schemes are time-symmetric, $\Psi^{[r]}(\tau)\,\Psi^{[r]}(-\tau)=I$, and of
even order. Here, however, since we are dealing with complex coefficients, we
also analyze the case $r=3$ for completeness.

#### Order 3.

The first five order conditions in (4) admit solutions with all $a_{j}$ real
and positive and $b_{j}\in\mathbb{C}$ with positive real part if one considers
a composition of the form

$\Psi_{SC,r}^{[3]}(\tau)=(\bar{b}_{1},a_{1},\bar{b}_{2},a_{2},b_{2},a_{1},b_{1})$
(4.2)

involving 6 parameters. Then one gets a 1-parametric family of solutions
(+c.c.) with the required properties. Among them, we choose

$a_{1}=\frac{3}{10},\qquad a_{2}=\frac{2}{5},\qquad
b_{1}=\frac{13}{126}-i\frac{\sqrt{59/2}}{63},\qquad
b_{2}=\frac{25}{63}+i\frac{5\sqrt{59/2}}{126}.$

Composition (1.11) constitutes of course another symmetric-conjugate method of
order 3, denoted here by

$\Psi_{SC,c}^{[3]}(\tau)=(\bar{b}_{1},\bar{a}_{1},b_{2},a_{1},b_{1})$ (4.3)

and involving less maps, although in this case $a_{1}\in\mathbb{C}$.

#### Order 4.

The simplest approach to construct a palindromic scheme with
$a_{j}\in\mathbb{R}$ and $b_{j}\in\mathbb{C}$ consists in taking all the
$a_{j}$ equal. In that case, with $s=4$, one has enough parameters to solve
the required four order conditions (at odd orders). Only two solutions
(complex conjugate to each other) are obtained, as shown in [10], thus
resulting in the scheme

$\Psi_{P,r}^{[4]}(\tau)=(b_{1},a_{1},b_{2},a_{2},b_{3},a_{2},b_{2},a_{1},b_{1})$
(4.4)

with

$a_{1}=a_{2}=\frac{1}{4},\qquad b_{1}=\frac{1}{10}-i\frac{1}{30},\qquad
b_{2}=\frac{4}{15}+i\frac{2}{15},\qquad b_{3}=\frac{4}{15}-i\frac{1}{5}.$

Although more efficient schemes can be obtained if one allows for different
$a_{j}$’s [5], since we are interested here mainly in the qualitative behavior
of the different methods, we limit ourselves to (4.4) as representative of
palindromic splitting methods with real $a_{j}$’s and complex $b_{j}$’s,
whereas we can take scheme (1.9)

$\Psi_{P,c}^{[4]}(\tau)=(b_{1},a_{1},b_{2},a_{2},b_{2},a_{1},b_{1}),$ (4.5)

as representative of palindromic methods with both $a_{j}\in\mathbb{C}$ and
$b_{j}\in\mathbb{C}$.

Symmetric-conjugate splitting methods with real $a_{j}$’s require at least
$s=5$ stages, in which case one has a free parameter. If we fix this as
$a_{1}=1/8$, we get the scheme

$\Psi_{SC,r}^{[4]}(\tau)=(\bar{b}_{1},a_{1},\bar{b}_{2},a_{2},\bar{b}_{3},a_{3},b_{3},a_{2},b_{2},a_{1},b_{1})$
(4.6)

with

$\displaystyle a_{2}=0.23670501659941197298,$ $\displaystyle
a_{3}=0.27658996680117605403,$ $\displaystyle
b_{1}=0.03881396214419327198-0.045572109263923104872\,i,$ $\displaystyle
b_{2}=0.19047619047619047619+0.115462072300408741306\,i,$ $\displaystyle
b_{3}=0.27070984737961625182-0.148322245509626403888\,i$

It is worth noticing that one can obtain symmetric-conjugate methods from
palindromic schemes and vice versa. Thus, in particular, by composing the
palindromic scheme (4.4) with its complex conjugate we can form a symmetric-
conjugate splitting method with 8 stages and $a_{j}\in\mathbb{R}$,
$b_{j}\in\mathbb{C}$:

$\Xi_{SC,r}^{[4]}(\tau)=\Psi_{P,r}^{[4]}(\tau/2)\,\overline{\Psi}_{P,r}^{[4]}(\tau/2),$
(4.7)

whereas doing the same with the 3rd-order symmetric-conjugate method (4.2)
results in the 4th-order palindromic scheme with 6 stages and
$a_{j}\in\mathbb{R}$, $b_{j}\in\mathbb{C}$:

$\Xi_{P,r}^{[4]}(\tau)=\Psi_{SC,r}^{[3]}(\tau/2)\,\overline{\Psi}_{SC,r}^{[3]}(\tau/2).$
(4.8)

This is possible because the adjoint of $(\Psi_{SC,r}^{[3]}(\tau))^{*}$
verifies

$(\Psi_{SC,r}^{[3]}(\tau))^{*}=\overline{\Psi}_{SC,r}^{[3]}(\tau).$

In our numerical experiments we shall also use for comparison one of the best
4th-order splitting methods with real coefficients designed specifically for
systems verifying (3.7). It reads

$\Psi_{RKN}^{[4]}(\tau)=(b_{1},a_{1},b_{2},a_{2},b_{3},a_{3},b_{4},a_{3},b_{3},a_{2},b_{2},a_{1},b_{1})$
(4.9)

and the coefficients can be found in [8]. The scheme has three additional
parameters that are used to minimize error terms at higher orders, and
provides by construction unitary approximations.

##  5 Numerical experiments

We next report on some numerical tests we have carried out with the splitting
methods presented in section 4 applied to the Scrödinger equation in one
dimension. Since many different schemes are tested and compared, it is
convenient to classify them into the following categories:

  * •

symmetric-conjugate methods with $a_{j}\in\mathbb{R}$, $b_{j}\in\mathbb{C}$

    * –

Order 3: $\Psi_{SC,r}^{[3]}$, Eq. (4.2);

    * –

Order 4: $\Psi_{SC,r}^{[4]}$, Eq. (4.6);

  * •

symmetric-conjugate with $a_{j}\in\mathbb{C}$, $b_{j}\in\mathbb{C}$: method
$\Psi_{SC,c}^{[3]}$, Eq. (4.3), order 3;

  * •

palindromic with $a_{j}\in\mathbb{R}$, $b_{j}\in\mathbb{C}$: method
$\Psi_{P,r}^{[4]}$, Eq. (4.4), order 4;

  * •

palindromic with $a_{j}\in\mathbb{C}$, $b_{j}\in\mathbb{C}$: method
$\Psi_{P,c}^{[4]}$, Eq. (4.5), order 4;

For completeness, we also consider the following schemes of order 4 with
$a_{j}\in\mathbb{R}$, $b_{j}\in\mathbb{C}$:

  * •

symmetric-conjugate obtained from a palindromic method: $\Xi_{SC,r}^{[4]}$,
Eq. (4.7);

  * •

palindromic obtained from a symmetric-conjugate method: $\Xi_{P,r}^{[4]}$, Eq.
(4.8).

#### Quartic potential.

As the first example we take the quartic oscillator

$V(x)=-\frac{1}{2}x^{2}+\frac{1}{20}x^{4}$ (5.1)

and the initial condition $\psi_{0}(x)=\sigma\,\mathrm{e}^{-x^{2}/2}$, where
$\sigma$ is a normalization constant. As usual, and since the exact solution
decays rapidly, we truncate the infinite spatial domain to the periodic domain
$[-L,L]$, provided $L$ is sufficiently large and use Fourier spectral methods.
We take $L=8$ and set up a uniform grid on the interval with $N=128$
subdivisions. Finally, we apply the different schemes to integrate in time the
resulting equation (3.4) in the interval $t\in[0,t_{f}]$, with $t_{f}=8000$.
As in the case of the example in $\mathrm{SU}(2)$, we check the behavior of
each scheme with respect to unitarity by computing
$|\|u_{\mathrm{app}}(t)\|-1|$ along the integration, where
$u_{\mathrm{app}}(t)$ denotes the numerical approximation obtained by each
method.

In addition, we also compute the expected value of the energy,
$u_{\mathrm{app}}^{*}(t)\cdot Hu_{\mathrm{app}}(t)$ and measure the error as
the difference with respect to the exact value:

$\mbox{energy
error:}\quad|u_{\mathrm{app}}^{*}(t)\cdot(Hu_{\mathrm{app}}(t))-u_{0}^{*}\cdot(Hu_{0})|.$
(5.2)

In each case, the time step is adjusted so that the number of FFTs (and their
inverses) are the same for all methods (specifically, 1572864), so that the
computational cost of all schemes is similar.

Figure 3 shows the corresponding results obtained by palindromic schemes with
the coefficients $a_{j}$ real, $\Psi_{P,r}^{[4]}$, and $a_{j}$ complex,
$\Psi_{P,c}^{[4]}$, together with the symmetric-conjugate method
$\Psi_{SC,c}^{[3]}$ with $a_{j}\in\mathbb{C}$. We notice that the qualitative
behavior of all of them is similar: after some point, depending on the
particular step size, the unitarity is lost and the error in energy grows
rapidly.

Figure 3:  Error in norm of the approximate solution (left) and error in
energy (5.2) (right) for the quartic potential (5.1) obtained by the
palindromic schemes $\Psi_{P,r}^{[4]}$ (magenta, dashed line),
$\Psi_{P,c}^{[4]}$ (blue dotted line) and the symmetric-conjugate method
$\Psi_{SC,c}^{[3]}$ (black solid line) along the integration interval. The
step size is chosen so that all methods have the same computational cost.

We notice here the same type of behavior observed in the case of the group
$\mathrm{SU}(2)$: palindromic schemes with both real and complex coefficients
$a_{j}$ are unable to preserve unitarity. On the other hand, symmetric-
conjugate schemes with $a_{j}\in\mathbb{C}$ lead also to unbounded errors,
according with the comments formulated at the beginning of section 4.

We collect in Figure 4 the corresponding results achieved by the palindromic
method $\Xi_{P,r}^{[4]}$ (blue dotted line), and the symmetric-conjugate
schemes $\Psi_{SC,r}^{[3]}$ (black solid line) and $\Xi_{SC,r}^{[4]}$ (magenta
dashed line), all of them with real parameters $a_{j}$. It is worth noticing
that both the norm of the solution and the expected value of the energy are
preserved for very long times by symmetric-conjugate methods with
$a_{j}\in\mathbb{R}$, and this happens even if the method is obtained by
composing a palindromic scheme (with a poor behavior) with its complex
conjugate. By contrast, a symmetric-conjugate method looses its good
preservation properties when composed to form a palindromic scheme, even if
all $a_{j}$ are real.

We have carried out the same experiment, but with the roles of $T$ and $V$
interchanged. In other words, the complex coefficients $b_{j}$ are now
multiplying the discretized Laplacian. In that case, the errors obtained by
all the previous schemes grow unbounded. This indicates that, at least for
this example, one needs both symmetric-conjugate schemes and real coefficients
multiplying the Laplacian to get bounded errors in the preservation of
unitarity and energy.

Figure 4:  Error in norm of the approximate solution (left) and error in
energy (5.2) (right) for the quartic potential (5.1) obtained by the
palindromic scheme $\Xi_{P,r}^{[4]}$ (blue dotted line), and the symmetric-
conjugate schemes $\Psi_{SC,r}^{[3]}$ (black solid line) and
$\Xi_{SC,r}^{[4]}$ (magenta dashed line) along the integration interval. The
step size is chosen so that all methods have the same computational cost.

#### Pöschl–Teller potential.

The next set of simulations is carried out with the well known one-dimensional
Pöschl–Teller potential,

$V(x)=-\frac{\lambda(\lambda+1)}{2}\mbox{sech}^{2}(x),$ (5.3)

with $\lambda(\lambda+1)=10$. It has been used in polyatomic molecular
simulations and admits an analytic treatment [12]. We take again as initial
condition $\psi_{0}(x)=\sigma\,\mathrm{e}^{-x^{2}/2}$, with $\sigma$ a
normalizing constant, then apply Fourier spectral methods on the interval
$x\in[-8,8]$ and integrate until the final time $t_{f}=8000$ with the previous
numerical splitting methods. For this potential we take $N=512$ subdivisions
of the space interval to better visualize the behavior of the different
methods. Figure 5 is the analogous of Fig. 3), and only displays the results
obtained by $\Psi_{SC,c}^{[3]}$ (black solid line) and $\Psi_{P,c}^{[4]}$
(blue dotted line), since the output corresponding to $\Psi_{P,r}^{[4]}$ is
out of the scale (the errors are greater than $10^{87}$). On the other hand,
Figure 6 shows the same pattern as Figure 4: only symmetric-conjugate schemes
with $a_{j}\in\mathbb{R}$ provide bounded errors in the norm and in the energy
of the solution.

Figure 5:  Error in norm of the approximate solution (left) and error in
energy (5.2) (right) for the Pöschl–Teller potential (5.3) obtained by the
palindromic scheme $\Psi_{P,c}^{[4]}$ (blue dotted line) and the symmetric-
conjugate method $\Psi_{SC,c}^{[3]}$ (black solid line) along the integration
interval. The result achieved by $\Psi_{P,r}^{[4]}$ is out of the scale.

Figure 6:  Error in norm of the approximate solution (left) and error in
energy(5.2) (right) for the Pöschl–Teller potential (5.3) obtained by the
palindromic scheme $\Xi_{P,r}^{[4]}$ (blue dotted line), and the symmetric-
conjugate schemes $\Psi_{SC,r}^{[3]}$ (black solid line) and
$\Xi_{SC,r}^{[4]}$ (magenta dashed line) along the integration interval.

Next, we take a shorter final time $t_{f}=100$ and compute the maximum error
in the energy along the time interval for several step sizes $h=\Delta t$ and
integration schemes. The corresponding results are displayed in a log-log
diagrama in Figure 7 (left). The order of each method is clearly visible, as
well as the values of $h$ where instabilities take place. Finally, in Figure 7
(right) we depict the same results but in terms of the computational cost as
measured by the number of FFTs necessary to carry out the calculations. Notice
that, for this range of times, the efficiency of the 4th-order symmetric-
conjugate methods is not far away from the optimized scheme (4.9) that takes
into account the special property (3.7).

Figure 7:  Maximum of error in the expected value of the energy in the
interval $t\in[0,100]$ as a function of the time step (left) and the
computational cost (number of FFTs, right) for several splitting schemes.
Pöschl–Teller potential.

##  6 Concluding remarks

Splitting and composition methods with complex coefficients have shown to be
an appropriate tool in the numerical time integration of differential
equations of parabolic type, when one or more pieces of the equations are only
defined in semi-groups and the aim is to get high accuracy. Since it is
possible to design methods of this class with positive real part, one is thus
able to circumvent the existing order barrier for methods with real
coefficients. In addition, these methods involve smaller truncation errors
than their real counterparts and also exhibit relatively large stability
thresholds. On the other hand, their computational cost notably increases, due
to the use of complex arithmetic.

More recently, it has been shown that the particular class of symmetric-
conjugate methods still exhibits remarkable preservation properties when
applied to differential equations defined by real vector fields and the
solution is projected on the real axis at each integration step. Here we have
extended the analysis to problems evolving in the $\mathrm{SU}(2)$ and more
generally to the numerical integration of the Schrödinger equation, where
preservation of unitarity is a physical requirement. In the former case we
have shown explicitly that symmetric-conjugate splitting methods are indeed
conjugate to unitary methods for sufficiently small time step sizes, and thus
there is not a secular component in the unitarity error propagation.

With respect to the Schrödinger equation, the examples we collect here
indicate that methods of this class (with real coefficients $a_{j}$) could
safely be applied just as other schemes involving only real coefficients for
sufficiently small step sizes, although a general theoretical analysis similar
to the one developed here for problems defined in $\mathrm{SU}(2)$ is lacking
at present. Such analysis is clearly more involved, since one has to take into
account the effect of the space discretization, the possible introduction of
artificial cut-off bounds for unbounded potentials, etc. In this sense, this
paper should be considered as a preliminary step for such analysis. In any
case, we should remark that the use of methods with complex coefficients in
this setting does not imply any extra computational cost, since the problem
has to be treated in the complex domain anyway. Our results show that even
some of the simplest methods within this class provide efficiencies close to
the best standard splitting schemes specifically designed for the integration
of the Schrödinger equation. Although we have limited ourselves here to
methods of order 3 and 4, it is clear that higher order integrators can also
be designed, just by solving the corresponding order conditions [6, 15], and
more efficient schemes can be obtained by taking into account property (3.7)
and the processing technique. It is also worth noticing that, in contrast with
the time integration of parabolic differential equations, here schemes with
real and _negative_ coefficients $a_{j}$ still provide unitary approximations,
and so more efficient schemes with $a_{j}<0$ and $b_{j}\in\mathbb{C}$ might be
possible. All these issues will be treated in a forthcoming paper.

### Acknowledgements

This work has been supported by Ministerio de Ciencia e Innovación (Spain)
through project PID2019-104927GB-C21/AEI/10.13039/501100011033. A.E.-T. has
been additionally funded by the predoctoral contract BES-2017-079697 (Spain).

## References

  * [1] A. Bandrauk, E. Dehghanian, and H. Lu, Complex integration steps in decomposition of quantum exponential evolution operators, Chem. Phys. Lett., 419 (2006), pp. 346–350.
  * [2] S. Blanes and F. Casas, On the necessity of negative coefficients for operator splitting schemes of order higher than two, Appl. Numer. Math., 54 (2005), pp. 23–37.
  * [3] S. Blanes and F. Casas, A Concise Introduction to Geometric Numerical Integration, CRC Press, 2016.
  * [4] S. Blanes, F. Casas, P. Chartier, and A. Escorihuela-Tomàs, On symmetric-conjugate composition methods in the numerical integration of differential equations, Tech. Rep. 2101.04100, arXiv, 2021.
  * [5] S. Blanes, F. Casas, P. Chartier, and A. Murua, Optimized high-order splitting methods for some classes of parabolic equations, Math. Comput., 82 (2013), pp. 1559–1576.
  * [6] S. Blanes, F. Casas, and A. Murua, Splitting and composition methods in the numerical integration of differential equations, Bol. Soc. Esp. Mat. Apl., 45 (2008), pp. 89–145.
  * [7] S. Blanes, F. Casas, and A. Murua, Splitting methods with complex coefficients, Bol. Soc. Esp. Mat. Apl., 50 (2010), pp. 47–61.
  * [8] S. Blanes and P. Moan, Practical symplectic partitioned Runge–Kutta and Runge–Kutta–Nyström methods, J. Comput. Appl. Math., 142 (2002), pp. 313–330.
  * [9] F. Casas, P. Chartier, A. Escorihuela-Tomàs, and Y. Zhang, Compositions of pseudo-symmetric integrators with complex coefficients for the numerical integration of differential equations, J. Comput. Appl. Math., 381 (2021), p. 113006.
  * [10] F. Castella, P. Chartier, S. Descombes, and G. Vilmart, Splitting methods with complex times for parabolic equations, BIT Numer. Math., 49 (2009), pp. 487–508.
  * [11] J. Chambers, Symplectic integrators with complex time steps, Astron. J., 126 (2003), pp. 1119–1126.
  * [12] S. Flügge, Practical Quantum Mechanics, Springer, 1971.
  * [13] A. Galindo and P. Pascual, Quantum Mechanics, Springer, 1990.
  * [14] F. Goth, Higher order auxiliary field quantum Monte Carlo methods, Tech. Rep. 2009.0449, arXiv, 2020.
  * [15] E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations, Springer-Verlag, Second ed., 2006.
  * [16] E. Hansen and A. Ostermann, Exponential splitting for unbounded operators, Math. Comput., 78 (2009), pp. 1485–1496.
  * [17] E. Hansen and A. Ostermann, High order splitting methods for analytic semigroups exist, BIT Numer. Math., 49 (2009), pp. 527–542.
  * [18] C. Lubich, From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis, European Mathematical Society, 2008.
  * [19] R. McLachlan and R. Quispel, Splitting methods, Acta Numerica, 11 (2002), pp. 341–434.
  * [20] T. Prosen and I. Pizorn, High order non-unitary split-step decomposition of unitary operators, J. Phys. A: Math. Gen., 39 (2006), pp. 5957–5964.
  * [21] Q. Sheng, Solving Partial Differential Equations by Exponential Splitting, PhD thesis, Cambridge University, 1989.
  * [22] M. Suzuki, Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations, Phys. Lett. A, 146 (1990), pp. 319–323.
  * [23] M. Suzuki, General theory of fractal path integrals with applications to many-body theories and statistical physics, J. Math. Phys., 32 (1991), pp. 400–407.


T: n→(

T_1: n→(

T_2: n→(

C_4: n→(

M: n→(

§ CLASSICAL REGULAR MAGNETIC ORDERS

Ref. <cit.> studied regular magnetic orders, magnetic orders that respect all the lattice symmetries modulo global $O(3)^T$ spin transformations (rotations and/or
spin flips). In particular, on triangular, kagome, honeycomb and square lattices, all classical regular magnetic orders are classified. These classical orders can all be realized by a product state, where each spin moment on the lattice can be assigned a definite orientation. In this appendix, we explicitly write down the spin configurations of these classical regular magnetic orders, and the lattice symmetry actions on the order parameters.

In terms of the symmetry breaking pattern of the spin $O(3)^T$ symmetry, there are three types of magnetic orders: collinear, coplanar and non-coplanar. The order parameter of a collinear magnetic order is a three-component vector, $\vec n$, which transforms in the spin-1 representation of the $O(3)^T$ spin symmetry. The order parameters of a coplanar magnetic order consists of two orthonormal three-component vectors, $\vec n_{1,2}$, both transforming in the spin-1 representation of the $O(3)^T$ spin symmetry. The order parameters of a non-coplanar magnetic order consists of three orthonormal three-component vectors, $\vec n_{1,2,3}$, all transforming in the spin-1 representation of the $O(3)^T$ spin symmetry.

We start from the triangular lattice. We will denote the position $\vec r$ of a site on a triangular lattice by its coordinates in the basis of translation vectors of $T_{1,2}$ (see Fig. <ref>), such that $\vec r=x\vec T_1+y\vec T_2$, where $\vec T_{1,2}$ is the translation vector of $T_{1,2}$. Under the $p6m$ symmetry,

T_1: (x, y)→(x+1, y)

T_2: (x, y)→(x, y+1)

C_6: (x, y)→(x-y, x)

M: (x, y)→(y, x)

* There is a single collinear classical regular magnetic order, the ferromagnetic order, where $\vec S(x, y)=\vec n$. Under the $p6m$ symmetry, $\vec n$ is invariant.

* There is a single coplanar classical regular magnetic order, the 120$^\circ$ order, where $\vec S(x, y)=(-1)^{x+y}\cos\frac{\pi(x+y)}{3}\vec n_1+(-1)^{x+y}\sin\frac{\pi(x+y)}{3}\vec n_2$. Under the $p6m$ symmetry,

T_1,2: n⃗_1→-1/2n⃗_1-√(3)/2n⃗_2,

C_6: n⃗_1→n⃗_1,

M: n⃗_1,2→n⃗_1,2

* There are two non-coplanar classical regular magnetic order. The first is the tetrahedral order, where $\vec S(x, y)=(-1)^x\vec n_1+(-1)^y\vec n_2+(-1)^{x+y}\vec n_3$. Under the $p6m$ symmetry,

T_1: n⃗_1→-n⃗_1,

T_2: n⃗_1→n⃗_1,

C_6: n⃗_1→n⃗_2,

M: n⃗_1→n⃗_2,

* The second non-coplanar classical regular magnetic order is the F-umbrella order, where $\vec S(x, y)=(-1)^{x+y}\cos\frac{\pi(x+y)}{3}\sin\theta\vec n_1+(-1)^{x+y}\sin\frac{\pi(x+y)}{3}\sin\theta\vec n_2+\cos\theta\vec n_3$, with $\theta$ a free parameter. Under the $p6m$ symmetry,

T_1,2: n⃗_1→-1/2n⃗_1+√(3)/2n⃗_2,

C_6: n⃗_1→n⃗_1,

M: n⃗_1,2,3→n⃗_1,2,3

Next we turn to the kagome lattice. Each unit cell in a kagome lattice includes three sites, so the spin configuration will be written as $\vec S_i(x, y)$, where $(x, y)$ labels the position of the unit cell in the same way as the triangular lattice, and $i=1, 2, 3$ represents the site obtained by applying a half translation $\vec T_1/2$, $\vec T_2/2$ and $(\vec T_1+\vec T_2)/2$ to the $C_6$ center of the unit cell, respectively.

* There is a single collinear classical regular magnetic order, the ferromagnetic order, where $\vec S_i(x, y)=\vec n$ for $i=1,2,3$. Under the $p6m$ symmetry, $\vec n$ is invariant.

* There are two coplanar classical regular magnetic orders. The first is the $\vec q=0$ order, where $\vec S_1(x, y)=\vec n_1$, $\vec S_2(x, y)=-\frac{1}{2}\vec n_1+\frac{\sqrt{3}}{2}\vec n_2$, and $\vec S_3(x. y)=-\frac{1}{2}\vec n_1-\frac{\sqrt{3}}{2}\vec n_2$. Under the $p6m$ symmetry,

T_1,2: n⃗_1, 2→n⃗_1, 2

C_6: n⃗_1→-1/2n⃗_1-√(3)/2

M: n⃗_1→-1/2n⃗_1+√(3)/2n⃗_2,

* The second coplanar classical regular magnetic order is the $\vec q=\sqrt{3}\times\sqrt{3}$ order, where $\vec S_1(x, y)=(-1)^{x+y}\cos\frac{\pi(x+y)}{3}\vec n_1+(-1)^{x+y}\sin\frac{\pi(x+y)}{3}\vec n_2$, $\vec S_2(x, y)=\vec S_1(x, y)$, and $\vec S_3(x, y)=(-1)^{x+y}\cos\frac{\pi(x+y+2)}{3}\vec n_1+(-1)^{x+y}\sin\frac{\pi(x+y+2)}{3}\vec n_2$. Under the $p6m$ symmetry,

T_1, 2: n⃗_1→-1/2n⃗_1-√(3)/2n⃗_2,

C_6: n⃗_1→-1/2n⃗_1+√(3)/2n⃗_2,

M: n⃗_1, 2→n⃗_1, 2

* There are five non-coplanar classical regular magnetic orders. The first is the octahedral order, where $\vec S_1(x, y)=(-1)^y\vec n_1$, $\vec S_2(x, y)=(-1)^x\vec n_2$ and $\vec S_3(x, y)=(-1)^{x+y}\vec n_3$. Under the $p6m$ symmetry,

T_1: n⃗_1→n⃗_1,

T_2: n⃗_1→-n⃗_1,

C_6: n⃗_1→n_3,

M: n⃗_1→n⃗_2,

* The second non-coplanar classical regular magnetic order is the cuboc1 order, where $\vec S_1(x, y)=(-1)^x\vec n_2+(-1)^{x+y}\vec n_3$, $\vec S_2(x, y)=(-1)^y\vec n_1+(-1)^{x+y}\vec n_3$ and $\vec S_3(x, y)=-(-1)^x\vec n_2-(-1)^y\vec n_1$. Under the $p6m$ symmetry,

T_1: n⃗_1→n⃗_1,

T_2: n⃗_1→-n⃗_1,

C_6: n⃗_1→-n⃗_3,

M: n⃗_1→n⃗_2,

* The third non-coplanar classical regular magnetic order is the cuboc2 order, where $\vec S_1(x, y)=(-1)^x\vec n_2-(-1)^{x+y}\vec n_3$, $\vec S_2(x, y)=(-1)^y\vec n_1+(-1)^{x+y}\vec n_3$ and $\vec S_3(x, y)=(-1)^x\vec n_2+(-1)^y\vec n_1$. Under the $p6m$ symmetry,

T_1: n⃗_1→n⃗_1,

T_2: n⃗_1→-n⃗_1,

C_6: n⃗_1→n⃗_3,

M: n⃗_1→n⃗_2,

* The fourth non-coplanar is the $\vec q=0$ umbrella order, where $\vec S_1(x, y)=\sin\theta\vec n_1+\cos\theta\vec n_3$, $\vec S_2(x, y)=-\frac{1}{2}\sin\theta\vec n_1+\frac{\sqrt{3}}{2}\sin\theta\vec n_2+\cos\theta\vec n_3$, and $\vec S_3(x, y)=-\frac{1}{2}\sin\theta\vec n_1-\frac{\sqrt{3}}{2}\sin\theta\vec n_2+\cos\theta\vec n_3$, with $\theta$ a free parameter. Under the $p6m$ symmetry,

T_1, 2: n⃗_1,2,3→n⃗_1,2,3

C_6: n⃗_1→-1/2n⃗_1-√(3)/2n⃗_2,

M: n⃗_1→-1/2n⃗_1+√(3)/2n⃗_2,

* The last non-coplanar classical regular magnetic order is the $\vec q=\sqrt{3}\times\sqrt{3}$ umbrella order, where $\vec S_1(x, y)=(-1)^{x+y}\cos\frac{\pi(x+y)}{3}\sin\theta\vec n_1+(-1)^{x+y}\sin\frac{\pi(x+y)}{3}\sin\theta\vec n_2+\cos\theta\vec n_3$, $\vec S_2(x, y)=\vec S_1(x, y)$, and $\vec S_3(x, y)=-(-1)^{x+y}\cos\frac{\pi(x+y-1)}{3}\sin\theta\vec n_1-(-1)^{x+y}\sin\frac{\pi(x+y-1)}{3}\sin\theta\vec n_2+\cos\theta\vec n_3$. Under the $p6m$ symmetry,

T_1, 2: n⃗_1→-1/2n⃗_1-√(3)/2n⃗_2,

C_6: n⃗_1→-1/2n⃗_1+√(3)/2n⃗_2,

M: n⃗_1,2,3→n⃗_1,2,3

Now we turn to the honeycomb lattice. Each unit cell of a honeycomb lattice includes two sites, so the spin configuration will be written in terms of $\vec S_A(x, y)$ and $\vec S_B(x, y)$, where the $A$ and $B$ sublattice can be obtained by translating by $\frac{2\vec T_1+\vec T_2}{3}$ and $\frac{\vec T_1-\vec T_2}{3}$ from the $C_6$ center, respectively.

* There are two collinear classical regular magnetic orders. The first is the ferromagnetic order, where $\vec S_A(x, y)=\vec S_B(x, y)=\vec n$. Under the $p6m$ symmetry, $\vec n$ is invariant.

* The second collinear classical regular magnetic order is the anti-ferromagnetic order, where $\vec S_A(x, y)=-\vec S_B(x, y)=\vec n$. Under the $p6m$ symmetry,

T_1, 2: n⃗→n⃗

C_6: n⃗→-n⃗

M: n⃗→-n⃗

* There is a single coplanar classical regular magnetic order, the $V$ order, where $\vec S_A(x, y)=\cos\theta\vec n_1-\sin\theta\vec n_2$ and $\vec S_B(x, y)=\cos\theta\vec n_1+\sin\theta\vec n_2$, with $\theta$ a free parameter. Under the $p6m$ symmetry,

T_1,2: n⃗_1, 2→n⃗_1, 2

C_6: n⃗_1→n⃗_1,

M: n⃗_1→n⃗_1,

* There are two non-coplanar classical regular magnetic orders. The first is the cubic order, where $\vec S_A(x, y)=(-1)^x\vec n_1+(-1)^y\vec n_2+(-1)^{x+y}\vec n_3$ and $\vec S_B(x, y)=(-1)^x\vec n_1-(-1)^y\vec n_2+(-1)^{x+y}\vec n_3$. Under the $p6m$ symmetry,

T_1: n⃗_1→-n⃗_1,

T_2: n⃗_1→n⃗_1,

C_6: n⃗_1→n⃗_2,

M: n⃗_1→n⃗_2,

* The second non-coplanar classical regular magnetic order is the tetrahedral order, where $\vec S_A(x, y)=(-1)^x\vec n_1+(-1)^y\vec n_2+(-1)^{x+y}\vec n_3$ and $\vec S_B(x, y)=-(-1)^x\vec n_1+(-1)^y\vec n_2-(-1)^{x+y}\vec n_3$. Under the $p6m$ symmetry,

T_1: n⃗_1→-n⃗_1,

T_2: n⃗_1→n⃗_1,

C_6: n⃗_1→-n⃗_2,

M: n⃗_1→-n⃗_2,

Finally, we discuss the square lattice. We will denote the position of a site by its coordinates in the basis of translation vectors of $T_{1, 2}$ (see Fig. <ref>). such that $\vec r=x\vec T_1+y\vec T_2$, where $\vec T_{1,2}$ is the translation vector of $T_{1,2}$. Under the $p4m$ symmetry,

T_1: (x, y)→(x+1, y)

T_2: (x, y)→(x, y+1)

C_4: (x, y)→(-y, x)

M: (x, y)→(-x, y)

* There are two collinear classical regular magnetic orders. The first is the ferromagnetic order, where $\vec S(x, y)=\vec n$. Under the $p4m$ symmetry, $\vec n$ is invariant.

* The second collinear classical regular magnetic order is the anti-ferromagnetic order, where $\vec S(x, y)=(-1)^{x+y}\vec n$. Under the $p4m$ symmetry,

T_1, 2: n⃗→-n⃗

C_4: n⃗→-n⃗

M: n⃗→n⃗

* There are two coplanar classical regular magnetic orders. The first is the orthogonal order, where $\vec S(x, y)=\frac{(-1)^x+(-1)^y}{2}\vec n_1+\frac{-(-1)^x+(-1)^y}{2}\vec n_2$. Under the $p4m$ symmetry,

T_1: n⃗_1→n⃗_2,

T_2: n⃗_1→-n⃗_2,

C_4: n⃗_1→n⃗_1,

M: n⃗_1,2→n⃗_1,2

* The second coplanar classical regular magnetic order is the $V$ order, where $\vec S(x, y)=\cos\theta\vec n_1-(-1)^{x+y}\sin\theta\vec n_2$, where $\theta$ a free parameter. Under the $p4m$ symmetry,

T_1, 2: n⃗_1→n⃗_1,

C_4: n⃗_1, 2→n⃗_1, 2

M: n⃗_1, 2→n⃗_1, 2

* There are two non-coplanar classical regular magnetic orders. The first is the tetrahedral umbrealla order (also known as the AF umbrella order), where $\vec S(x, y)=\frac{(-1)^x\sin\theta}{\sqrt{2}}\vec n_1-\frac{(-1)^y\sin\theta}{\sqrt{2}}\vec n_2-(-1)^{x+y}\cos\theta\vec n_3$, with $\theta$ a free parameter. Under the $p4m$ symmetry,

T_1: n⃗_1→-n⃗_1,

T_2: n⃗_1→n⃗_1,

C_4: n⃗_1→-n⃗_2,

M: n⃗_1, 2 ,3→n⃗_1, 2, 3

* The second non-coplanar classical regular magnetic order is the umbrella order (also known as the F umbrella order), where $\vec S(x, y)=\cos\theta\vec n_1+\frac{(-1)^x\sin\theta}{\sqrt{2}}\vec n_2+\frac{(-1)^y\sin\theta}{\sqrt{2}}\vec n_3$. Under the $p4m$ symmetry,

T_1: n⃗_1→n⃗_1,

T_2: n⃗_1→n⃗_1,

C_4: n⃗_1→n_1,

M: n⃗_1, 2, 3→n⃗_1, 2, 3

§ STABILITY OF DSL REALIZATIONS ON NAYBO$_2$ AND TWISTED BILAYER WSE$_2$

In this appendix, we discuss the stability of a few more examples of DSL realizations on systems with spin-orbit coupling (SOC). The specific systems we have in mind are NaYbO$_2$ and twisted bilayer WSe$_2$ (tWSe$_2$). Recently, it was pointed out that tWSe$_2$ is a good quantum simulator of triangular lattice Hubbard model, which can be effectively described by a triangular lattice spin-1/2 system in the strong coupling regime <cit.>.

The symmetries of NaYbO$_2$ are given in Eq. (<ref>). The symmetries of tWSe$_2$ are
T_1,2, C_3≡C_6^2, SO(2), T
where $SO(2)$ is a reduced spin rotational symmetry [These are the symmetries of tWSe$_2$ in the presence of a displacement field, which is satisfied in the generic experimental setting. If there is no displacement field, there is an extra mirror symmetry.].

On triangular lattice spin-1/2 systems with the full $p6m\times O(3)^T$ symmetry, our exhaustive search finds 3 realizations of DSL, given by Eqs. (<ref>), (<ref>) and (<ref>). Using these symmetry actions, it is straightforward to see that for all three realizations, the remaining symmetries of NaYbO$_2$ are sufficient to forbid all relevant operators of DSL listed in Sec. <ref>. However, for the symmetry setting of tWSe$_2$ and for all three realizations, the $(A_L, A_R)$ operator (the fermion mass that transforms in the adjoint representation of the flavor symmetry) is always symmetry-allowed and will destablize the DSL. This means if a DSL is stably realized in tWSe$_2$, that realization cannot be compatible with a full $p6m\times O(3)^T$ symmetry.