# Symmetric-conjugate splitting methods for evolution equations of parabolic
type

###### Abstract.

The present work provides a comprehensive study of symmetric-conjugate
operator splitting methods in the context of linear parabolic problems and
demonstrates their additional benefits compared to symmetric splitting
methods. Relevant applications include nonreversible systems and ground state
computations for linear Schrödinger equations based on the imaginary time
propagation. Numerical examples confirm the favourable error behaviour of
higher-order symmetric-conjugate splitting methods and illustrate the
usefulness of a time stepsize control, where the local error estimation relies
on the computation of the imaginary parts and thus requires negligible costs.

###### Key words and phrases:

Linear evolution equations, Parabolic problems, Schrödinger equations,
Operator splitting methods, Fourier spectral methods, Stability, Convergence,
Efficiency.

###### 1991 Mathematics Subject Classification:

65J10, 65L04, 65M12.

∗Corresponding author: Mechthild Thalhammer

S. Blanes${}^{{\href mailto:serblaza@imm.upv.es}\,1}$, F. Casas${}^{{\href
mailto:fernando.casas@uji.es}\,2}$, C. González${}^{{\href
mailto:cesareo.gonzalez@uva.es}\,3}$, M. Thalhammer${}^{{\href
mailto:mechthild.thalhammer@uibk.ac.at}\,*4}$

1Universitat Politècnica de València, Instituto de Matemática
Multidisciplinar, 46022 Valencia, Spain.

2Universitat Jaume I, IMAC and Departament de Matemàtiques, 12071 Castellón,
Spain.

3Universidad de Valladolid, Departamento de Matemática Aplicada, 47011
Valladolid, Spain.

4Universität Innsbruck, Institut für Mathematik, 6020 Innsbruck, Austria.

##  1\. Introduction

A wide range of mathematical models for dynamical processes involves initial
value problems for ordinary and partial differential equations. Specifically
designed space and time discretisation methods are of major importance in view
of their effective simulation. Over the last decades, a variety of
contributions has established theoretical and numerical evidence that the
class of operator splitting methods leads to favourable time integration
methods and additionally preserves structural properties of linear and
nonlinear evolution equations.

_General references._ The monographs [35, 41] give comprehensive overviews of
applications in quantum physics. Expositions of approved functional analytical
frameworks within the scope of parabolic and Schrödinger equations are found
in [24, 33, 39, 42]. For detailed information on splitting and composition
methods, we refer to [29, 40, 46].

_Objectives._ In the present work, we are concerned with the study of higher-
order symmetric-conjugate splitting methods in comparison with symmetric
splitting methods. For both classes of time integration methods, certain
positivity conditions on the real parts of the complex coefficients are
crucial to ensure stability for nonreversible systems. Moreover, provided that
the problem data satisfy suitable regularity and consistency requirements, it
is ensured that the nonstiff orders of convergence are retained. We provide
theoretical results and numerical illustrations, which confirm the reliability
and efficiency of symmetric-conjugate splitting methods for relevant model
problems including parabolic counterparts of linear Schrödinger equations. We
exemplify the general observation that the numerical evolution operator
associated with a symmetric-conjugate splitting method inherits the
fundamental property of self-adjointness, in contrast to complex splitting
methods with a strict symmetry in the configuration of their coefficients. We
conclude that a particular benefit of symmetric-conjugate splitting methods is
the possibility to base local error estimations on the computation of
imaginary parts.

_Related works._ Our investigations are inspired by a series of works on real
and complex splitting methods as well as modified splitting methods. As an
excerpt, we mention [1, 3, 4, 5, 6, 8, 9, 18, 19, 20, 21, 22, 23, 25, 27, 31,
32, 34, 36, 37, 38, 43, 44, 47, 48, 51] and our former contributions [2, 7,
10, 11, 12, 14, 15, 16, 17, 49, 50], where further references are given.

_Outline._ The present work has the following structure. In Section 2, we
introduce operator splitting methods and give a brief summary of fundamental
concepts for their convergence analysis and practical implementation.
Moreover, we introduce a unifying formulation for relevant model problems such
as parabolic analogues of Schrödinger equations. In Section 3, we demonstrate
the benefits of symmetric-conjugate splitting methods over symmetric splitting
methods. For this purpose, in order to reduce the amount of technicalities, it
is useful to employ the generally understandable setting of real symmetric
matrices. Further explanations on the construction of higher-order splitting
methods by composition and numerical comparisons of symmetric-conjugate
splitting methods with standard and modified splitting methods are finally
given in Section 4.

##  2\. Fundamental concepts and model problems

The starting point of our considerations is the linear evolution equation

$\begin{cases}u^{\prime}(t)=F\big{(}u(t)\big{)}=A\,u(t)+B\,u(t)\,,\quad
t\in[t_{0},T]\,,\\\ u(t_{0})\text{ given}\,.\end{cases}$ (1)

Throughout, we denote by $(X,\|\cdot\|_{X})$ the underlying Banach space and
assume that the domains of the operators $A:D(A)\subset X\to X$ and
$B:D(B)\subset X\to X$ have a non-empty intersection. Under the presumption of
suitably restricted domains, the commutator of linear operators is defined by

$\big{[}A,B\big{]}=A\,B-BA\,.$

More generally, iterated commutators are determined recursively

$\begin{gathered}\text{ad}_{A}^{\,\ell}(B)=\begin{cases}B\,,&\ell=0\,,\\\
\big{[}A,\text{ad}_{A}^{\,\ell-1}(B)\big{]}\,,&\ell\in\mathbb{N}_{\geq
1}\,.\end{cases}\end{gathered}$

_Splitting approach._ For simplicity, we restrict ourselves to uniform time
grids

$t_{n}=t_{0}+nh\,,\quad n\in\\{0,1,\dots,N\\}\,,\quad
h=\tfrac{T-t_{0}}{N}>0\,,$ (2a) defined by positive integer numbers
$N\in\mathbb{N}_{\geq 1}$. As usual in a time-stepping approach, our aim is to
compute approximations to the exact solution values through a recurrence
relation of the form $u_{n}=\mathcal{S}^{(F)}_{h}(u_{n-1})\approx
u(t_{n})=\mathcal{E}^{(F)}_{h}\big{(}u(t_{n-1})\big{)}\,,\quad
n\in\\{1,\dots,N\\}\,,$ (2b) where $\mathcal{E}^{(F)}$ and $\mathcal{S}^{(F)}$
represent the exact and numerical evolution operators, respectively.

Operator splitting methods rely on the idea to treat the subproblems that
arise from the natural decomposition in (1) separately and to compose their
solutions in a favourable manner. With regard to (2b), we consider a single
subinterval and denote the evolution operators associated with

$\begin{split}&\begin{cases}v^{\prime}(t)=A\,u(t)\,,\quad
t\in[t_{n-1},t_{n}]\,,\\\ v(t_{n-1})\text{ given}\,,\quad
v(t_{n})=\mathcal{E}^{(A)}_{h}\big{(}v(t_{n-1})\big{)}\,,\end{cases}\\\
&\begin{cases}w^{\prime}(t)=B\,w(t)\,,\quad t\in[t_{n-1},t_{n}]\,,\\\
w(t_{n-1})\text{ given}\,,\quad
w(t_{n})=\mathcal{E}^{(B)}_{h}\big{(}w(t_{n-1})\big{)}\,,\end{cases}\end{split}$

by $\mathcal{E}^{(A)}$ and $\mathcal{E}^{(B)}$, respectively. Incorporating
suitably chosen real or more generally complex coefficients

$a_{j},b_{j}\in\mathbb{C}\,,\quad j\in\\{1,\dots,s\\}\,,$ (2c)

splitting methods for (1) can be cast into the format (2b) with

$\mathcal{S}^{(F)}_{h}=\mathcal{E}^{(b_{s}B)}_{h}\circ\mathcal{E}^{(a_{s}A)}_{h}\circ\dots\circ\mathcal{E}^{(b_{1}B)}_{h}\circ\mathcal{E}^{(a_{1}A)}_{h}\,.$
(2d)

Throughout, we denote by $p\in\mathbb{N}$ the classical order of a splitting
method and tacitly assume that the coefficients of the considered schemes
fulfill the elementary consistency condition

$\sum_{j=1}^{s}a_{j}=1\,,\quad\sum_{j=1}^{s}b_{j}=1\,.$ (2e)

On occasion, when the structural characteristics are essential, we use the
compact symbolic notation

$\mathcal{S}^{(F)}_{h}=h\,\big{(}b_{s},a_{s},\dots,b_{1},a_{1}\big{)}\,.$ (2f)

_Symmetric-conjugate methods._ We primarily focus on symmetric-conjugate
operator splitting methods of the form

$\begin{gathered}s=2\,r\,,\quad a_{1}=0\,,\\\
a_{s+2-j}=\overline{a_{j}}\,,\quad j\in\\{2,\dots,r\\}\,,\quad
b_{s+1-j}=\overline{b_{j}}\,,\quad j\in\\{1,\dots,r\\}\,,\\\
\mathcal{S}^{(F)}_{h}=h\,\big{(}\overline{b_{1}},\overline{a_{2}},\overline{b_{2}},\dots,\overline{a_{r}},\overline{b_{r}},a_{r+1},b_{r},a_{r},\dots,b_{2},a_{2},b_{1},0\big{)}\,,\end{gathered}$
(3a) and impose the positivity conditions
$\begin{gathered}a_{j}\in\mathbb{C}\,,\quad\Re(a_{j})>0\,,\quad
j\in\\{2,\dots,r+1\\}\,,\\\ b_{j}\in\mathbb{C}\,,\quad\Re(b_{j})>0\,,\quad
j\in\\{1,\dots,r\\}\,,\end{gathered}$ (3b)

to ensure well-definedness and thus stability for evolution equations of
parabolic type. It is natural to contrast symmetric-conjugate with symmetric
splitting methods

$\begin{gathered}s=2\,r\,,\quad a_{1}=0\,,\\\ a_{s+2-j}=a_{j}\,,\quad
j\in\\{2,\dots,r\\}\,,\quad b_{s+1-j}=b_{j}\,,\quad j\in\\{1,\dots,r\\}\,,\\\
\mathcal{S}^{(F)}_{h}=h\,\big{(}b_{1},a_{2},b_{2},\dots,a_{r},b_{r},a_{r+1},b_{r},a_{r},\dots,b_{2},a_{2},b_{1},0\big{)}\,.\end{gathered}$
(4)

For notational simplicity, we henceforth omit coefficients that are equal to
zero and write for instance

$\begin{gathered}\mathcal{S}^{(F)}_{h}=h\,\big{(}\overline{b_{1}},\overline{a_{2}},\overline{b_{2}},\dots,\overline{a_{r}},\overline{b_{r}},a_{r+1},b_{r},a_{r},\dots,b_{2},a_{2},b_{1}\big{)}\,,\\\
\mathcal{S}^{(F)}_{h}=h\,\big{(}b_{1},a_{2},b_{2},\dots,a_{r},b_{r},a_{r+1},b_{r},a_{r},\dots,b_{2},a_{2},b_{1}\big{)}\,,\end{gathered}$

for short, see (3) and (4). The different behaviour exhibited by both classes
of splitting methods and the particularly favourable performance of symmetric-
conjugate schemes in the time integration of linear ordinary differential
equations that are defined by real symmetric matrices deserves a detailed
analysis, which is carried out in this work.

_Elementary splitting methods._ As elementary instances, we introduce the
famous Lie–Trotter and Strang splitting methods and a third-order symmetric-
conjugate splitting method. We recall that the positive integers
$s\in\mathbb{N}_{\geq 1}$ and $p\in\mathbb{N}_{\geq 1}$ denote the number of
stages and the classical order of a splitting method (2).

  1. (i)

The simplest first-order scheme

$\begin{gathered}p=1\,,\quad s=1\,,\quad a_{1}=1\,,\quad b_{1}=1\,,\\\
\mathcal{S}^{(F)}_{h}=\mathcal{E}^{(B)}_{h}\circ\mathcal{E}^{(A)}_{h}\approx\mathcal{E}^{(F)}_{h}\,,\end{gathered}$

is known as Lie–Trotter splitting method. Evidently, it fulfils the positivity
condition $a_{1},b_{1}>0$, but it does not fit into the classes (3) or (4),
respectively.

  2. (ii)

The second-order Strang splitting method, which comprises two stages and the
symmetric composition

$\begin{gathered}p=2\,,\quad s=2\,,\quad a_{1}=0\,,\quad a_{2}=1\,,\quad
b_{1}=\tfrac{1}{2}=b_{2}\,,\\\
\mathcal{S}^{(F)}_{h}=\mathcal{E}^{(\frac{1}{2}B)}_{h}\circ\mathcal{E}^{(A)}_{h}\circ\mathcal{E}^{(\frac{1}{2}B)}_{h}\approx\mathcal{E}^{(F)}_{h}\,,\end{gathered}$
(5)

is contained in both classes (3) and (4) with $r=1$.

  3. (iii)

The simplest symmetric-conjugate splitting method of order three

$\begin{gathered}p=3\,,\quad s=3\,,\quad a_{1}=0\,,\quad
a_{2}=\tfrac{1}{2}\,\big{(}1+\mathrm{i}\,\tfrac{1}{\sqrt{3}}\big{)}\,,\quad
b_{1}=\tfrac{1}{2}\,a_{2}\,,\quad b_{2}=\tfrac{1}{2}\,,\\\
\mathcal{S}^{(F)}_{h}=\mathcal{E}^{(\overline{b_{1}}B)}_{h}\circ\mathcal{E}^{(\overline{a_{2}}A)}_{h}\circ\mathcal{E}^{(b_{2}B)}_{h}\circ\mathcal{E}^{(a_{2}A)}_{h}\circ\mathcal{E}^{(b_{1}B)}_{h}\approx\mathcal{E}^{(F)}_{h}\,,\end{gathered}$

was proposed in [3]. Alternatively, it is retained as a special double jump
composition of the Strang splitting method

$\mathcal{S}^{(F)}_{h}=\mathcal{S}^{(F,\,\text{Strang})}_{\overline{a_{2}}\,h}\circ\mathcal{S}^{(F,\,\text{Strang})}_{a_{2}\,h}=h\,\big{(}\tfrac{1}{2}\,\overline{a_{2}},\overline{a_{2}},\tfrac{1}{2}\,(a_{2}+\overline{a_{2}}),a_{2},\tfrac{1}{2}\,a_{2}\big{)}\,.$

_Higher-order splitting methods._ In our numerical tests, detailed in Section
4, we compare higher-order standard and modified operator splitting methods
involving real coefficients with complex symmetric and symmetric-conjugate
splitting methods, see Figure 1. For the convenience of the readers, we
display the method coefficients of the symmetric-conjugate schemes with
corresponding denominations in Figures 9 and 10. Besides, a link to a Matlab
code is provided in Section 4.

A fourth-order symmetric scheme by Yoshida [51] comprises negative
coefficients, wherefore instabilities arise for evolution equations of
parabolic type, see (22a). Its complex analogue given in (22b) leads to a
stable alternative. A non-standard scheme is Chin’s fourth-order modified
potential operator splitting method [20] involving positive coefficients and
double commutators, see also [16, 17].

Moreover, we apply different complex symmetric and symmetric-conjugate
splitting methods, which satisfy the stability condition that all coefficients
have non-negative real parts. Optimised symmetric schemes of order four and a
symmetric scheme of order six are found in [12].111See also
www.gicas.uji.es/Research/splitting-complex.html. Amongst the symmetric-
conjugate splitting methods, we highlight sixth-order schemes, recently
proposed in [7] for linear unitary problems. They set aside the second- and
fourth-order barriers for standard and modified splitting methods, see [1, 10,
47, 48, 25] and references given therein. We point out that these schemes are
characterised by positive coefficients $a_{j}>0$ for $j\in\\{2,\dots,s\\}$.
Consequently, they are suitable for the time integration of different classes
of evolution equations including parabolic and Schrödinger equations, see (6)
as well as (7) below and Table 1.

_Model problems._ As prototype models for linear evolution equations of
parabolic type, we study partial differential equations that involve the
Laplacian and a potential. This kind of nonreversible systems in particular
arises in ground and excited state computations based on the imaginary time
propagation, see [35, 41] for detailed information on the theoretical
foundations.

In the context of these relevant applications, it is often reasonable to
assume that solutions are localised and highly regular. Thus, it is justified
to restrict the originally unbounded space domain to a Cartesian product of
sufficiently large intervals and to impose periodic boundary conditions.
Conveniences of such settings are that fast Fourier techniques can be applied
and that higher-order splitting methods do not suffer from severe order
reductions. Throughout, we base our considerations on these assumptions.

In the following, we denote by $\Omega\subseteq\mathbb{R}^{d}$ the underlying
space domain, by $\Delta=\partial_{x_{1}}^{2}+\dots+\partial_{x_{d}}^{2}$ the
Laplacian with respect to the spatial variables
$x=(x_{1},\dots,x_{d})\in\Omega$, and by $V:\Omega\to\mathbb{R}$ a space-
dependent real-valued potential acting as multiplication operator. For
notational simplicity, we omit scaling constants and signs, unless they are
significant with regard to classifications as parabolic or Schrödinger
equations.

As a test problem, we consider the linear parabolic problem

$\begin{cases}\partial_{t}U(x,t)=\tfrac{1}{2}\,\Delta\,U(x,t)-V(x)\,U(x,t)\,,\\\
U(x,t_{0})\text{ given}\,,\quad(x,t)\in\Omega\times[t_{0},T]\,,\end{cases}$
(6)

for a real-valued solution $U:\Omega\times[t_{0},T]\to\mathbb{R}$. This
special choice is justified by the imaginary time propagation of the linear
Schrödinger equation

$\begin{cases}\mathrm{i}\,\partial_{t}\Psi(x,t)=-\,\tfrac{1}{2}\,\Delta\Psi(x,t)+V(x)\,\Psi(x,t)\,,\\\
\Psi(x,t_{0})\text{ given}\,,\quad(x,t)\in\Omega\times[t_{0},T]\,,\end{cases}$
(7)

with complex-valued wave function $\Psi:\Omega\times[t_{0},T]\to\mathbb{C}$,
see Section 4 for further explanations.

Evidently, the above model problems (6) and (7) can be cast into the unifying
formulation

$\begin{cases}&\partial_{t}U(x,t)=\alpha\,\Delta\,U(x,t)+\beta\,V(x)\,U(x,t)\,,\\\
&U(x,t_{0})\text{ given}\,,\quad(x,t)\in\Omega\times[t_{0},T]\,,\end{cases}$
(8a) with $\alpha\,,\beta\in\mathbb{C}$ denoting certain constants and
$U:\Omega\times[t_{0},T]\to\mathbb{C}$ the solution. Specifically, we chose
the arising quantities as
$\displaystyle\alpha=\tfrac{1}{2}\,,\quad\beta=-\,1\,,$ (8b)
$\displaystyle\alpha=\tfrac{1}{2}\,\mathrm{i}\,,\quad\beta=-\,\mathrm{i}\,.$
(8c)

Additional scaling constants and signs that are insignificant with regard to
classifications as parabolic or Schrödinger equations are again neglected.

_Practicable presumptions._ In this work, we are primarily interested in the
time integration of parabolic initial-boundary value problems by high-order
operator splitting methods. Hence, to ensure that the nonstiff orders of
convergence are retained in our numerical comparisons, we restrict ourselves
to situations, where the problem data satisfy suitable regularity and
consistency requirements. Otherwise, even though well-tailored higher-order
schemes are generally more favourable than lower-order schemes, we have to
expect substantial order reductions.

Furthermore, to implement the physically relevant and numerically challenging
case of three space dimensions by means of Fourier spectral space
discretisations, we simplify the general setting in this respect.
Specifically, we presume that it is appropriate to replace the underlying
space domain with a Cartesian product of sufficiently large intervals

$\Omega=[-\,a,a]^{d}\,,\quad a>0\,,\quad d\in\\{1,2,3\\}\,,$ (9a) and
prescribe Gaussian-like initial states
$U(x,t_{0})=c_{1}\,\mathrm{e}^{\,-\,c_{2}\,|x-c_{3}|^{2}}\,,\quad
x=(x_{1},\dots,x_{d})\in\Omega\,,$ (9b) which fulfil intrinsic periodicity
conditions on $\Omega$ with high accuracy.

In the following, we sketch the employed means that permit efficient
implementations of operator splitting methods for the model problem (8a) by
fast Fourier techniques. For detailed descriptions, we refer the readers to
our former works [16, 17].

_Fourier series representations._ We denote by
$(\mathcal{F}_{m})_{m\in\mathbb{Z}^{d}}$ the Fourier functions with
periodicity domain (9a) and by $(\lambda_{m})_{m\in\mathbb{Z}^{d}}$ the
corresponding real eigenvalues of the Laplace operator

$\begin{gathered}\mathcal{F}_{m}(x)=(2\,a)^{-\,\frac{d}{2}}\;\mathrm{e}^{\,\mathrm{i}\,\pi\,m_{1}\,(\frac{x_{1}}{a}+1)}\cdots\,\mathrm{e}^{\,\mathrm{i}\,\pi\,m_{d}\,(\frac{x_{d}}{a}+1)}\,,\\\
\Delta\,\mathcal{F}_{m}=\lambda_{m}\,\mathcal{F}_{m}\,,\quad\lambda_{m}=-\,\tfrac{\pi^{2}|m|^{2}}{a^{2}}\in\mathbb{R}\,,\\\
x=(x_{1},\dots,x_{d})\in\Omega\,,\quad
m=(m_{1},\dots,m_{d})\in\mathbb{Z}^{d}\,.\end{gathered}$ (9c)

Realisations of representations by Fourier series

$\begin{gathered}v=\sum_{m\in\mathbb{Z}^{d}}v_{m}\,\mathcal{F}_{m}\,,\quad\Delta
v=\sum_{m\in\mathbb{Z}^{d}}\lambda_{m}\,v_{m}\,\mathcal{F}_{m}\,,\\\
v_{m}=\int_{[-\,a,a]^{d}}v(x)\,\mathcal{F}_{-m}(x)\;\mathrm{d}x\,,\quad
m\in\mathbb{Z}^{d}\,,\end{gathered}$ (9d)

are based on suitable truncations of the infinite index sets
$\mathcal{M}\subset\mathbb{Z}^{d}$ such that $|\mathcal{M}|=M\in\mathbb{N}$ as
well as quadrature approximations by the trapezoidal rule.

_Stiff and nonstiff subproblems._ We rewrite the partial differential equation
in (8a) as abstract evolution equation of the form (1) with $u(t)=U(\cdot,t)$
for $t\in[t_{0},T]$. We employ the natural decomposition into a stiff and a
nonstiff partial differential equation

$\begin{cases}&\partial_{t}U(x,t)=\alpha\,\Delta\,U(x,t)\,,\\\
&\partial_{t}U(x,t)=\beta\,V(x)\,U(x,t)\,.\end{cases}$

Here, it should be noted that the considered Fourier spectral space
discretisation affects the definition of the operators.

Accordingly, we assign the first unbounded linear operator with the Laplacian

$A=\alpha\,\Delta\,,\quad\alpha\in\mathbb{R}\,,\quad\alpha>0\,.$ (10a) For any
complex coefficient $a\in\mathbb{C}$ with non-negative real part $\Re(a)\geq
0$, it is ensured that the corresponding subproblem
$\begin{cases}v^{\prime}(t)=a\,A\,v(t)\,,\quad t\in[t_{n-1},t_{n}]\,,\\\
\displaystyle
v(t_{n-1})=\sum_{m\in\mathbb{Z}^{d}}v_{m}(t_{n-1})\,\mathcal{F}_{m}\text{
given}\,,\end{cases}$ (10b) is well-posed and its solution formally given by
the Fourier series representation
$v(t_{n})=\mathcal{E}^{(aA)}_{h}\,v(t_{n-1})=\sum_{m\in\mathbb{Z}^{d}}\mathrm{e}^{\,a\,h\,\alpha\,\lambda_{m}}\,v_{m}(t_{n-1})\,\mathcal{F}_{m}\,,$
(10c)

see also (9).

The second nonstiff subproblem is defined by the potential, which acts as a
multiplication operator, and resolved by pointwise products.

_Evolution equations of Schrödinger type._ It is instructive to observe that
linear Schrödinger equations (7) are included in (8a), see also (8c). For this
reason, we may expect that the convergence analysis provided in [49, 50] for
Fourier spectral space discretisations combined with high-order operator
splitting methods based on real coefficients can be transferred to evolution
equations of parabolic type and complex symmetric-conjugate splitting methods
satisfying a positivity condition. For the convenience of the readers, we next
recall fundamental notions and concepts.

_Local and global errors._ We employ a standard argument based on the
telescopic identity to conclude that the validity of stability bounds combined
with local error expansions implies global error estimates of the form

$\big{\|}u_{n}-u(t_{n})\big{\|}_{X}\leq
C\,\big{(}\big{\|}u_{0}-u(t_{0})\big{\|}_{X}+h^{p}\big{)}\,,\quad
n\in\\{1,\dots,N\\}\,.$

This fundamental principle serves as guideline for our convergence analysis of
operator splitting methods applied to evolution equations of parabolic type.

The general approach is most comprehensible within the context of evolution
equations that are defined by bounded linear operators

$u^{\prime}(t)=F\big{(}u(t)\big{)}=(A+B)\,u(t)\,,\quad t\in[t_{0},T]\,,\quad
A,B:X\longrightarrow X\,,$

since then the exact and splitting solutions can be represented by products of
matrix exponentials

$\mathcal{E}^{(F)}_{h}=\mathrm{e}^{\,h\,(A+B)}\,,\quad\mathcal{S}^{(F)}_{h}=\mathrm{e}^{\,b_{s}hB}\,\mathrm{e}^{\,a_{s}hA}\cdots\mathrm{e}^{\,b_{1}hB}\,\mathrm{e}^{\,a_{1}hA}\,,$

and thus stability bounds follow at once from

$\begin{gathered}\big{\|}\mathcal{E}^{(F)}_{h}\big{\|}_{X\leftarrow
X}\leq\mathrm{e}^{\,Ch}\,,\quad\big{\|}\mathcal{S}^{(F)}_{h}\big{\|}_{X\leftarrow
X}\leq\mathrm{e}^{\,Ch}\,,\\\ C=\max\Bigg{\\{}\|A\|_{X\leftarrow
X}+\|B\|_{X\leftarrow X},\sum_{j=1}^{s}\Big{(}|a_{j}|\|A\|_{X\leftarrow
X}+|b_{j}|\|B\|_{X\leftarrow X}\Big{)}\bigg{\\}}\,.\end{gathered}$

Provided that the coefficients satisfy certain order conditions such that the
local error expansion

$\mathcal{L}^{(F)}_{h}=\mathcal{S}^{(F)}_{h}-\mathcal{E}^{(F)}_{h}=\mathcal{O}\big{(}h^{p+1}\big{)}$

holds, the desired global error estimate follows by means of the relations

$\begin{gathered}u_{n}-u(t_{n})=\big{(}\mathcal{S}^{(F)}_{h}\big{)}^{n}\,\big{(}u_{0}-u(t_{0})\big{)}+\sum_{k=0}^{n-1}\big{(}\mathcal{S}^{(F)}_{h}\big{)}^{n-1-k}\mathcal{L}^{(F)}_{h}\big{(}\mathcal{E}^{(F)}_{h}\big{)}^{k}\,u(t_{0})\,,\\\
\big{\|}u_{n}-u(t_{n})\big{\|}_{X}\leq\mathrm{e}^{\,Ct_{n}}\,\Big{(}\big{\|}u_{0}-u(t_{0})\big{\|}_{X}+n\,\big{\|}\mathcal{L}^{(F)}_{h}\big{\|}_{X\leftarrow
X}\big{\|}u(t_{0})\big{\|}_{X}\Big{)}\,,\\\
n\in\\{1,\dots,N\\}\,.\end{gathered}$

##  3\. Symmetric-conjugate versus symmetric methods

In this section, we contrast the favourable properties of symmetric-conjugate
splitting methods comprising complex coefficients with those of symmetric
splitting methods. Essential ingredients are series expansions that
characterise local errors and the spectral theorem. The construction of
higher-order schemes by composition and numerical illustrations are described
in Section 4.

_Restrictions and generalisations._ We point out that the following analysis
based on infinite series expansions is powerful regarding the treatment of
high-order splitting methods, but it also has some restrictions.

Our main conclusions concerning the accumulation of inaccurate imaginary parts
rely on the assumption that the defining operators correspond to real
symmetric matrices.

To a certain extent, our setting is associated with spatial semi-
discretisations of partial differential equations, but it should be noted that
the constants $C_{A},C_{B}$ in (11) below increase when the space grids are
refined. Furthermore, for parabolic evolution equations, the inverses of
$\mathcal{E}$ and $\mathcal{S}$ in (12) below are not well-defined, since they
involve evaluations at negative times, e.g.

$\big{(}\mathcal{E}(t)\big{)}^{-1}=\mathcal{E}(-\,t)\,,\quad
t\in(0,T-t_{0}]\,.$

We indicate this issue for the simplest splitting method, the Lie–Trotter
splitting method. In situations, where $A$ represents the Laplacian and $B$ a
multiplication operator related to a real-valued potential, the numerical
evolution operator and its time derivative

$\mathcal{S}(t)=\mathrm{e}^{\,tA}\,\mathrm{e}^{\,tB}\,,\quad\tfrac{\mathrm{d}}{\mathrm{d}t}\,\mathcal{S}(t)=\mathrm{e}^{\,tA}\,(A+B)\,\mathrm{e}^{\,tB}\,,\quad
t\in(0,T-t_{0}]\,,$

are well-defined. However, the inverse

$\big{(}\mathcal{S}(t)\big{)}^{-1}=\mathrm{e}^{-\,tB}\,\mathrm{e}^{-\,tA}\,,\quad
t\in(0,T-t_{0}]\,,$

and likewise the composition

$G(t)=\tfrac{\mathrm{d}}{\mathrm{d}t}\,\mathcal{S}(t)\,\big{(}\mathcal{S}(t)\big{)}^{-1}=\mathrm{e}^{\,tA}\,(A+B)\,\mathrm{e}^{-\,tA}\,,\quad
t\in(0,T-t_{0}]\,,$

involve negative times and hence necessitate the restriction to finite
dimensional cases, which result from spatial semi-discretisations.

A rigorous analysis in the lines of [49, 50] is intended for future work. It
will be based on suitable adaptations of the arguments and will require
alternative solution representations as well as specifications of the employed
expansions and arising remainders. Fundamental means for stepwise expansions
of the exact and numerical evolution operators will be provided by the
variation-of-constants formula and Taylor series expansions. The
characterisation of the resulting regularity requirements will be linked to
the identification of iterated commutators.

_Simplified setting of matrices._ We consider the initial value problem for a
linear ordinary differential equation

$\begin{cases}u^{\prime}(t)=F\big{(}u(t)\big{)}=A\,u(t)+B\,u(t)\,,\quad
t\in[t_{0},T]\,,\\\ u(t_{0})\in\mathbb{R}^{M}\text{ given}\,,\end{cases}$
(11a) under the additional assumption that the defining matrices are real and
symmetric $\begin{gathered}A\in\mathbb{R}^{M\times M}\,,\quad\|A\|\leq
C_{A}\,,\quad A^{*}=A^{T}=A\,,\\\ B\in\mathbb{R}^{M\times M}\,,\quad\|B\|\leq
C_{B}\,,\quad B^{*}=B^{T}=B\,.\end{gathered}$ (11b)

The exact and numerical evolution operators

$\begin{gathered}\mathcal{E}_{t}=\mathcal{E}(t)=\mathrm{e}^{\,t\,(A+B)}\in\mathbb{R}^{M\times
M}\,,\\\
\mathcal{S}_{t}=\mathcal{S}(t)=\mathrm{e}^{\,b_{s}\,t\,B}\,\mathrm{e}^{\,a_{s}\,t\,A}\cdots\,\mathrm{e}^{\,b_{1}\,t\,B}\,\mathrm{e}^{\,a_{1}\,t\,A}\in\mathbb{C}^{M\times
M}\,,\\\ t\in[0,T-t_{0}]\,,\end{gathered}$ (12)

are given by exponential series, that is

$\mathrm{e}^{\,t\,L}=\sum_{\ell=0}^{\infty}\tfrac{1}{\ell!}\,t^{\ell}L^{\ell}\,,\quad\big{\|}\mathrm{e}^{\,t\,L}\big{\|}\leq\mathrm{e}^{\,|t|\,\|L\|}\,,\quad
t\in\mathbb{R}\,.$

We note that a splitting method involving complex coefficients yields complex-
valued approximations to the real-valued solution and recall the consistency
condition (2e) ensuring order $p\in\mathbb{N}_{\geq 1}$. For the ease of
notation, we only indicate the dependence on the time increment.

_Series representations._ In the context of matrices, we may use formal
infinite series expansions of decisive components to deduce substantive
results for splitting methods with complex coefficients.

On the one hand, the exact evolution operator

$\mathcal{E}(t)=\mathrm{e}^{\,t\,(A+B)}\,,\quad t\in[0,T-t_{0}]\,,$

satisfies the initial value problem

$\begin{cases}\mathcal{E}^{\prime}(t)=(A+B)\,\mathcal{E}(t)\,,\quad
t\in[0,T-t_{0}]\,,\\\ \mathcal{E}(0)=I\,,\end{cases}$

where $I\in\mathbb{R}^{M\times M}$ denotes the identity matrix.

On the other hand, we make use of the fact that the evolution operator
associated with a splitting method is given as the exponential of a time-
dependent operator

$\mathcal{S}(t)=\mathrm{e}^{\,\Omega(t)}=\mathrm{e}^{\,t\,H(t)}\,,\quad
t\in[0,T-t_{0}]\,,$ (13)

and fulfills a related nonautonomous linear differential equation

$\begin{cases}\mathcal{S}^{\prime}(t)=G(t)\,\mathcal{S}(t)\,,\quad
t\in[0,T-t_{0}]\,,\\\ \mathcal{S}(0)=I\,,\end{cases}$

where formally $G=\mathcal{S}^{\prime}\,\mathcal{S}^{-1}$ as well as
$\Omega^{\prime}=(\exp(\text{ad}_{\Omega})-I)^{-1}\,\text{ad}_{\Omega}\,G$.
Important findings are that statements on the difference $H-(A+B)$ allow to
draw conclusions on $G-(A+B)$ and hence on $\mathcal{S}-\mathcal{E}$.

More precisely, formal representations of $H-(A+B)$ as infinite series can be
found by applying recursively the Baker–Campbell–Hausdorff formula to the
numerical evolution operator (12). Considering the Lie algebra
$\mathcal{L}(A,B)$ generated by $\\{A,B\\}$ with the commutator

$[A,B]=A\,B-BA$

as Lie bracket and denoting by $\mathcal{L}_{\ell}(A,B)$ the homogeneous
subspace of degree $\ell\in\mathbb{N}_{\geq 2}$ with $k$th basis element
$E_{\ell k}(A,B)$ for $k\in K_{\ell}$, e.g.

$\begin{gathered}E_{21}(A,B)=[A,B]\,,\quad K_{2}=\\{1\\}\,,\\\
E_{31}(A,B)=\big{[}A,[A,B]\big{]}\,,\quad
E_{32}(A,B)=\big{[}B,[A,B]\big{]}\,,\quad K_{3}=\\{1,2\\}\,.\end{gathered}$

we obtain the formal series expansion

$H(t)=A+B+\sum_{\ell=2}^{\infty}t^{\ell-1}\sum_{k\in K_{\ell}}e_{\ell
k}(a,b)\,E_{\ell k}(A,B)\,,\quad t\in[0,T-t_{0}]\,,$ (14)

where $e_{\ell k}(a,b)$ represents a polynomial of degree $k$ with respect to
the complex coefficients $(a_{j},b_{j})_{j=1}^{s}$.

_Particular structures of symmetric-conjugate schemes._ The above formal
representation (14) is valid for arbitrary real matrices
$A,B\in\mathbb{R}^{M\times M}$ and splitting methods with complex
coefficients. Provided that $A$ and $B$ are symmetric and the considered
splitting methods are symmetric-conjugate, it turns out that $H$ has a
particular structure.

We henceforth fix the time increment $h=\tfrac{T-t_{0}}{N}$ for some positive
integer $N\in\mathbb{N}_{\geq 1}$ and use again the convenient notation
$\mathcal{E}_{h}=\mathcal{E}(h)$, $\mathcal{S}_{h}=\mathcal{S}(h)$, and
$H_{h}=H(h)$. For brevity, we do not indicate the dependencies on the defining
operators $E_{\ell k}=E_{\ell k}(A,B)$ and the splitting coefficients $e_{\ell
k}=e_{\ell k}(a,b)$.

Evidently, the properties real and symmetric are inherited by the exact
evolution operator

$\mathcal{E}_{h}=\mathcal{E}_{h}^{T}\in\mathbb{R}^{M\times M}\,.$

Moreover, for symmetric-conjugate splitting methods, we conclude that

$\mathcal{S}_{h}=\mathcal{S}_{h}^{*}\in\mathbb{C}^{M\times M}\,,\quad
H_{h}=H_{h}^{*}\in\mathbb{C}^{M\times M}\,,$ (15)

see also (3) and (11). Note, however, that this is in general not true for
symmetric splitting methods. Observing that the iterated commutators of real
symmetric matrices satisfy the relations

$E_{\ell k}^{T}=(-1)^{\ell+1}\,E_{\ell k}\,,\quad k\in
K_{\ell}\,,\quad\ell\in\mathbb{N}_{\geq 2}\,,$

e.g. $E_{21}^{T}=[A,B]^{T}=BA-AB=-\,E_{21}$, the difference of the formal
series expansion (14) and its adjoint

$\begin{split}0&=H_{h}-H_{h}^{*}=\sum_{\ell=2}^{\infty}h^{\ell-1}\sum_{k\in
K_{\ell}}\big{(}e_{\ell k}\,E_{\ell k}-\overline{e_{\ell k}}\,E_{\ell
k}^{T}\big{)}\\\ &=\sum_{\ell=2}^{\infty}h^{\ell-1}\sum_{k\in
K_{\ell}}\Big{(}\big{(}1+(-1)^{\ell}\big{)}\,\Re(e_{\ell
k})+\mathrm{i}\,\big{(}1+(-1)^{\ell+1}\big{)}\,\Im(e_{\ell k})\Big{)}\,E_{\ell
k}\,,\end{split}$

imply that the arising coefficients are either real or purely imaginary
quantities

$\begin{cases}e_{\ell k}=r_{\ell k}\,,&\text{$\ell$ odd}\,,\\\ e_{\ell
k}=\mathrm{i}\,r_{\ell k}\,,&\text{$\ell$ even}\,,\end{cases}\quad r_{\ell
k}\in\mathbb{R}\,,\quad k\in K_{\ell}\,,\quad\ell\in\mathbb{N}_{\geq 2}\,.$

In consequence, we obtain a decomposition of the complex matrix $H_{h}$ into
real symmetric and skew-symmetric contributions

$\begin{gathered}H_{h}=A+B+H_{h}^{(\mathbb{R},\text{sym})}+\mathrm{i}\,H_{h}^{(\mathbb{R},\text{skew})}\,,\\\
H_{h}^{(\mathbb{R},\text{sym})}=\underset{\ell\text{
odd}}{\sum_{\ell=3}^{\infty}}h^{\ell-1}\sum_{k\in K_{\ell}}r_{\ell k}\,E_{\ell
k}\in\mathbb{R}^{M\times M}\,,\\\
H_{h}^{(\mathbb{R},\text{sym})}=\big{(}H_{h}^{(\mathbb{R},\text{sym})}\big{)}^{T}\,,\\\
H_{h}^{(\mathbb{R},\text{skew})}=\underset{\ell\text{
even}}{\sum_{\ell=2}^{\infty}}h^{\ell-1}\sum_{k\in K_{\ell}}r_{\ell
k}\,E_{\ell k}\in\mathbb{R}^{M\times M}\,,\\\
H_{h}^{(\mathbb{R},\text{skew})}=-\,\big{(}H_{h}^{(\mathbb{R},\text{skew})}\big{)}^{T}\,.\end{gathered}$
(16)

_Errors in imaginary parts._ The subsequent considerations will explain the
favourable behaviour of symmetric-conjugate splitting methods in comparison
with symmetric splitting methods for linear ordinary differential equations
that are defined by real symmetric matrices. We first restate the above
decomposition (16) in terms of local errors and next analyse the accumulation
of errors over time.

  1. (i)

A splitting method has order $p\in\mathbb{N}_{\geq 1}$, i.e.
$\mathcal{S}_{h}-\mathcal{E}=\mathcal{O}(h^{p+1})$, if
$H_{h}-(A+B)=\mathcal{O}(h^{p})$. Specifically, if a symmetric-conjugate
operator is of even order $p\in\mathbb{N}_{\geq 2}$, then

$H_{h}^{(\mathbb{R},\text{sym})}=\mathcal{O}\big{(}h^{p}\big{)}\,,\quad
H_{h}^{(\mathbb{R},\text{skew})}=\mathcal{O}\big{(}h^{p+1}\big{)}\,.$

Otherwise, if it is of odd order $p\in\mathbb{N}_{\geq 1}$, then

$H_{h}^{(\mathbb{R},\text{sym})}=\mathcal{O}\big{(}h^{p+1}\big{)}\,,\quad
H_{h}^{(\mathbb{R},\text{skew})}=\mathcal{O}\big{(}h^{p}\big{)}\,.$

  2. (ii)

By means of the spectral theorem applied to the self-adjoint matrix
$H_{h}=H_{h}^{*}\in\mathbb{C}^{M\times M}$, see (15), we conclude that there
exist a unitary matrix $U_{h}\in\mathbb{C}^{M\times M}$ and a real diagonal
matrix $D_{h}\in\mathbb{R}^{M\times M}$ such that

$H_{h}=U_{h}\,D_{h}\,U_{h}^{*}\,.$

Notice that evaluation at $h=0$ implies that the corresponding transformation
matrix is real

$A+B=H_{0}=U_{0}\,D_{0}\,U_{0}^{T}\,,\quad U_{0}\in\mathbb{R}^{M\times M}\,.$

Due to $H_{h}=H_{0}+\mathcal{O}(h^{p})$, the matrices depend smoothly on the
time increment

$U_{h}=U_{0}+\mathcal{O}\big{(}h^{p}\big{)}\,,\quad
D_{h}=D_{0}+\mathcal{O}\big{(}h^{p}\big{)}\,.$

In particular, the errors in the imaginary parts fulfill

$\Im(U_{h})=\mathcal{O}\big{(}h^{p}\big{)}\,.$

In consequence, we obtain the following identities for the exact and numerical
evolution operators and multiple compositions

$\begin{gathered}\mathcal{E}_{h}=\mathrm{e}^{\,h\,(A+B)}=U_{0}\,\mathrm{e}^{\,h\,D_{0}}\,U_{0}^{T}\,,\\\
\mathcal{E}_{h}^{n}=\big{(}U_{0}\,\mathrm{e}^{\,h\,D_{0}}\,U_{0}^{T}\big{)}^{n}=U_{0}\,\mathrm{e}^{\,n\,h\,D_{0}}\,U_{0}^{T}\,,\\\
\mathcal{S}_{h}=\mathrm{e}^{\,h\,H_{h}}=U_{h}\,\mathrm{e}^{\,h\,D_{h}}\,U_{h}^{*}\,,\\\
\mathcal{S}_{h}^{n}=\big{(}U_{h}\,\mathrm{e}^{\,h\,D_{h}}\,U_{h}^{*}\big{)}^{n}=U_{h}\,\mathrm{e}^{\,n\,h\,D_{h}}\,U_{h}^{*}\,,\\\
\quad n\in\\{0,1,\dots,N\\}\,.\end{gathered}$ (17)

As $\mathrm{e}^{\,n\,h\,D_{h}}$ is bounded by $\mathrm{e}^{\,T\,\|D_{h}\|}$,
we finally conclude that the relation

$\mathcal{S}_{h}^{n}=U_{0}\,\mathrm{e}^{\,n\,h\,D_{h}}\,U_{0}^{T}+\mathcal{O}\big{(}h^{p}\big{)}\,,\quad
n\in\\{0,1,\dots,N\\}\,,$ (18)

is valid, where the implicit constant in the $\mathcal{O}$ term does not
depend on $n$.

_Conclusions._ The identity in (18) has several remarkable consequences, which
we confirm and complement by numerical illustrations for model problems with
real-valued solutions in Section 4.

  1. (i)

When applied to the initial value problem

$\begin{cases}u^{\prime}(t)=(A+B)\,u(t)\,,\quad t\in[t_{0},T]\,,\\\
u(t_{0})=u_{0}\,,\end{cases}$ (19)

with time increment $h=\tfrac{T-t_{0}}{N}$, any symmetric-conjugate splitting
method of order $p\in\mathbb{N}_{\geq 1}$ is conjugate to the exact solution
to the initial value problem

$\begin{cases}v^{\prime}(t)=D_{h}\,v(t)\,,\quad t\in[t_{0},T]\,,\\\
v(t_{0})=u_{0}\,,\end{cases}$

where the real diagonal matrix $D_{h}=D_{0}+\mathcal{O}(h^{p})$ is a
perturbation of the same order of the matrix diagonalising the real symmetric
matrix $A+B$, see (3) and (11).

  2. (ii)

The numerical approximation to the real-valued solution has an imaginary part
of the same order $p$, when $p$ is odd, or of order $p+1$, when $p$ is even,
respectively, see also Table 2.

  3. (iii)

This error does not accumulate, and, hence, it does not affect the global
performance. More precisely, for symmetric-conjugate splitting methods, the
relative errors in the imaginary part

$\frac{\|\Im(u_{n})\|}{\|u_{n}\|}\,,\quad
u_{n}=\mathcal{S}_{h}^{n}u_{0}\,,\quad n\in\\{0,1,\dots,N\\}\,,$

remain bounded over time, since they are only due to the transformation by
$U_{h}=U_{0}+\mathcal{O}(h^{p})$. For a comparison of fourth-order symmetric-
conjugate versus symmetric splitting methods, we refer to Figure 2. The
considered situation is related to a ground state computation by the imaginary
time propagation for a linear Schrödinger equation under a quartic potential.
We mirror the errors in the imaginary parts and the ground state energy.

Altogether, we conclude that symmetric-conjugate splitting methods are
particularly favourable for the numerical approximation of linear ordinary
differential equations that are defined by real symmetric matrices (11). The
incorporation of complex coefficients requires the use of complex arithmetics
and increases the computational effort, but it permits the design of high-
order schemes, that satisfy the stability condition (3) and thus overcome the
second-order and fourth-order barriers for standard and modified splitting
methods. Illustrative numerical results are given in Section 4.

##  4\. Test equations and numerical comparisons

In this section, we provide numerical evidence confirming and complementing
our theoretical analysis of symmetric-conjugate splitting methods for spatial
semi-discretisations of parabolic equations.

A Matlab code, which illustrates the practical implementation of operator
splitting methods combined with Fourier spectral space discretisations for
three-dimensional model problems is available at

doi.org/10.5281/zenodo.8238819.

It in particular reproduces numerical results presented in the sequel.

In connection with our first illustration and the numerical computation of
global errors, we focus on the case of a single space dimension. However, for
schemes comprising negative coefficients, severe stability issues have to be
expected when the spatial grid width is refined or in higher dimensions, since
then the problems become significantly stiffer.

_Symmetric-conjugate versus symmetric methods._ A first numerical illustration
is related to the relevant issue of ground state computations for quantum-
mechanical systems. The provided comparisons for different splitting methods
of order four verify and complement our observations in Section 3. In
particular, they show that symmetric-conjugate splitting methods possess
distinctive features for evolution equations that are defined by self-adjoint
operators and have real-valued solutions. Contrary to symmetric counterparts,
it is ensured that the numerical evolution operators inherit the property of
self-adjointness, which results in favourable approximations and is also
reflected in the errors in the imaginary parts and the ground state energies.

_Test equation._ We consider the linear Schrödinger equation (7) for a quartic
potential and a localised Gaussian-like initial state. Moreover, we impose
periodic boundary conditions on a sufficiently large space interval such that
the effect of the truncation error is insignificant. A related parabolic
problem is obtained by integration in imaginary time, i.e. by replacing the
time variable with $-\,\mathrm{i}\,t$, see also (6) and (8a). For the
convenience of the readers, we restate the test equation and specify the data
of the problem

$\begin{gathered}\begin{cases}&\partial_{t}U(x,t)=\tfrac{1}{2}\,\Delta\,U(x,t)-V(x)\,U(x,t)\,,\\\
&U(x,t_{0})=U_{0}(x)\,,\quad(x,t)\in\Omega\times[0,T]\,.\end{cases}\\\
V:\mathbb{R}\longrightarrow\mathbb{R}:x\longmapsto
5-\tfrac{1}{2}\,x^{2}+\tfrac{1}{80}\,x^{4}\,,\\\
U_{0}:\mathbb{R}\longrightarrow\mathbb{R}:x\longmapsto\tfrac{1}{\sqrt[4]{\pi}}\,\mathrm{e}^{-\frac{1}{2}(x-1)^{2}}\,,\\\
\Omega=[-\,a,a]\,,\quad a=10\,,\quad M=256\,.\end{gathered}$ (20)

Within our setting, the eigenvalues of the discrete Laplace operator are non-
positive, see (9). Moreover, due to the fact that the potential takes non-
negative values, it is ensured that the discretisation of the defining
operator $\frac{1}{2}\,\Delta-V$ has negative eigenvalues. The time
propagation of (20) combined with suitable projection thus yields
approximations to stationary states of the quantum-mechanical system. The
ground state is linked to the lowest energy level in modulus, and excited
states correspond to higher energies in modulus.

Basically, the spatially discretised system is defined by real and symmetric
matrices $A,B\in\mathbb{R}^{M\times M}$, see (11) and (19), respectively, and
hence the following considerations are appropriate. The eigenvalues of $A+B$
are negative numbers and the solution value at the final time can formally be
written as a linear combination of associated normalised eigenvectors

$\begin{gathered}(A+B)\,v_{m}=E_{m}\,v_{m}\,,\\\
E_{m}\in\mathbb{R}_{<0}\,,\quad
v_{m}\in\mathbb{R}^{M}\,,\quad\|v_{m}\|=1\,,\quad
m\in\\{0,1,\dots,M-1\\}\,,\\\
u(T)=\sum_{m=0}^{M-1}c_{m}\,\mathrm{e}^{\,T\,E_{m}}\,v_{m}\,.\end{gathered}$

Provided that the dominant eigenvalue is simple, i.e. $E_{m}<E_{0}$ for
$m\in\\{1,\dots,M-1\\}$, and the corresponding coefficient nonzero $c_{0}\neq
0$ such that

$u(T)=c_{0}\,\mathrm{e}^{\,T\,E_{0}}\bigg{(}v_{0}+\sum_{m=1}^{M-1}\tfrac{c_{m}}{c_{0}}\,\mathrm{e}^{-\,T\,(E_{0}-E_{m})}\,v_{m}\bigg{)}\approx
c_{0}\,\mathrm{e}^{\,T\,E_{0}}\,v_{0}\,,$

this allows to determine a numerical approximation to the first eigenvector
related to the ground state

$\tfrac{1}{\|u(T)\|}\,u(T)\approx v_{0}\,.$

We point out that an appropriate choice of the final time has to be adjusted
to the location of the eigenvalues, i.e. the smaller the difference
$E_{0}-E_{1}$, the larger $T>0$. The computation of the ground state energy
then relies on the identity

$E_{0}=v_{0}^{T}(A+B)\,v_{0}\,.$

The space discretisation and in particular the features of the time
discretisation method will affect the quality of the obtained numerical
results.

_Fourth-order splitting methods._ For the time integration of (20), we apply
symmetric and symmetric-conjugate splitting methods of order four,
respectively, see Figures 1 and 9.

For the purpose of illustration, we comment on the construction of the fourth-
order schemes comprising four stages from the second-order Strang splitting
method

$\mathcal{S}^{(F)}_{h}=\mathcal{S}^{[2,\mathbb{R}]}_{h}=h\,\big{(}\tfrac{1}{2},1,\tfrac{1}{2}\big{)}\,,$

by means of the composition technique, see also (2) and (5). In order to
distinguish different schemes, we adapt our former notation and indicate the
numbers of stages as well as structural characteristics, but omit the defining
function. Specifically, we use the triple jump composition

$\begin{split}\mathcal{S}^{[4]}_{h}&=\mathcal{S}_{\alpha_{3}\,h}^{[2]}\circ\mathcal{S}_{\alpha_{2}\,h}^{[2]}\circ\mathcal{S}_{\alpha_{1}\,h}^{[2]}\\\
&=h\,\big{(}\tfrac{1}{2}\,\alpha_{3},\alpha_{3},\tfrac{1}{2}\,(\alpha_{2}+\alpha_{3}),\alpha_{2},\tfrac{1}{2}\,(\alpha_{2}+\alpha_{1}),\alpha_{1},\tfrac{1}{2}\,\alpha_{1}\big{)}\,.\end{split}$

Provided that the arising coefficients fulfil the conditions

$\sum_{j=1}^{3}\alpha_{j}=1\,,\quad\sum_{j=1}^{3}\alpha_{j}^{3}=0\,,$ (21)

this leads to splitting methods of nonstiff order four. These algebraic
equations have the real-valued solutions

$\begin{gathered}\gamma=2-2^{\frac{1}{3}}\,,\quad\alpha_{1}=\gamma^{-1}\,,\quad\alpha_{2}=1-2\,\alpha_{1}\,,\quad\alpha_{3}=\alpha_{1}\,,\\\
\mathcal{S}^{[4,\mathbb{R},\text{sym}]}_{h}=h\,\big{(}\tfrac{\alpha_{1}}{2},\alpha_{1},\tfrac{1-\alpha_{1}}{2},1-2\,\alpha_{1},\tfrac{1-\alpha_{1}}{2},\alpha_{1},\tfrac{\alpha_{1}}{2}\big{)}\,,\end{gathered}$
(22a) which define the standard symmetric Yoshida splitting method, see [51].
In addition, they admit a complex-valued solution, which correspond to the
symmetric splitting method
$\begin{gathered}\gamma=2-2^{\frac{1}{3}}\,\mathrm{e}^{\,\frac{2}{3}\,\mathrm{i}\,\pi}\,,\quad\alpha_{1}=\gamma^{-1}\,,\quad\alpha_{2}=1-2\,\alpha_{1}\,,\quad\alpha_{3}=\alpha_{1}\,,\\\
\mathcal{S}^{[4,\mathbb{C},\text{sym}]}_{h}=h\,\big{(}\tfrac{\alpha_{1}}{2},\alpha_{1},\tfrac{1-\alpha_{1}}{2},1-2\,\alpha_{1},\tfrac{1-\alpha_{1}}{2},\alpha_{1},\tfrac{\alpha_{1}}{2}\big{)}\,,\end{gathered}$
(22b) as well as a complex-valued solution, which correspond to the symmetric-
conjugate splitting method
$\begin{gathered}\alpha_{1}=\tfrac{1}{4}+\mathrm{i}\,\tfrac{\sqrt{15}}{12}\,,\quad\alpha_{2}=\tfrac{1}{2}\,,\quad\alpha_{3}=\overline{\alpha_{1}}\,,\\\
\mathcal{S}^{[4,\mathbb{C},\text{sym-
conj}]}_{h}=h\,\big{(}\tfrac{\overline{\alpha_{1}}}{2},\overline{\alpha_{1}},\tfrac{\overline{\alpha_{1}}+\alpha_{2}}{2},\alpha_{2},\tfrac{\alpha_{1}+\alpha_{2}}{2},\alpha_{1},\tfrac{\alpha_{1}}{2}\big{)}\,.\end{gathered}$
(22c)

Further schemes are obtained from (21) by complex conjugation. The symmetric
and symmetric-conjugate splitting methods have in common that they contain the
same number of exponentials and provide fourth-order approximations. Moreover,
the sizes of the main error terms at order five, measured as

$\text{err}=\bigg{|}\sum_{j=1}^{3}\alpha_{j}^{5}\bigg{|}\,,$

are about 200 times smaller than the error of the triple jump composition with
real coefficients

$\text{err}^{[4,\mathbb{R},\text{sym}]}_{h}\approx
5.3\,,\quad\text{err}^{[4,\mathbb{C},\text{sym}]}_{h}\approx
0.024\,,\quad\text{err}^{[4,\mathbb{C},\text{sym-conj}]}_{h}\approx 0.028\,,$

see also [14].

The additionally considered symmetric and symmetric-conjugate splitting
methods of order four comprise a higher number of stages $s>4$ and include
positive coefficients $(a_{j})_{j=1}^{s}$, which distinguishes them from the
schemes with four stages and ensures stability for parabolic problems as well
as Schrödinger equations. The selected optimised symmetric-conjugate splitting
method shows a favourable accuracy behaviour and leads to relatively small
errors.

_Numerical results._ In order to confirm the expected qualitative and
quantitative differences between the above described fourth-order symmetric-
conjugate and symmetric splitting methods, we prescribe a certain time
increment $h>0$ and perform the integration until a sufficiently large final
time $T=Nh$ is reached. In addition, a reference solution
$u(T)\in\mathbb{R}^{M\times M}$ with real values up to machine precision and
the associated ground state energy $E_{0}=u(T)^{T}(A+B)\,u(T)$ are computed
numerically. Then, at each time step, the relative sizes of the imaginary
parts with respect to the solution values and the relative errors of the
ground state energies

$\begin{gathered}\frac{\|\Im(u_{n})\|}{\|u_{n}\|}\,,\\\
\frac{|E_{0}-\Re(u_{n})^{T}(A+B)\,\Re(u_{n})|}{|E_{0}|}\,,\\\
n\in\\{0,1,\dots,N\\}\,,\end{gathered}$

are determined.

As explained in Section 3, for the symmetric-conjugate splitting methods, it
is expected that the errors in the imaginary parts remain bounded, whereas for
symmetric splitting methods involving complex coefficients a significant error
growth over time may occur. Concerning the errors in the approximation of the
ground state energy, the quantities $\Re(u_{n})^{T}(A+B)\,\Re(u_{n})$ for
$n\in\\{0,1,\dots,N\\}$ constitute approximations to the ground state of a
perturbed matrix that depends on the splitting method and the time increment.
Hence, it is expected that the errors decrease up to a certain time, which
depends on the order of the method and the time increment. Beyond that time,
it is again assumed that the error accumulation in the imaginary part will
lead to a significant error growth for symmetric splitting methods and bounded
errors for symmetric-conjugate schemes.

The numerical results displayed in Figure 2 illustrate how the different
characters of symmetric-conjugate versus symmetric splitting methods manifest
in practice. In the left panel, we depict the relative errors in the imaginary
parts of the numerical solutions over time. In the right panel, we display the
corresponding relative errors in the ground state energy. We indeed observe
that the symmetric schemes introduce errors that grow linearly in a log-log
scale, which may lead to unphysical effects. Contrary, after a transient time,
the errors committed by the symmetric-conjugate schemes are nearly constant.

_Local and global errors._

In the lines of our first illustration, we next consider the linear test
equation (8a) subject to periodic boundary conditions on a suitably chosen
space interval and apply the operator splitting methods listed in Figure 1.
Specifically, we set $\alpha=1=\beta$ and study the quadratic and quartic
potentials $V:\mathbb{R}\longrightarrow\mathbb{R}:x\longmapsto x^{2}\,,\quad
V:\mathbb{R}\longrightarrow\mathbb{R}:x\longmapsto\tfrac{1}{24}\,x^{4}\,.$
(23a) The space grid points and the time interval are given by
$\begin{gathered}\Omega=[-\,a,a]\,,\quad a=10\,,\quad M=100\,,\quad
t_{0}=0\,,\quad T\in\\{1,10\\}\,.\end{gathered}$ (23b) For the problem, when
considered on the whole real line with quadratic potential and initial state
$U_{0}:\mathbb{R}\longrightarrow\mathbb{R}:x\longmapsto\mathrm{e}^{-\frac{1}{2}\,x^{2}}\,,$
(23c) the exact solution is known. In this special situation, we thus have the
possibility to verify the correctness of the implementation and conclude that
the errors caused by the truncation of the space domain, the implicitely
imposed periodic boundary conditions, and the application of the Fourier
spectral discretisation method are insignificant. In the general case, we
determine the local and global time discretisation errors with respect to
numerical reference solutions.

Altogether, the results displayed in Figures 3 and 4 show that the nonstiff
orders of convergence are retained. For the errors in the imaginary parts of
symmetric-conjugate schemes of even orders, we indeed observe
superconvergence, see Table 2 and Section 3.

_Adaptive time integration._ The design and theoretical analysis of an
automatic stepsize control algorithm for the time integration of evolution
equations is a complex subject, see for instance [28, 30, 45]. In order to
demonstrate its practical usefulness, we reconsider the linear test equation
(8a) for the above specified situation, where the exact solution is available.
A particular benefit of higher-order symmetric-conjugate splitting methods,
which distinguishes them from symmetric counterparts, is the possibility to
base the local error estimation on the computation of the imaginary parts with
negligible additional costs, see Section 3. More precisely, for the solution
value at the current time $t_{n}$, we determine

$\text{local error
estimate}=\big{\|}\Im\big{(}u(t_{n})\big{)}\big{\|}_{q}\,,\quad
q\in\\{2,\infty\\}\,.$

As standard, for a prescribed tolerance, the optimal time stepsize is then
adjusted, in essence, through

$\tau_{\text{optimal}}=\tau\,\sqrt[p+1]{\frac{\text{tolerance}}{\text{local
error estimate}}}\,.$

The total number of time steps depends on the choice of the norm, since

$\big{\|}\Im\big{(}u(t_{n})\big{)}\big{\|}_{\infty}\leq\big{\|}\Im\big{(}u(t_{n})\big{)}\big{\|}_{2}\,.$

Our tests confirm this dependence, yielding smaller time stepsizes for the
Euclidean norm ($q=2$) compared to the maximum norm ($q=\infty$). The results
obtained for third- and sixth-order schemes with respect to the Euclidean norm
are displayed in Figures 5 and 6. On the one hand, it is observed that the
third-order symmetric-conjugate splitting method performs 47 and 997 time
steps to fullfill the requirements that the local errors remain below the
prescribed tolerances $10^{-6}$ and $10^{-10}$, respectively. On the other
hand, for the sixth-order symmetric-conjugate splitting method and the
prescribed tolerances $10^{-10}$ and $10^{-12}$, significantly reduced numbers
of time steps, namely 6 and 14, are performed. The corresponding results with
respect to the maximum norm are displayed in Figures 7 and 8. Overall, we
observe good correlations of the prescribed tolerances for the local errors
and the resulting global errors.

##  5\. Conclusions and future investigations

The present work is dedicated to a comprehensive analysis of symmetric-
conjugate operator splitting methods for the time integration of linear
evolution equations. It is seen that the natural approach to incorporate
complex coefficients with non-negative real parts permits the design of high-
order schemes that remain stable in the context of parabolic problems and
thereby overcome the order barriers for standard and modified splitting
methods with real coefficients.

Moreover, it is demonstrated that symmetric-conjugate splitting methods are
particularly favourable in the numerical integration of nonreversible systems
defined by real and symmetric matrices. The main reasons are that the errors
in the imaginary parts and energies remain bounded and hence do not lead to
unphysical perturbations. Typically, this kind of problems arises in ground
and excited state computations for Schrödinger equations by the imaginary time
propagation, fractal path integrals with applications to many-body theories
and statistical physics as well as Monte Carlo simulations of quantum systems.

Future theoretical and numerical investigations will concern extensions to
nonlinear evolution equations and a rigorous convergence analysis. Special
attention will be given to complex operator splitting methods applied to
complex Ginzburg–Landau equations, since this relevant type of problems
interlinks certain characteristics of parabolic as well as Schrödinger-type
equations.

## Acknowledgements

Part of this work was developed during a research stay at the Wolfgang Pauli
Institute Vienna; the authors are grateful to the director Norbert Mauser and
the staff members for their support and hospitality. This work has been
supported by Ministerio de Ciencia e Innovación (Spain) through projects
PID2019-104927GB-C21 and PID2019-104927GB-C22, MCIN/AEI/10.13039/501100011033,
ERDF (_A way of making Europe_). Sergio Blanes and Fernando Casas acknowledge
the support of the Conselleria d’Innovació, Universitats, Ciència i Societat
Digital from the Generalitat Valenciana (Spain) through project
CIAICO/2021/180.

## References

  * [1] W. Auzinger, H. Hofstätter, O. Koch. Non-existence of generalized splitting methods with positive coefficients of order higher than four. Appl. Math. Lett. 97 (2019) 48–52.
  * [2] P. Bader, S. Blanes, F. Casas. Solving the Schrödinger eigenvalue problem by the imaginary time propagation technique using splitting methods with complex coefficients. J. Chem. Phys. 139 (2013) 124117.
  * [3] A. Bandrauk, H. Shen. Improved exponential split operator method for solving the time-dependent Schrödinger equation. Chem. Phys. Lett. 176 (1991) 428–432.
  * [4] W. Bao. Ground states and dynamics of multicomponent Bose–Einstein condensates. Multiscale Model. Simul. 2 (2004) 210–236.
  * [5] W. Bao, Q. Du. Computing the ground state solution of Bose–Einstein condensates by a normalized gradient flow. SIAM J. Sci. Comput. 25 (2004) 1674–1697.
  * [6] W. Bao, S. Jin, P. Markowich. On time-splitting spectral approximations for the Schrödinger equation in the semiclassical regime. J. Comput. Phys. 175 (2002) 487–524.
  * [7] J. Bernier, S. Blanes, F. Casas, A. Escorihuela-Tomàs. Symmetric-conjugate splitting methods for linear unitary problems. Submitted (2023).
  * [8] G. Bertoli, C. Besse, G. Vilmart. Superconvergence of the Strang splitting when using the Crank–Nicolson scheme for parabolic PDEs with Dirichlet and oblique boundary conditions. Math. Comp. 90 (2021) 2705–2729.
  * [9] G. Bertoli, G. Vilmart. Strang splitting method for semilinear parabolic problems with inhomogeneous boundary conditions: a correction based on the flow of the nonlinearity. SIAM J. Sci. Comput. 42 (2020), A1913–A1934.
  * [10] S. Blanes, F. Casas. On the necessity of negative coefficients for operator splitting schemes of order higher than two. Appl. Numer. Math. 54 (2005) 23–37.
  * [11] S. Blanes, F. Casas, A. Murua. Splitting methods with complex coefficients. Bol. Soc. Esp. Mat. Apl. 50 (2010) 47–61.
  * [12] S. Blanes, F. Casas, P. Chartier, A. Murua. Optimized high-order splitting methods for some classes of parabolic equations. Math. Comp. 82 (2013) 1559–1576.
  * [13] S. Blanes, F. Casas, M. Thalhammer. Convergence analysis of high-order commutator-free quasi-Magnus exponential integrators for nonautonomous linear evolution equations of parabolic type. IMA J. Numer. Anal. 38 (2018) 743–778.
  * [14] S. Blanes, F. Casas, A. Escorihuela-Tomàs. Applying splitting methods with complex coefficients to the numerical integration of unitary problems. J. Comput. Dyn. 9 (2022) 85–101.
  * [15] S. Blanes, F. Casas, P. Chartier, A. Escorihuela-Tomàs. On symmetric-conjugate composition methods in the numerical integration of differential equations. Math. Comp. 91 (2022) 1739–1761.
  * [16] S. Blanes, F. Casas, C. González, M. Thalhammer. Efficient splitting methods based on modified potentials: numerical integration of linear parabolic problems and imaginary time propagation of the Schrödinger equation. Commun. Comput. Phys. 33 (2023) 937–961.
  * [17] S. Blanes, F. Casas, C. González, M. Thalhammer. Generalisation of splitting methods based on modified potentials to nonlinear evolution equations of parabolic and Schrödinger type. Comput. Phys. Commun. 295 (2024) 109007.
  * [18] M. Caliari, S. Zuccher. A fast time splitting finite difference approach to Gross–Pitaevskii equations. Commun. Comput. Phys. 29 (2021) 1336–1364.
  * [19] F. Castella, P. Chartier, S. Decombes, G. Vilmart. Splitting methods with complex times for parabolic equations. BIT 49 (2009) 487–508.
  * [20] S. Chin. Symplectic integrators from composite operator factorizations. Phys. Lett. A 226 (1997) 344–348.
  * [21] I. Danaila, B. Protas. Computation of ground states of the Gross–Pitaevskii functional via Riemannian optimization. SIAM J. Sci. Comput. 39 (2017) 1102–1129.
  * [22] S. Descombes, M. Duarte, T. Dumont, F. Laurent, V. Louvet, M. Massot. Analysis of operator splitting in the nonasymptotic regime for nonlinear reaction-diffusion equations. Application to the dynamics of premixed flames. SIAM J. Numer. Anal. 52/3 (2014) 1311–1334.
  * [23] L. Einkemmer, A. Ostermann. Overcoming order reduction in diffusion-reaction splitting. Part 2: Oblique boundary conditions. SIAM J. Sci. Comput. 38/6 (2016) A3741–A3757.
  * [24] K. J. Engel, R. Nagel. One-Parameter Semigroups for Linear Evolution Equations. Springer, 2000.
  * [25] D. Goldman, T. Kaper. $n$th-order operator splitting schemes and nonreversible systems. SIAM J. Numer. Anal. 33 (1996) 349–367.
  * [26] C. González, M. Thalhammer. A second-order Magnus-type integrator for quasi-linear parabolic problems. Math. Comp. 76/257 (2007) 205–231.
  * [27] F. Goth. Higher order auxiliary field quantum Monte Carlo methods. J. Phys. Conf. Ser. 2207 (2022) 012029.
  * [28] K. Gustafsson, M. Lundh, G. Söderlind. API stepsize control for the numerical solution of ordinary differential equations. BIT 28 (1988) 270–-287.
  * [29] E. Hairer, Ch. Lubich, G. Wanner. Geometric Numerical Integration. Springer, 2006.
  * [30] E. Hairer, S. P. Nørset, G. Wanner. Solving Ordinary Differential Equations II. Springer, 2002.
  * [31] E. Hansen, A. Ostermann. Exponential splitting for unbounded operators. Math. Comp. 78/267 (2009) 1485–1496.
  * [32] E. Hansen, A. Ostermann. High order splitting methods for analytic semigroups exist. BIT 49 (2009) 527–542.
  * [33] E. Hille, R. Phillips. Functional Analysis and Semi-Groups. American Mathematical Society, 1974.
  * [34] T. Jahnke, C. Lubich. Error bounds for exponential operator splittings. BIT Numer. Math. 40/4 (2000) 735–744.
  * [35] H. Kleinert. Gauge Fields in Condensed Matter. World Scientific, 1989.
  * [36] E. Kieri. Stiff convergence of force-gradient operator splitting methods. Appl. Numer. Math. 94 (2015) 33–45.
  * [37] R. Kozlov, A. Kværnø, B. Owren. The behaviour of the local error in splitting methods applied to stiff problems. J. Comput. Phys. 195/2 (2004) 576–593.
  * [38] A. Lukassen, M. Kiehl. Operator splitting for chemical reaction systems with fast chemistry. J. Comput. Appl. Math. 344 (2019) 495–511.
  * [39] A. Lunardi. Analytic Semigroups and Optimal Regularity in Parabolic Problems. Birkhäuser/Springer, 1995, reprint 2013.
  * [40] R. McLachlan, G. Quispel. Splitting methods. Acta Numerica 11 (2002) 341–434.
  * [41] A. Messiah. Quantum Mechanics. Dover, 1999.
  * [42] A. Pazy. Semigroups of linear operators and applications to partial differential equations. Springer, 1983.
  * [43] I. Omelyan, I. Mryglod, R. Folk. On the construction of high order force gradient algorithms for integration of motion in classical and quantum systems. Phys. Rev. E 66 (2002) 026701.
  * [44] I. Omelyan, I. Mryglod, R. Folk. Symplectic analytically integrable decomposition algorithms: classification, derivation, and application to molecular dynamics, quantum and celestial mechanics simulations. Comput. Phys. Comm. 151 (2003) 272–314.
  * [45] W. Press, S. Teukolsky, W. Vetterling, B. Flannery. Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press, 2007.
  * [46] J. Sanz-Serna, M. Calvo. Numerical Hamiltonian Problems. Dover, 2018.
  * [47] Q. Sheng. Solving linear partial differential equations by exponential splitting. IMA J. Numer. Anal. 9 (1989) 199–212.
  * [48] M. Suzuki. General theory of fractal path integrals with applications to many-body theories and statistical physics. J. Math. Phys. 32 (1991) 400–407.
  * [49] M. Thalhammer. High-order exponential operator splitting methods for time-dependent Schrödinger equations. SIAM J. Numer. Anal. 46 (2008) 2022–2038.
  * [50] M. Thalhammer. Convergence analysis of high-order time-splitting pseudo-spectral methods for nonlinear Schrödinger equations. SIAM J. Numer. Anal. 50 (2012) 3231–3258.
  * [51] H. Yoshida. Construction of higher order symplectic integrators. Phys. Lett. A 150 (1990) 262–268.

Figure 1. Real and complex splitting methods applied in numerical tests. Denominations and characteristics (nonstiff order $p$, number of stages $s$). The coefficients of the symmetric-conjugate schemes are given in Figures 9 and 10. Lie–Trotter (real, $p=1$, $s=1$) |  Stability ($a_{1}>0$)
---|---
Strang (real, $p=2$, $s=2$) |  Stability ($a_{1},a_{2}\geq 0$)
Yoshida (real, $p=4$, $s=4$) |  Instability ($a_{3}<0$)
Complex (symmetric, $p=4$, $s=4$) |  Stability ($\Re\,a_{1},\dots,\Re\,a_{s}\geq 0$)
Complex (symmetric, $p=4$, $s=5$) |  Stability ($a_{1},\dots,a_{s}\geq 0$)
Complex (symmetric, $p=4$, $s=6$) |  Stability ($a_{1},\dots,a_{s}\geq 0$)
Complex (symmetric, $p=6$, $s=17$) |  Stability ($a_{1},\dots,a_{s}\geq 0$)
Complex (symmetric-conj., $p=3$, $s=3$) |  Stability ($\Re\,a_{1},\dots,\Re\,a_{s}\geq 0$)
Complex (symmetric-conj., $p=3$, $s=4$) |  Stability ($a_{1},\dots,a_{s}\geq 0$)
Complex (symmetric-conj., $p=4$, $s=4$) |  Stability ($\Re\,a_{1},\dots,\Re\,a_{s}\geq 0$)
Complex (symmetric-conj., $p=4$, $s=6$) |  Stability ($a_{1},\dots,a_{s}\geq 0$)
Complex (symmetric-conj., $p=6$, $s=12$) |  Stability ($a_{1},\dots,a_{s}\geq 0$)
Complex (symmetric-conj., $p=6$, $s=16$) |  Stability ($a_{1},\dots,a_{s}\geq 0$)

Table 1. Stability properties of real and complex splitting methods in the context of parabolic equations. Schemes with non-negative coefficients $(a_{j})_{j=1}^{s}$ remain stable for Schrödinger equations. Lie–Trotter (real, $p=1$) | $p_{\,\text{num}}=p=1$ | —
---|---|---
Strang (real, $p=2$) | $p_{\,\text{num}}=p=2$ | —
Yoshida (real, $p=4$) | $p_{\,\text{num}}=p=4$ | —
Complex (symmetric, $p=4$) | $p_{\,\text{num}}=p=4$ | $p_{\,\text{num},\Im}=p=4$
Complex (symmetric, $p=6$) | $p_{\,\text{num}}=p=6$ | $p_{\,\text{num},\Im}=p=6$
Complex (symmetric-conj., $p=3$) | $p_{\,\text{num}}=p=3$ | $p_{\,\text{num},\Im}=p=3$
Complex (symmetric-conj., $p=4$) | $p_{\,\text{num}}=p=4$ | $p_{\,\text{num},\Im}=p+1=5$
Complex (symmetric-conj., $p=6$) | $p_{\,\text{num}}=p=6$ | $p_{\,\text{num},\Im}=p+1=7$

Table 2. Application of real and complex splitting methods to parabolic model
problems with real-valued solutions. List of classical orders (first column),
numerically observed orders of convergence for solution values (second
column), and numerically observed orders for imaginary parts (third column).
Figure 2. Time integration of the parabolic model problem (20) by non-
optimised and optimised fourth-order operator splitting methods involving
complex coefficients with time increments $h=\frac{T}{40}$ (top curves) and
$h=\frac{T}{400}$ (bottom curves). Relative errors in the imaginary parts of
the numerical solutions over time (left) and corresponding errors in the
ground state energy (right). Symmetric schemes comprising $s=4$ (thin black
dashed line) and $s>4$ (thin black solid line) stages with increasing errors
in a log-log scale versus symmetric-conjugate schemes comprising $s=4$ (thick
red dashed line) and $s>4$ (thick red solid line) stages with bounded errors.

Figure 3. Time integration of the linear parabolic model problem with real-
valued solution by real and complex splitting methods, see also (8) and Figure
1. For the considered quadratic potential, the exact solution is known. Left:
Local and global errors. Right: Corresponding errors in the imaginary parts.

Figure 4. Time integration of the linear parabolic model problem with real-
valued solution by real and complex splitting methods, see also (8) and Figure
1. For the considered quartic potential, a numerical reference solution is
computed. Left: Local and global errors. Right: Corresponding errors in the
imaginary parts.

Figure 5. Adaptive time integration of a linear parabolic model problem with
known real-valued solution for $t_{0}=0$ and $T=1$ by a third-order symmetric-
conjugate splitting method, see also (8) and Figure 1. The local error
estimation is based on the computation of the imaginary parts with respect to
the Euclidean norm and requires negligible additional costs. The total numbers
of time steps 47 and 997 are adjusted in accordance with the prescribed
tolerances $10^{-6}$ and $10^{-10}$, respectively. Left: Sequences of time
grid points. Right: Associated sequences of local errors determined with
respect to the exact solution values.

Figure 6. Corresponding results for a sixth-order symmetric-conjugate
splitting method and lower tolerances $10^{-10}$ and $10^{-12}$, respectively.

Figure 7. Corresponding results for the third-order symmetric-conjugate
splitting method with respect to the maximum norm.

Figure 8. Corresponding results for the sixth-order symmetric-conjugate
splitting method with respect to the maximum norm. Figure 9. Coefficients of
symmetric-conjugate operator splitting methods applied in numerical tests.
Figure 10. Coefficients of symmetric-conjugate operator splitting methods
applied in numerical tests.