# Generalisation of splitting methods based
on modified potentials to nonlinear evolution equations of parabolic and
Schrödinger type

S. Blanes, F. Casas, C. González, M. Thalhammer

###### Abstract

The present work is concerned with the extension of modified potential
operator splitting methods to specific classes of nonlinear evolution
equations. The considered partial differential equations of Schrödinger and
parabolic type comprise the Laplacian, a potential acting as multiplication
operator, and a cubic nonlinearity. Moreover, an invariance principle is
deduced that has a significant impact on the efficient realisation of the
resulting modified operator splitting methods for the Schrödinger case.

Numerical illustrations for the time-dependent Gross–Pitaevskii equation in
the physically most relevant case of three space dimensions and for its
parabolic counterpart related to ground state and excited state computations
confirm the benefits of the proposed fourth-order modified operator splitting
method in comparison with standard splitting methods.

The presented results are novel and of particular interest from both, a
theoretical perspective to inspire future investigations of modified operator
splitting methods for other classes of nonlinear evolution equations and a
practical perspective to advance the reliable and efficient simulation of
Gross–Pitaevskii systems in real and imaginary time.

###### keywords:

Nonlinear evolution equations , Parabolic problems , Schrödinger equations ,
Gross–Pitaevskii systems , Geometric time integration , Operator splitting
methods , Fourier spectral method , Stability , Convergence , Efficiency.

111Addresses: Sergio Blanes, Universitat Politècnica de València, Instituto de
Matemática Multidisciplinar, 46022 Valencia, Spain. Fernando Casas,
Universitat Jaume I, IMAC and Departament de Matemàtiques, 12071 Castellón,
Spain. Cesáreo González, Universidad de Valladolid, Departamento de Matemática
Aplicada, 47011 Valladolid, Spain. Mechthild Thalhammer, Leopold–Franzens-
Universität Innsbruck, Institut für Mathematik, 6020 Innsbruck, Austria.
Websites: personales.upv.es/serblaza, www.gicas.uji.es/fernando.html,
www.imuva.uva.es/en/investigadores/33, techmath.uibk.ac.at/mecht. Email
addresses<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Scope of applications. A wide range of relevant applications in sciences
includes the numerical integration of initial value problems for nonlinear
evolution equations. In many cases, the function defining the right-hand side
comprises two or more parts

$\begin{cases}\tfrac{\text{d}}{\text{d}t}\,u(t)=F_{1}\big{(}u(t)\big{)}+F_{2}\big{(}u(t)\big{)}\,,\\\
u(0)=u_{0}\,,\quad t\in[0,T]\,.\end{cases}$ (1)

As prominent instances, we highlight nonlinear Schrödinger equations, more
specifically, time-dependent Gross–Pitaevskii equations that arise in the
description of Bose–Einstein condensation, see [17, 25]. For a comprehensive
overview of the underlying principles of quantum theory, we refer to [24].

Nonlinear Schrödinger equation (Gross–Pitaevskii equation). A fundamental
model for the nonlinear dynamics of a single Bose–Einstein condensate reads as

$\begin{cases}\text{i}\,\partial_{t}\Psi(x,t)=-\,\Delta\Psi(x,t)+V(x)\,\Psi(x,t)+\vartheta\,|\Psi(x,t)|^{2}\,\Psi(x,t)\,,\\\
\Psi(x,0)=\Psi_{0}(x)\,,\quad(x,t)\in\Omega\times[0,T]\,,\end{cases}$ (2a)
where $\Delta=\partial_{x_{1}}^{2}+\dots+\partial_{x_{d}}^{2}$ denotes the
Laplacian with respect to the spatial variables
$x=(x_{1},\dots,x_{d})\in\mathbb{R}^{d}$, $V:\mathbb{R}^{d}\to\mathbb{R}$ a
real-valued potential, $\vartheta\in\mathbb{R}$ the coupling constant, and
$\Psi:\Omega\times[0,T]\subset\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{C}$ the
space-time-dependent complex-valued macroscopic wave function. Assigning for a
regular function $v:\Omega\to\mathbb{C}$ the linear differential and nonlinear
multiplication operators
$\begin{gathered}\big{(}F_{1}(v)\big{)}(x)=c\,\Delta\,v(x)\,,\quad
c=\text{i}\,,\\\
\big{(}F_{2}(v)\big{)}(x)=\bar{c}\,\big{(}V(x)+\vartheta\,|v(x)|^{2}\big{)}\,v(x)\,,\quad\bar{c}=-\,\text{i}\,,\\\
x\in\Omega\,,\end{gathered}$ (2b)

and setting $u(t)=\Psi(\cdot,t)$ for $t\in[0,T]$, we retain the general
formulation (1).

Nonlinear parabolic equation. By analogy to the time-dependent
Gross–Pitaevskii equation (2a), we consider the parabolic problem

$\begin{cases}\partial_{t}U(x,t)=\Delta
U(x,t)+V(x)\,U(x,t)+\vartheta\,|U(x,t)|^{2}\,U(x,t)\,,\\\
U(x,0)=U_{0}(x)\,,\quad(x,t)\in\Omega\times[0,T]\,,\end{cases}$ (3a) for a
real-valued solution
$U:\Omega\times[0,T]\subset\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{R}$.
Accordingly, it corresponds to (2b) with different constant
$\begin{gathered}\big{(}F_{1}(v)\big{)}(x)=c\,\Delta\,v(x)\,,\quad c=1\,,\\\
\big{(}F_{2}(v)\big{)}(x)=\bar{c}\,\big{(}V(x)+\vartheta\,|v(x)|^{2}\big{)}\,v(x)\,,\quad\bar{c}=1\,,\\\
x\in\Omega\,,\end{gathered}$ (3b)

and, setting $u(t)=U(\cdot,t)$ for $t\in[0,T]$, we obtain again the general
form (1). It is noteworthy that the parabolic equation (3a) arises in ground
state and excited state computations, see for instance [5, 16].

Splitting approach. In essence, operator splitting methods rely on the
presumption that the numerical approximation of the subproblems

$\begin{gathered}\tfrac{\text{d}}{\text{d}t}\,u_{1}(t)=F_{1}\big{(}u_{1}(t)\big{)}\,,\quad\tfrac{\text{d}}{\text{d}t}\,u_{2}(t)=F_{2}\big{(}u_{2}(t)\big{)}\,,\end{gathered}$

is significantly simpler compared to the numerical approximation of the
original problem (1). Then, within multiple scopes, for ordinary differential
equations and time-dependent partial differential equations, for linear
problems as well as nonlinear problems, a variety of works has confirmed the
benefits of operator splitting methods regarding desirable features that are
subsumed under the central concepts stability, efficiency, and preservation of
conserved quantities. For general information, we refer to [23, 29]. Specific
studies in the context of Schrödinger equations are given, e.g., in [6, 9, 12,
34].

Alternative approach. In this work, we propose an approach that provides a
favourable alternative to standard operator splitting methods in situations,
where the operator $F_{2}$ and an iterated commutator of $F_{2}$ and $F_{1}$,
given by

$\begin{split}G_{2}(v)&=F_{1}^{\prime\prime}(v)\,F_{2}(v)\,F_{2}(v)+F_{1}^{\prime}(v)\,F_{2}^{\prime}(v)\,F_{2}(v)+F_{2}^{\prime}(v)\,F_{2}^{\prime}(v)\,F_{1}(v)\\\
&\qquad-
F_{2}^{\prime\prime}(v)\,F_{1}(v)\,F_{2}(v)-2\,F_{2}^{\prime}(v)\,F_{1}^{\prime}(v)\,F_{2}(v)\,,\end{split}$
(4)

have a similar structure. As relevant nonlinear partial differential equations
with this property, we identify Schrödinger and parabolic equations such as
(2) and (3) that comprise the Laplacian, a potential acting as multiplication
operator, and a cubic nonlinearity. For specifications concerning (4), we in
particular refer to Sections 3 and 4.

Formal means and objectives. Our educated guess that leads us to modified
operator splitting methods relies on a formal generalisation of the linear
case, which we briefly sketch next and describe in further detail in the
subsequent sections.

  1. (i)

_Linear ordinary differential equations._ The starting point is a linear
ordinary differential equation defined by non-commuting square matrices

$\tfrac{\text{d}}{\text{d}t}\,u(t)=A\,u(t)+B\,u(t)\,,\quad t\in[0,T]\,.$

The corresponding solution value at the final time is given by the matrix
exponential, that is

$u(T)=\text{e}^{T(A+B)}\,u(0)=\Big{(}\text{e}^{\tau(A+B)}\Big{)}^{N}\,u(0)\,,\quad\tau=\tfrac{T}{N}\,,\quad
N\in\mathbb{N}\,.$

Standard splitting methods are built on compositions of the factors
$\text{e}^{a\,\tau A}$ and $\text{e}^{b\,\tau B}$ with suitably chosen
coefficients $a$ and $b$. Beyond that, components of the form

$\text{e}^{b\,\tau
B+c\,\tau^{3}\,[B,[B,A]]}\,,\quad\big{[}B,[B,A]\big{]}=B^{2}A-2\,BAB+AB^{2}\,,$
(5)

with certain coefficients $b$ and $c$ are incorporated in modified potential
operator splitting methods. The underlying idea of this approach is to gain
freedom in the adjustment of the method coefficients and, amongst others, to
overcome an order barrier valid for standard splitting methods.

  2. (ii)

_Linear partial differential equations._ Advantages of this approach become
apparent in the context of the imaginary time integration of linear
Schrödinger equations comprising the Laplacian and a potential. There, the
operator arising in (5) reduces to a multiplication operator, which is defined
by the potential and its gradient.

  3. (iii)

_Nonlinear partial differential equations._ The guide line for the extension
to nonlinear evolution equations (1) is provided by the calculus of Lie
derivatives, see [35] for a detailed exposition. In order to make our
contribution accessible to a broader readership, we do not presume the
knowledge of this formal calculus and explain the required elementary means on
occasion. Basically, the operators $F_{1}$ and $F_{2}$ take the roles of the
matrices $A$ and $B$, and (5) is replaced by the solution to

$\tfrac{\text{d}}{\text{d}t}\,u(t)=b\,F_{2}\big{(}u(t)\big{)}+c\,\tau^{2}\,G_{2}\big{(}u(t)\big{)}\,,\quad
t\in[0,\tau]\,,$

see also (4). This formalism is expedient with regard to the design of novel
higher-order time integration methods for nonlinear partial differential
equations. Though, it is then equally of importance to concretise formal
considerations and to confirm that the resulting modified operator splitting
methods are indeed practicable and beneficial.

  4. (iv)

_Main objectives._ In this work, for the sake of concreteness, we focus on the
extension of a famous fourth-order modified potential operator splitting
method by Chin [14] to the Gross–Pitaevskii equation (2) and the parabolic
equation (3). So far, this scheme has been introduced and studied merely for
the linear case.

Outline. The present manuscript is organised as follows. In Section 2, we
review fundamental concepts for operator splitting methods. In Sections 3 and
4, we state the formal generalisation of a fourth-order modified potential
operator splitting method to the nonlinear case and substantiate it for the
Gross–Pitaevskii equation (2) and its parabolic analogue (3). A fundamental
invariance principle that includes a known result for standard splitting
methods as a special case is deduced in Section 5. In Section 6, we detail the
implementation of the novel modified operator splitting method based on a
Fourier spectral space discretisation and provide numerical comparisons with
standard splitting methods. Additional information on a publicly accessible
Matlab code is found in Appendix A. The observed order reduction of Yoshida’s
fourth-order complex splitting method (8d) is analysed in Appendix B.

##  2 Survey of standard and modified potential splitting methods

Linear case. As an illustrative example, we state the simplest representative
of standard splitting methods, the first-order Lie–Trotter splitting method,
for a system of linear differential equations

$\begin{cases}\tfrac{\text{d}}{\text{d}t}\,u(t)=A\,u(t)+B\,u(t)\,,\\\
u(0)=u_{0}\,,\quad t\in[0,T]\,,\end{cases}$ (6)

defined by non-commuting time-independent complex matrices
$A,B\in\mathbb{C}^{M\times M}$. For a positive integer number $N\in\mathbb{N}$
with associated time increment and equidistant grid points

$\tau=\tfrac{T}{N}\,,\quad t_{n}=n\,\tau\,,\quad n\in\\{0,1,\dots,N\\}\,,$

numerical approximations to the exact solution values are obtained by the
recurrence

$u_{n+1}=\text{e}^{\tau B}\,\text{e}^{\tau A}\,u_{n}\approx u(t_{n+1})\,,\quad
n\in\\{0,1,\dots,N-1\\}\,.$

Higher-order splitting methods for (6) involve the action of several matrix
exponentials on the current approximation and can be cast into the format

$u_{n+1}=\text{e}^{b_{s}\tau B}\,\text{e}^{a_{s}\tau
A}\,\cdots\,\text{e}^{b_{1}\tau B}\,\text{e}^{a_{1}\tau A}\,u_{n}\approx
u(t_{n+1})\,,\quad n\in\\{0,1,\dots,N-1\\}\,,$

with real or complex coefficients $(a_{j},b_{j})_{j=1}^{s}$, respectively.

Nonlinear case. Their generalisation to nonlinear evolution equations (1) is
based on the composition of the solutions to the subproblems defined by
$F_{1}$ and $F_{2}$. We henceforth employ the compact notation

$\begin{gathered}\tfrac{\text{d}}{\text{d}t}\,u_{1}(t)=\alpha
F_{1}\big{(}u_{1}(t)\big{)}\,,\quad\mathcal{E}_{\tau,\alpha
F_{1}}\big{(}u_{1}(t_{n})\big{)}=u_{1}(t_{n}+\tau)\,,\\\
\tfrac{\text{d}}{\text{d}t}\,u_{2}(t)=\beta
F_{2}\big{(}u_{2}(t)\big{)}\,,\quad\mathcal{E}_{\tau,\beta
F_{2}}\big{(}u_{2}(t_{n})\big{)}=u_{2}(t_{n}+\tau)\,,\\\
\alpha,\beta\in\mathbb{C}\,,\quad t\in[t_{n},t_{n}+\tau]\,,\end{gathered}$
(7a) so that a higher-order splitting method applied to (1) reads as
$\begin{gathered}u_{n+1}=\big{(}\mathcal{E}_{\tau,b_{s}F_{2}}\circ\mathcal{E}_{\tau,a_{s}F_{1}}\circ\dots\circ\mathcal{E}_{\tau,b_{1}F_{2}}\circ\mathcal{E}_{\tau,a_{1}F_{1}}\big{)}(u_{n})\approx
u(t_{n+1})\,,\\\ n\in\\{0,1,\dots,N-1\\}\,.\end{gathered}$ (7b)

Schemes. In view of numerical comparisons, we introduce the coefficients of
the first-order Lie–Trotter splitting method

$s=1\,,\quad a_{1}=1\,,\quad b_{1}=1\,,$ (8a) and the second-order Strang
splitting method $s=2\,,\quad a_{1}=0\,,\quad a_{2}=1\,,\quad
b_{1}=\tfrac{1}{2}=b_{2}\,.$ (8b) The well-known fourth-order splitting method
by Yoshida [36] involves four stages $\begin{gathered}s=4\,,\quad
a_{1}=0\,,\quad a_{2}=1-2\,b_{2}=a_{4}\,,\quad a_{3}=4\,b_{2}-1\,,\\\
b_{1}=\tfrac{1}{2}-b_{2}=b_{4}\,,\quad
b_{2}=\tfrac{1}{6}\,\big{(}1-\sqrt[3]{2}-\tfrac{1}{2}\sqrt[3]{4}\big{)}=b_{3}\,.\end{gathered}$
(8c) Reconsidering its construction based on a triple jump composition of the
Strang splitting method, a corresponding fourth-order splitting method with
complex coefficients is obtained $\begin{gathered}s=4\,,\quad a_{1}=0\,,\quad
a_{2}=1-2\,b_{2}=a_{4}\,,\quad a_{3}=4\,b_{2}-1\,,\\\
b_{1}=\tfrac{1}{2}-b_{2}=b_{4}\,,\quad
b_{2}=\tfrac{1}{6}\,\big{(}1+\tfrac{1}{2}\sqrt[3]{2}+\tfrac{1}{4}\sqrt[3]{4}\big{)}+\text{i}\,\tfrac{\sqrt{3}}{12}\,\big{(}\tfrac{1}{2}\sqrt[3]{4}-\sqrt[3]{2}\big{)}=b_{3}\,,\end{gathered}$
(8d) see also [8, 11, 18].

Stability issues. In connection with the time integration of dissipative
systems and parabolic equations as well as the imaginary time propagation of
Schrödinger equations, operator splitting methods are subject to additional
stability constraints. In order to explain this matter of fact, we recall that
the application of a splitting method with real coefficients
$(a_{j},b_{j})_{j=1}^{s}$ to the Gross–Pitaevskii equation (2) in imaginary
time and to the parabolic equation (3) involves the subproblems

$\partial_{t}U(x,t)=a_{j}\,\Delta
U(x,t)\,,\quad(x,t)\in\Omega\times[t_{n},t_{n}+\tau]\,,\quad
j\in\\{1,\dots,s\\}\,.$

Evidently, requiring well-posedness of these subproblems or stability of the
resulting splitting method, respectively, implies

$a_{j}\geq 0\,,\quad j\in\\{1,\dots,s\\}\,.$

This positivity condition, however, excludes higher-order schemes, since any
splitting method that exceeds a second-order barrier necessarily comprises
negative coefficients, see for example [7, 30, 31]. Specifically, this holds
true for the fourth-order splitting method by Yoshida

$a_{3}\approx-\,1.7<0\,,$

see (8c). A feasible remedy to this issue is the design of splitting methods
with complex coefficients $(a_{j},b_{j})_{j=1}^{s}$ such that

$\Re(a_{j})\geq 0\,,\quad j\in\\{1,\dots,s\\}\,,$

see, e.g., [13, 19]. The fourth-order scheme (8d) indeed fulfills these
constraints

$a_{1}=0\,,\quad\Re(a_{2})=\Re(a_{4})\approx 0.3>0\,,\quad\Re(a_{3})\approx
0.4>0\,.$

For further considerations in the context of the imaginary time propagation of
the linear Schrödinger equation by complex splitting methods, we refer to [4].

Modified potential splitting methods. Reviving former work by Ruth and Suzuki,
see for instance [28, 31, 32], a favourable alternative to standard operator
splitting methods was proposed by Chin. In a seminal work [14], he developed a
famous fourth-order scheme of splitting type that comprises positive
coefficients and hence overcomes the second-order barrier for standard
splitting methods with real coefficients. Expressed in his own words, the
basic idea is to incorporate _an additional higher order composite operator_
so that the implementation of _one algorithm requires only one evaluation of
the force and one evaluation of the force and its gradient_. For linear
evolution equations (6), the resulting scheme takes the form

$\begin{gathered}u_{n+1}=\text{e}^{\frac{1}{6}\tau
B}\,\text{e}^{\frac{1}{2}\tau A}\,\text{e}^{\frac{2}{3}\tau
B-\frac{1}{72}\tau^{3}[B,[B,A]]}\,\text{e}^{\frac{1}{2}\tau
A}\,\text{e}^{\frac{1}{6}\tau B}\,u_{n}\approx u(t_{n+1})\,,\\\
n\in\\{0,1,\dots,N-1\\}\,.\end{gathered}$ (9)

Here, the iterated commutator of complex matrices $A,B\in\mathbb{C}^{M\times
M}$ is given by

$[B,A]=BA-AB\,,\quad\big{[}B,[B,A]\big{]}=B^{2}A-2\,BAB+AB^{2}\,,$ (10a) see
also (5). More generally, for linear differential and multiplication operators
$\begin{gathered}(A\,v)(x)=c\,\Delta\,v(x)\,,\\\
(B\,v)(x)=\bar{c}\,V(x)\,v(x)\,,\\\ x\in\Omega\,,\quad
c\in\mathbb{C}\,,\end{gathered}$ (10b) retained from (2) and (3) for the
special choice $\vartheta=0$, a straightforward calculation yields a linear
multiplication operator that depends on the Euclidean norm of the gradient of
the potential $\begin{split}&\big{[}B,[B,A]\big{]}v(x)\\\
&\quad=\bar{c}\,|c|^{2}\,\bigg{(}\big{(}V(x)\big{)}^{2}\Delta\,v(x)-2\,V(x)\,\Delta\big{(}V(x)\,v(x)\big{)}\\\
&\quad\qquad+\Delta\Big{(}\big{(}V(x)\big{)}^{2}\,v(x)\Big{)}\bigg{)}\\\
&\quad=2\,\bar{c}\,|c|^{2}\,\big{(}\nabla V(x)\big{)}^{T}\,\nabla
V(x)\,v(x)\,,\quad x\in\Omega\,.\end{split}$ (10c)

This explains the common notion _force-gradient operator splitting method_ or
_modified potential operator splitting method_ for the scheme (9) and related
splitting methods in the context of classical or quantum many-body problems
and beyond. More recent contributions that exploit (9) for linear ordinary and
partial differential equations are, e.g., [2, 15, 26, 27].

Generalisations to specific classes of nonlinear evolution equations. We point
out that the operator $B$ defined in (10b) and the iterated commutator (10c)
are of the same nature. This structural similarity explains the efficiency of
Chin’s scheme (9) for partial differential equations comprising the Laplacian
and a potential acting as multiplication operator. Anticipating the detailed
expositions in Sections 3 and 4, we stress that the generalisations to the
time-dependent Gross–Pitaevskii equation (2) and its parabolic analogue (3)
maintain this feature, even though the iterated commutators are more involved.

##  3 Modified operator splitting method

Formal generalisation. Our educated guess to formally generalise the modified
potential operator splitting method (9) to the significantly more involved
case of a nonlinear evolution equation (1) is

$\begin{gathered}u_{n+1}=\big{(}\mathcal{E}_{\tau,\frac{1}{6}F_{2}}\circ\mathcal{E}_{\tau,\frac{1}{2}F_{1}}\circ\mathcal{E}_{\tau,\frac{2}{3}F_{2}-\frac{1}{72}\tau^{2}G_{2}}\circ\mathcal{E}_{\tau,\frac{1}{2}F_{1}}\circ\mathcal{E}_{\tau,\frac{1}{6}F_{2}}\big{)}\,u_{n}\approx
u(t_{n+1})\,,\\\ n\in\\{0,1,\dots,N-1\\}\,.\end{gathered}$ (11a) That is, we
replace the matrices $A$ and $B$ by the nonlinear operators $F_{1}$ and
$F_{2}$ as well as the iterated commutator of matrices by the following
operator
$\begin{gathered}G_{1}(v)=F_{2}^{\prime}(v)\,F_{1}(v)-F_{1}^{\prime}(v)\,F_{2}(v)\,,\\\
G_{2}(v)=F_{2}^{\prime}(v)\,G_{1}(v)-G_{1}^{\prime}(v)\,F_{2}(v)\,,\end{gathered}$
(11b)

see (4) and (5). We determine the Gâteaux derivatives generalising directional
derivatives through

$H^{\prime}(v)\,w=\lim_{\varepsilon\to
0}\tfrac{1}{\varepsilon}\,\big{(}H(v+\varepsilon\,w)-H(v)\big{)}\,,$ (12)

see also [20], and presume well-definedness of the arising operators on
suitably chosen domains. According to (7), the decisive operator is associated
with a nonlinear evolution equation that comprises the time increment as
parameter

$\begin{gathered}\tfrac{\text{d}}{\text{d}t}\,\widetilde{u}_{2}(t)=\beta_{1}F_{2}\big{(}\widetilde{u}_{2}(t)\big{)}+\beta_{2}\,\tau^{2}\,G_{2}\big{(}\widetilde{u}_{2}(t)\big{)}\,,\quad
t\in[t_{n},t_{n}+\tau]\,,\\\
\mathcal{E}_{\tau,\,\beta_{1}F_{2}+\,\beta_{2}\tau^{2}G_{2}}\big{(}\widetilde{u}_{2}(t_{n})\big{)}=\widetilde{u}_{2}(t_{n}+\tau)\,,\quad\beta_{1},\beta_{2}\in\mathbb{R}\,.\end{gathered}$
(13)

Specification and implementation. In the subsequent sections, we complete the
remaining tasks. We first specify the iterated commutator (11b) for the time-
dependent Gross–Pitaevskii equation (2) and contrast it to the result obtained
for the parabolic counterpart (3). Then, we deduce an invariance principle
that has a substantial impact on the efficiency of the modified operator
splitting method (11) when applied to the Gross–Pitaevskii equation (2).
Implementation issues as well as strategies to reduce the computational cost
for parabolic equations are finally discussed in Section 6.

##  4 Iterated commutators

Generally speaking, the appropriate framework for the extension of iterated
commutators for matrices or linear operators, see (10a), to nonlinear
operators is provided by the formal calculus of Lie-derivatives. In the
present work, we concretise and verify the heuristic characterisation (11) for
relevant applications, the time-dependent Gross–Pitaevskii equation (2) and
the related parabolic equation (3). In the context of Schrödinger equations,
the arising functions $v,w:\Omega\to\mathbb{C}$ take complex values, whereas
it suffices to consider real-valued functions $v,w:\Omega\to\mathbb{R}$ for
parabolic problems. In both cases, suitable regularity requirements apply. For
notational simplicity, we omit the dependence on the spatial variable.

Derivatives. The Gâteaux derivatives of the linear differential and nonlinear
multiplication operators defined in (2) and (3) are given by

$\begin{gathered}F_{1}(v)=c\,\Delta v\,,\\\ F_{1}^{\prime}(v)\,w=c\,\Delta
w\,,\\\
F_{2}(v)=\bar{c}\,\big{(}V+\vartheta\,|v|^{2}\big{)}\,v=\bar{c}\,\big{(}Vv+\vartheta\,v^{2}\,\bar{v}\big{)}\,,\\\
F_{2}^{\prime}(v)\,w=\bar{c}\,\big{(}Vw+2\,\vartheta\,\bar{v}\,v\,w+\vartheta\,v^{2}\,\overline{w}\big{)}\,,\end{gathered}$

see also (12).

First commutators. On the one hand, by performing differentiation twice, we
obtain

$\begin{split}&F_{1}^{\prime}(v)\,F_{2}(v)\\\
&\quad=|c|^{2}\,\Delta\big{(}Vv+\vartheta\,v^{2}\,\bar{v}\big{)}\\\
&\quad=|c|^{2}\,\big{(}\Delta Vv+2\,(\nabla V)^{T}\,\nabla v+V\Delta
v\big{)}\\\ &\quad\quad\,+|c|^{2}\,\vartheta\,\big{(}2\,\Delta
v\,\bar{v}\,v+\Delta\bar{v}\,v^{2}+2\,(\nabla v)^{T}\,\nabla
v\,\bar{v}+4\,(\nabla v)^{T}\,\nabla\bar{v}\,v\big{)}\,.\end{split}$

On the other hand, a simple replacement yields

$\begin{split}&F_{2}^{\prime}(v)\,F_{1}(v)=|c|^{2}\,\big{(}V\Delta
v+2\,\vartheta\,\bar{v}\,v\,\Delta
v\big{)}+\bar{c}^{2}\,\vartheta\,v^{2}\Delta\bar{v}\,.\end{split}$

As a consequence, the difference is given by

$\begin{split}&G_{1}(v)=F_{2}^{\prime}(v)\,F_{1}(v)-F_{1}^{\prime}(v)\,F_{2}(v)\\\
&\quad=-\,|c|^{2}\,\big{(}\Delta Vv+2\,(\nabla V)^{T}\,\nabla v\big{)}\\\
&\quad\quad\,+\big{(}\bar{c}^{2}-|c|^{2}\big{)}\,\vartheta\,\Delta\bar{v}\,v^{2}-2\,|c|^{2}\,\vartheta\,\big{(}(\nabla
v)^{T}\,\nabla v\,\bar{v}+2\,(\nabla
v)^{T}\,\nabla\bar{v}\,v\big{)}\,.\end{split}$

For the parabolic case (3) with $c=1$ this implies

$G_{1}(v)=-\,\Delta Vv-2\,(\nabla V)^{T}\,\nabla v-6\,\vartheta\,(\nabla
v)^{T}\,\nabla v\,v\,,$

and the analogous result for the Schrödinger case (2) with $c=\text{i}$ is

$G_{1}(v)=-\,\Delta Vv-2\,(\nabla V)^{T}\,\nabla
v-2\,\vartheta\,\big{(}\Delta\bar{v}\,v^{2}+(\nabla v)^{T}\,\nabla
v\,\bar{v}+2\,(\nabla v)^{T}\,\nabla\bar{v}\,v\big{)}\,.$

Iterated commutators. The iterated commutator associated with the parabolic
equation (3) results from straighforward but lengthy calculations

$\begin{gathered}G_{2}(v)=2\,\big{(}(\nabla V)^{T}(\nabla
V)+\vartheta\,\widetilde{G}_{2}(v)\big{)}\,v\,,\\\
\widetilde{G}_{2}(v)=-\,\Delta V\,v^{2}+6\,(\nabla V)^{T}(\nabla
v)\,v+6\,\big{(}V+2\,\vartheta\,v^{2}\big{)}\,(\nabla v)^{T}(\nabla
v)\,,\end{gathered}$ (14)

and for the Gross–Pitaevskii equation (2), we instead arrive at

$\begin{gathered}G_{2}(v)=-\,2\,\text{i}\,\Big{(}(\nabla V)^{T}(\nabla
V)-2\,\vartheta\,\big{(}\widetilde{G}_{21}(v)+\vartheta\,\widetilde{G}_{22}(v)\big{)}\Big{)}\,v\,,\\\
\widetilde{G}_{21}(v)=|v|^{2}\,\Delta V\,,\\\
\widetilde{G}_{22}(v)=|v|^{2}\,\big{(}2\,\Re(\bar{v}\,\Delta
v)+3\,(\nabla\bar{v})^{T}(\nabla v)\big{)}+\Re\big{(}\bar{v}^{2}\,(\nabla
v)^{T}(\nabla v)\big{)}\,.\end{gathered}$ (15)

For the special case of linear multiplication operators, i.e. $\vartheta=0$,
we indeed recover (10c). It is also noteworthy that the operator in (15)
involves $\Delta v$ and hence implies stronger regularity requirements on $v$
compared to the operator in (14), which only comprises the gradient $\nabla
v$.

##  5 Invariance principle

Invariance principle. For the purpose of illustration, we introduce the
nonlinear Schrödinger equation

$\begin{cases}\text{i}\,\partial_{t}\Psi(x,t)=-\,\Delta\Psi(x,t)+|\Psi(x,t)|^{2}\,\Psi(x,t)\,,\\\
\Psi(x,0)=\Psi_{0}(x)\,,\quad(x,t)\in\Omega\times[0,T]\,,\end{cases}$

retained from (2) for vanishing potential $V=0$ and normalised constant
$\vartheta=1$. Any standard operator splitting method (7) relies on the time
integration of the subproblem involving the Laplacian

$\begin{cases}\text{i}\,\partial_{t}\Psi_{1}(x,t)=-\,\alpha\,\Delta\Psi_{1}(x,t)\,,\quad\alpha\in\mathbb{R}\,,\\\
\Psi_{1}(x,t_{n})\text{
given}\,,\quad(x,t)\in\Omega\times[t_{n},t_{n}+\tau]\,,\end{cases}$

which we perform by means of a Fourier spectral space discretisation, see
Section 6. Moreover, we make use of the fact that the subproblem comprising
the cubic nonlinearity

$\begin{cases}\text{i}\,\partial_{t}\Psi_{2}(x,t)=\beta_{1}\,|\Psi_{2}(x,t)|^{2}\,\Psi_{2}(x,t)\,,\quad\beta_{1}\in\mathbb{R}\,,\\\
\Psi_{2}(x,t_{n})\text{
given}\,,\quad(x,t)\in\Omega\times[t_{n},t_{n}+\tau]\,,\end{cases}$ (16a)
satisfies the fundamental invariance principle
$|\Psi_{2}(x,t)|^{2}=|\Psi_{2}(x,t_{n})|^{2}\,,\quad(x,t)\in\Omega\times[t_{n},t_{n}+\tau]\,.$
This identity follows by means of differentiation with respect to time and
substitution of the time derivative by the nonlinearity
$\begin{gathered}\partial_{t}\,|\Psi_{2}(x,t)|^{2}=2\,\Re\big{(}\overline{\Psi_{2}(x,t)}\,\partial_{t}\Psi_{2}(x,t)\big{)}=-\,2\,\beta_{1}\,\Re\big{(}\text{i}\,|\Psi_{2}(x,t)|^{4}\big{)}=0\,,\\\
(x,t)\in\Omega\times[t_{n},t_{n}+\tau]\,.\end{gathered}$ It implies that the
exact solution to the nonlinear subproblem satisfies the reduced problem
$\begin{cases}\text{i}\,\partial_{t}\Psi_{2}(x,t)=\beta_{1}\,|\Psi_{2}(x,t_{n})|^{2}\,\Psi_{2}(x,t)\,,\quad\beta_{1}\in\mathbb{R}\,,\\\
\Psi_{2}(x,t_{n})\text{
given}\,,\quad(x,t)\in\Omega\times[t_{n},t_{n}+\tau]\,,\end{cases}$ (16b) and
hence can be determined by pointwise multiplication
$\Psi_{2}(x,t)=\text{e}^{-\,\text{i}\,\beta_{1}\,|\Psi_{2}(x,t_{n})|^{2}}\,\Psi_{2}(x,t_{n})\,,\quad(x,t)\in\Omega\times[t_{n},t_{n}+\tau]\,.$
(16c)

For the numerical realisation, equidistant grid points are chosen accordingly
to the Fourier spectral space discretisation, see also Section 6.

Generalisation. In this section, we establish a substantial extension of the
above stated result (16) that has important implications concerning the
efficient implementation of the modified operator splitting method (11) for
the time-dependent Gross–Pitaevskii equation (2), see also (13) and (15).
Furthermore, the obtained invariance principle for the evolution operator

$\mathcal{E}_{\tau,\,\beta_{1}F_{2}+\,\beta_{2}\tau^{2}G_{2}}\,,\quad\tau,\beta_{1},\beta_{2}\in\mathbb{R}\,,$

is connected with the significance of the modified operator splitting method
as a geometric integrator, see [18, 21] and references given therein.

Notation. With regard to a compact formulation as abstract evolution equation,
we again omit the dependence of the potential and a regular complex-valued
function $v:\Omega\to\mathbb{C}$ on the spatial variable. Besides, for
accomplishing relations of the form

$F_{2}(v)=\bar{c}\,f_{1}(v)\,v\,,\quad G_{2}(v)=\bar{c}\,f_{2}(v)\,v\,,$

it is convenient to introduce the abbreviations

$\begin{gathered}f(v)=\beta_{1}\,f_{1}(v)+\beta_{2}\,\tau^{2}f_{2}(v)\,,\\\
f_{1}(v)=V+\vartheta\,g_{1}(v)\,,\quad f_{2}(v)=2\,(\nabla V)^{T}(\nabla
V)-4\,\vartheta\,g_{6}(v)\,,\\\ g_{1}(v)=|v|^{2}\,,\quad
g_{2}(v)=\Re(\bar{v}\,\Delta v)\,,\\\ g_{3}(v)=(\nabla\bar{v})^{T}(\nabla
v)\,,\quad g_{4}(v)=\Re\big{(}\bar{v}^{2}\,(\nabla v)^{T}(\nabla
v)\big{)}\,,\\\
g_{5}(v)=\vartheta\,\big{(}2\,g_{2}(v)+3\,g_{3}(v)\big{)}\,,\quad
g_{6}(v)=g_{1}(v)\,\big{(}\Delta
V+g_{5}(v)\big{)}+\vartheta\,g_{4}(v)\,.\end{gathered}$ (17)

Theorem (Invariance principle). The solution to the subproblem

$\begin{cases}\tfrac{\text{d}}{\text{d}t}\,\psi(t)=-\,\text{i}\,f\big{(}\psi(t)\big{)}\,\psi(t)\,,\\\
\psi(0)=\psi_{0}\,,\quad t\in[0,\tau]\,,\end{cases}$

with defining function introduced in (17) satisfies the invariance principle

$f\big{(}\psi(t)\big{)}=f(\psi_{0})\,,\quad t\in[0,\tau]\,.$

Proof. In order to demonstrate that the invariance principle holds, we
determine the Gâteaux derivatives of the defining functions

$\begin{gathered}f^{\prime}(v)\,w=\beta_{1}\,f_{1}^{\prime}(v)\,w+\beta_{2}\,\tau^{2}f_{2}^{\prime}(v)\,w\,,\\\
f_{1}^{\prime}(v)\,w=\vartheta\,g_{1}^{\prime}(v)\,w\,,\quad
f_{2}^{\prime}(v)\,w=-\,4\,\vartheta\,g_{6}^{\prime}(v)\,w\,,\\\
g_{1}^{\prime}(v)\,w=2\,\Re(\bar{v}\,w)\,,\quad
g_{2}^{\prime}(v)\,w=\Re(\Delta\bar{v}\,w)+\Re(\bar{v}\,\Delta w)\,,\\\
g_{3}^{\prime}(v)\,w=2\,\Re\big{(}(\nabla\bar{v})^{T}(\nabla w)\big{)}\,,\\\
g_{4}^{\prime}(v)\,w=2\,\Re\big{(}\bar{v}\,\overline{w}\,(\nabla v)^{T}(\nabla
v)\big{)}+2\,\Re\big{(}\bar{v}^{2}\,(\nabla v)^{T}(\nabla w)\big{)}\,,\\\
g_{5}^{\prime}(v)\,w=2\,\vartheta\,g_{2}^{\prime}(v)\,w+3\,\vartheta\,g_{3}^{\prime}(v)\,w\,,\\\
g_{6}^{\prime}(v)\,w=\big{(}\Delta
V+g_{5}(v)\big{)}\,g_{1}^{\prime}(v)\,w+g_{1}(v)\,g_{5}^{\prime}(v)\,w+\vartheta\,g_{4}^{\prime}(v)\,w\,,\\\
\end{gathered}$

where suitable regularity requirements apply to $v,w:\Omega\to\mathbb{C}$.
Observing that the potential $V$ and the basic components
$g_{1},g_{2},g_{3},g_{4}$ define real-valued functions, we have

$\begin{gathered}\big{(}g_{j}(v)\big{)}(x)\in\mathbb{R}\,,\quad
j\in\\{1,\dots,6\\}\,,\\\ \big{(}f_{k}(v)\big{)}(x)\in\mathbb{R}\,,\quad
k\in\\{1,2\\}\,,\quad\big{(}f(v)\big{)}(x)\in\mathbb{R}\,,\quad
x\in\Omega\,.\end{gathered}$

Evidently, this implies that the following composition vanishes

$g_{1}^{\prime}(v)\,\big{(}\text{i}\,f(v)\,v\big{)}=2\,\Re\big{(}\text{i}\,|v|^{2}\,f(v)\big{)}=0\,.$

Certain contributions originating from the iterated commutator, however,
require a closer examination, namely

$\begin{split}&g_{2}^{\prime}(v)\,\big{(}\text{i}\,f(v)\,v\big{)}\\\
&\quad=\Re\Big{(}\text{i}\,f(v)\,\Re\big{(}\bar{v}\,\Delta
v\big{)}\Big{)}+\Re\big{(}\text{i}\,|v|^{2}\Delta
f(v)\big{)}+2\,\Re\Big{(}\text{i}\,\bar{v}\,\big{(}\nabla
f(v)\big{)}^{T}\,\nabla v\Big{)}\\\
&\quad=2\,\Re\Big{(}\text{i}\,\bar{v}\,\big{(}\nabla f(v)\big{)}^{T}\nabla
v\Big{)}\,,\\\ &g_{3}^{\prime}(v)\,\big{(}\text{i}\,f(v)\,v\big{)}\\\
&\quad=2\,\Re\Big{(}\text{i}\,v\,\big{(}\nabla
f(v)\big{)}^{T}\,\nabla\bar{v}\Big{)}+2\,\Re\Big{(}\text{i}\,f(v)\,(\nabla\bar{v})^{T}\,\nabla
v\Big{)}\\\ &\quad=2\,\Re\Big{(}\text{i}\,v\,\big{(}\nabla
f(v)\big{)}^{T}\,\nabla\bar{v}\Big{)}\,,\\\
&g_{4}^{\prime}(v)\,\big{(}\text{i}\,f(v)\,v\big{)}\\\
&\quad=-\,2\,\Re\big{(}\text{i}\,f(v)\,\bar{v}^{2}\,(\nabla v)^{T}\,\nabla
v\big{)}+2\,\Re\Big{(}\text{i}\,|v|^{2}\,\bar{v}\,\big{(}\nabla
f(v)\big{)}^{T}\nabla v\Big{)}\\\
&\qquad+2\,\Re\big{(}\text{i}\,f(v)\,\bar{v}^{2}\,(\nabla v)^{T}\,\nabla
v\big{)}\\\ &\quad=2\,\Re\Big{(}\text{i}\,|v|^{2}\,\bar{v}\,\big{(}\nabla
f(v)\big{)}^{T}\nabla v\Big{)}\,.\end{split}$

On the basis of these identities, we conclude

$\begin{split}&\Big{(}g_{1}(v)\,\big{(}2\,g_{2}^{\prime}(v)+3\,g_{3}^{\prime}(v)\big{)}+g_{4}^{\prime}(v)\big{)}\Big{)}\,\big{(}\text{i}\,f(v)\,v\big{)}\\\
&\quad=12\,\Re\Big{(}\text{i}\,|v|^{2}\,\big{(}\nabla
f(v)\big{)}^{T}\Re\big{(}\bar{v}\,\nabla v\big{)}\Big{)}=0\,.\end{split}$

This proves that any composition of the special form

$\begin{split}&f^{\prime}(v)\,\big{(}\text{i}\,f(v)\,v\big{)}\\\
&\quad=\vartheta\,\Big{(}\beta_{1}-4\,\beta_{2}\,\tau^{2}\,\big{(}\Delta
V+g_{5}(v)\big{)}\Big{)}\,g_{1}^{\prime}(v)\,\big{(}\text{i}\,f(v)\,v\big{)}\\\
&\qquad\,-4\,\vartheta^{2}\beta_{2}\,\tau^{2}\Big{(}g_{1}(v)\,\big{(}2\,g_{2}^{\prime}(v)+3\,g_{3}^{\prime}(v)\big{)}+g_{4}^{\prime}(v)\Big{)}\,\big{(}\text{i}\,f(v)\,v\big{)}\\\
&\quad=0\end{split}$

vanishes. As a consequence, the time-derivative of the decisive function is
equal to zero

$\tfrac{\text{d}}{\text{d}t}\,f\big{(}\psi(t)\big{)}=f^{\prime}\big{(}\psi(t)\big{)}\,\tfrac{\text{d}}{\text{d}t}\,\psi(t)=-\,f^{\prime}\big{(}\psi(t)\big{)}\,\Big{(}\text{i}\,f\big{(}\psi(t)\big{)}\,\psi(t)\Big{)}=0\,,\quad
t\in[0,\tau]\,,$

and hence, the desired identity follows

$f\big{(}\psi(t)\big{)}=f(\psi_{0})\,,\quad t\in[0,\tau]\,.$

Summary. In view of Section 6, we summarise our considerations for the
modified operator splitting method (11) applied to the time-dependent
Gross–Pitaevskii equation (2). On the one hand, the realisation of

$\begin{gathered}\psi_{n+1}=\big{(}\mathcal{E}_{\tau,\frac{1}{6}F_{2}}\circ\mathcal{E}_{\tau,\frac{1}{2}F_{1}}\circ\mathcal{E}_{\tau,\frac{2}{3}F_{2}-\frac{1}{72}\tau^{2}G_{2}}\circ\mathcal{E}_{\tau,\frac{1}{2}F_{1}}\circ\mathcal{E}_{\tau,\frac{1}{6}F_{2}}\big{)}\,\psi_{n}\,,\\\
n\in\\{0,1,\dots,N-1\\}\,,\end{gathered}$

relies on the numerical integration of the linear Schrödinger equation

$\begin{gathered}\tfrac{\text{d}}{\text{d}t}\,\psi(t)=\text{i}\,\alpha\,\Delta\psi(t)\,,\quad
t\in[t_{n},t_{n}+\tau]\,,\quad\mathcal{E}_{\tau,\alpha
F_{1}}\big{(}\psi(t_{n})\big{)}=\psi(t_{n}+\tau)\,,\end{gathered}$

for $\alpha\in\mathbb{R}$. On the other hand, it reduces to the pointwise
evaluation of the solution representation

$\mathcal{E}_{\tau,\,\beta_{1}F_{2}+\,\beta_{2}\tau^{2}G_{2}}(\psi_{0})=\text{e}^{-\,\text{i}\,\tau(\beta_{1}f_{1}(\psi_{0})+\beta_{2}\tau^{2}f_{2}(\psi_{0}))}\,\psi_{0}\,,\quad\tau\in\mathbb{R}\,,$
(18)

for appropriate choices of $\beta_{1},\beta_{2}\in\mathbb{R}$, see (17) for
the definitions of $f_{1}$ and $f_{2}$. We recall that the known identity (16)
for the particular case $\beta_{2}=0$ is established by a simplified argument.

##  6 Numerical results

In the following, we illustrate the stability and global error behaviour of
the novel modified operator splitting method (11) for the time-dependent
Gross–Pitaevskii equation (2) and its parabolic analogue (3). The numerical
tests, performed in one, two, and three space dimensions, in particular
confirm the theoretical considerations of Sections 3 to 5. For the purpose of
comparison, we in addition include the corresponding results for widely-used
standard splitting methods. Further information on a publicly accessible
Matlab code to reproduce Figures 1–4 is found in Appendix A.

Implementation. The practical realisation of standard operator splitting
methods such as (8) and of the modified operator splitting method (11),
respectively, requires the time integration of the subproblems involving the
Laplacian and the nonlinear multiplication operator. In our implementation, we
make use of fast Fourier techniques, which are based on the following
considerations.

  1. (i)

_Space grid._ With regard to the employed Fourier spectral space
discretisation, we replace the underlying unbounded domain by a Cartesian
product of sufficiently large intervals

$a=10\,,\quad x\in[-\,a,a]^{d}\subset\Omega\,,$

and choose the total number of equidistant spatial grid points according to
the dimension

$d=1:\;M=512\,,\quad d=2:\;M=128^{2}\,,\quad d=3:\;M=64^{3}\,.$

  2. (ii)

_Derivatives._ The iterated commutator arising in the modified operator
splitting method (11) for (2) and (3) involves the gradient $\nabla V$ and the
Laplacian $\Delta V$ of the space-dependent potential, which we may assume to
be known analytically. Otherwise, we employ the approach described
subsequently. The numerical computation of the spatial derivatives $\nabla v$
and $\Delta v$, where $v$ represents the current value of the time-discrete
solution, is traced back to a fast Fourier transform, pointwise
multiplications, and an inverse fast Fourier transform. Denoting by
$(\mathcal{F}_{m})_{m\in\mathbb{Z}^{d}}$ the Fourier functions with
periodicity domain $[-\,a,a]^{d}$, by $(\mu_{m})_{m\in\mathbb{Z}^{d}}$ the
purely imaginary eigenvalues associated with the first spatial derivatives,
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}}\;\text{e}^{\,\text{i}\,\pi\,m_{1}\,(\frac{x_{1}}{a}+1)}\cdots\text{e}^{\,\text{i}\,\pi\,m_{d}\,(\frac{x_{d}}{a}+1)}\,,\\\
\nabla\mathcal{F}_{m}=\mu_{m}\,\mathcal{F}_{m}\,,\quad\mu_{m}=\tfrac{\text{i}\,\pi\,m}{a}\in\mathbb{C}^{d\times
1}\,,\\\
\Delta\,\mathcal{F}_{m}=\lambda_{m}\,\mathcal{F}_{m}\,,\quad\lambda_{m}=-\,\tfrac{\pi^{2}|m|^{2}}{a^{2}}\in\mathbb{R}\,,\\\
m=(m_{1},\dots,m_{d})\in\mathbb{Z}^{d}\,,\quad
x=(x_{1},\dots,x_{d})\in[-\,a,a]^{d}\,,\end{gathered}$

the following formal representations hold

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

Their realisation relies on a suitable truncation of the infinite index set
$\mathcal{M}\subset\mathbb{Z}^{d}$ so that $|\mathcal{M}|=M$ and quadrature
approximations by the trapezoidal rule.

  3. (iii)

_Linear subproblem._ Formally, the solution to the linear subproblem is given
by a Fourier series

$\begin{gathered}\tfrac{\text{d}}{\text{d}t}\,u_{1}(t)=\alpha
F_{1}\big{(}u(t)\big{)}=c\,\alpha\,\Delta\,u_{1}(t)\,,\quad
t\in[t_{n},t_{n}+\tau]\,,\quad\alpha\in\mathbb{R}\,,\\\
\mathcal{E}_{\tau,\alpha F_{1}}\big{(}u(t_{n})\big{)}=\mathcal{E}_{\tau,\alpha
F_{1}}\bigg{(}\sum_{m\in\mathbb{Z}^{d}}u_{1,m}(t_{n})\,\mathcal{F}_{m}\bigg{)}=\sum_{m\in\mathbb{Z}^{d}}\text{e}^{\,c\,\alpha\,\tau\lambda_{m}}\,u_{1,m}(t_{n})\,\mathcal{F}_{m}\,,\end{gathered}$

see also (2), (3), and (7). Again, the application of fast Fourier techniques
permits the efficient computation of approximations to the spectral
coefficients and the evaluation of finite sums on the equidistant spatial grid
points covering the underlying domain

$\begin{gathered}\widetilde{u}_{1,m,n}\approx u_{1,m}(t_{n})\,,\quad
m\in\mathcal{M}\,,\\\
\sum_{m\in\mathcal{M}}\text{e}^{\,c\,\alpha\,\tau\lambda_{m}}\,\widetilde{u}_{1,m,n}\,\mathcal{F}_{m}(x)\,,\quad
x\in[-\,a,a]^{d}\,.\end{gathered}$

  4. (iv)

_Nonlinear subproblem (Schrödinger equation)._ In the context of the
Gross–Pitaevskii equation (2), we make use of the fact that the solution to
the nonlinear subproblem (13) satisfies the invariance principle deduced in
Section 5. Consequently, it simply remains to evaluate the representation (18)
on the equidistant grid.

  5. (v)

_Nonlinear subproblem (Parabolic equation)._ In the case of the parabolic
equation (3), we additionally apply an explicit Runge–Kutta method of order
four. Due to the stiffness of the problem, the time stepsize has to be
adjusted to the spatial grid width to ensure stability. Alternative approaches
with improved stability properties and reduced computational costs are
detailed below.

Numerical results. In our numerical tests, we perform the time integration of
the Gross–Pitaevskii equation (2) and the parabolic equation (3), comparing
the nonlinear case with $\vartheta=1$ to the simplified linear case with
$\vartheta=0$. We prescribe the Gaussian-shaped initial state

$u_{0}(x)=\text{e}^{-\frac{1}{2}\,(x_{1}^{2}+\dots+x_{d}^{2})}\,,\quad
x=(x_{1},\dots,x_{d})\in\mathbb{R}^{d}\,,$

as well as the two polynomial potentials

$\begin{gathered}V(x)=C_{0}\,C_{q}\sum_{j=1}^{d}x_{j}^{q}\,,\\\
C_{2}=1\,,\quad C_{4}=\tfrac{1}{24}\,,\quad q\in\\{2,4\\}\,,\quad
x=(x_{1},\dots,x_{d})\in\mathbb{R}^{d}\,,\end{gathered}$

for which the required first- and second-order derivatives are known
analytically. For the special case of a quadratic potential with prefactor
$C_{0}\in\\{1,-\,1\\}$ chosen accordingly to the type of the equation and
$\vartheta=0$, the knowledge of the exact solution

$u(x,t)=\text{e}^{d\,\bar{c}\,t}\,u_{0}(x)\,,\quad(x,t)\in\mathbb{R}^{d}\times[0,T]\,,$

permits the validation of the Fourier spectral space discretisation and the
time-splitting approach. In the general case, we instead determine a numerical
reference solution based on a refined time stepsize. The global errors of the
modified operator splitting method (11) at final time $T=1$, measured in a
discrete $L_{2}$-norm, are compared to those obtained by the standard
Lie–Trotter, Strang, and Yoshida splitting methods of non-stiff orders one,
two, and four, see also (8). The obtained results, displayed in Figures 1–4,
confirm the favourable behaviour of the modified operator splitting method.

Energy conservation. In order to complement our numerical comparisons
regarding the stability and global error behaviour of standard and modified
operator splitting methods, we perform the time integration of the one-
dimensional Gross–Pitaevskii equation (2) on the interval $[0,T]$. We
determine approximations to the values of the energy at equidistant time grid
points

$\begin{split}E(t_{n})&=\int_{\Omega}\big{(}-\Delta\Psi(x,t_{n})+V(x)\,\Psi(x,t_{n})+\vartheta\,|\Psi(x,t_{n})|^{2}\,\Psi(x,t_{n})\big{)}\\\
&\qquad\quad\times\overline{\Psi(x,t_{n})}\;\text{d}x\,,\quad
n\in\\{0,1,\dots,N\\}\,,\end{split}$ (19a) and their deviations with respect
to the minimal value
$E(t_{n})-\min\big{\\{}E(t_{\ell}):\ell\in\\{0,1,\dots,N\\}\big{\\}}\,,\quad
n\in\\{0,1,\dots,N\\}\,.$ (19b)

The obtained results confirm the favourable geometric properties of the
modified operator splitting method (11), see Figure 6.

Computational cost. In general, an expedient measure for the computational
cost of the modified operator splitting method (11) for the Gross–Pitaevskii
equation (2) and the parabolic equation (3), respectively, is the number of
fast Fourier transforms and their inverses. Evidently, the numerical solution
of two linear subproblems per time step amounts to two fast Fourier transforms
and two inverse fast Fourier transforms. Besides, the time integration of the
nonlinear subproblem requires the computation of space derivatives via Fourier
transforms.

  1. (i)

For nonlinear Schrödinger equations such as (2), the validity of the
invariance property permits to significantly reduce the cost related to the
evaluation of $\mathcal{E}_{\tau,\frac{2}{3}F_{2}-\frac{1}{72}\tau^{2}G_{2}}$.
Due to the fact that the spectral coefficients of the current numerical
solution $v$ are available, the computation of gradient $\nabla v$ and
Laplacian $\Delta v$ results in $d+1$ inverse fast Fourier transforms.

  2. (ii)

For the parabolic equation (3), a favourable approach is based on the
following considerations. The presence of the additional factor $\tau^{2}$ in
connection with the double commutator permits to use an approximation by means
of the the second-oder Strang splitting method, that is

$\mathcal{E}_{\tau,\frac{2}{3}F_{2}-\frac{1}{72}\tau^{2}G_{2}}\approx\mathcal{E}_{\frac{1}{2}\tau,\frac{2}{3}F_{2}}\circ\mathcal{E}_{\tau,-\frac{1}{72}\tau^{2}G_{2}}\circ\mathcal{E}_{\frac{1}{2}\tau,\frac{2}{3}F_{2}}\,.$

For linear ordinary differential equations of the form

$\tfrac{\text{d}}{\text{d}t}\,u(t)=B\,u(t)+\tau^{2}\,C\,u(t)\,,\quad
t\in[0,\tau]\,,$

defined by complex matrices $B,C\in\mathbb{C}^{M\times M}$, elementary
calculations based on series expansions of matrix exponentials confirm that
this approach leads to a fourth-order approximation

$\text{e}^{\tau(B+\tau^{2}C)}=\text{e}^{\frac{1}{2}\tau
B}\,\text{e}^{\tau^{3}C}\,\text{e}^{\frac{1}{2}\tau
B}+\mathcal{O}\big{(}\tau^{5}\big{)}\,.$

The rigorous generalisation of the argument to (3) is part of a future
convergence analysis. On the one hand, an explicit representation of the
evolution operator associated with the nonlinear ordinary differential
equation

$\tfrac{\text{d}}{\text{d}t}\,u(x,t)=\big{(}V(x)+\vartheta\,|u(x,t)|^{2}\big{)}\,u(x,t)\,,\quad
x\in\Omega\,,\quad t\in[0,\tau]\,,$

is known. More precisely, fixing $x\in\Omega$, the distinction of cases yields

$u(x,t)=\begin{cases}\frac{u(x,0)}{\sqrt{\text{e}^{-2\,t\,V(x)}+\vartheta\,(\text{e}^{-2\,t\,V(x)}-1)/V(x)\,(u(x,0))^{2}}}\,,&V(x)\neq
0\,,\\\ \frac{u(x,0)}{\sqrt{1-2\,t\,\vartheta\,(u(x,0))^{2}}}\,,&V(x)=0\,,\\\
\end{cases}\quad t\in[0,\tau]\,.$

On the other hand, for the time integration of the nonlinear subproblem
involving the double commutator, it suffices to apply the first-order explicit
Euler method. Numerical tests confirm the enhanced stability, accuracy, and
efficiency of the resulting approach, in particular, in case of a fourth-order
polynomial potential, see Figure 5.

##  7 Conclusions

In the present work, we have introduced a general framework for the extension
of Chin’s fourth-order modified potential operator splitting method (9) to
nonlinear evolution equations. To the best of our knowledge, this matter is
novel and of major interest from theoretical and practical perspectives.

Moreover, we have specified the resulting fourth-order modified operator
splitting method (11) for the time-dependent Gross–Pitaevskii equation (2) and
its parabolic counterpart (3). It seems likely that our approach and the drawn
conclusions also extend to Schrödinger equations involving nonlinear and
possibly nonlocal terms of the form $g(|\Psi(x,t)|)\,\Psi(x,t)$ .

Due to the fact that our numerical tests have confirmed the favourable
performance of the proposed fourth-order scheme in comparison with standard
real and complex splitting methods of order four, it is natural to extend our
considerations in various respects.

Proceeding our recent work [10] for linear evolution equations of parabolic
and Schrödinger type, we find it promising to design high-order modified
operator splitting methods for nonlinear evolution equations that are optimal
with regard to a preselected criterium such as efficiency. In contrast to the
linear case, we have to take into account the additional costs for the
evaluation of second iterated commutators and that higher-order iterated
commutators will not vanish, in general, see (10c) as well as (14) and (15).
Nonetheless, based on the successful strategies for the efficient
implementation of the fourth-order scheme (11), it suggests itself to address
an in-depth analysis of similarly structured modified operator splitting
methods.

A desirable feature of the proposed fourth-order modified operator splitting
method is the positivity of the coefficients. Concerning the design of high-
order schemes, in light of contributions on the linear case, see for instance
[3, 15], this will necessitate further deliberations on the appropriate
format. Evidently, positivity is intrinsically related to the issue of well-
posedness and stability for parabolic problems. Besides, it affects aliasing
effects on truncated space domains and the incorporation of artificial
boundary conditions, which is of particular interest in the context of
nonlinear Schrödinger equations, see for example [1] and references given
therein.

Furthermore, we intend to carry out a rigorous convergence analysis of
modified operator splitting methods applied to the time-dependent
Gross–Pitaevskii equation and its parabolic analogue. We point out that the
study for the linear case [22] together with the formal calculus of Lie-
derivatives provides a guiding principle. However, the accomplishment for
specific nonlinear evolution equations involving unbounded operators requires
careful calculations and investigations, see also [33].

## Acknowledgements

The authors dedicate this work to Arieh Iserles due to his seminal
contributions in the area of numerical analysis and geometric numerical
integration. The authors are grateful to the two anonymous reviewers for their
careful reading of the manuscript and valuable comments. 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] X. Antoine, A. Arnold, Ch. Besse, M. Ehrhardt, A. Schädle. _A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations_. Commun. Comput. Phys. 4 (2008) 729–796.
  * [2] J. Auer, E. Krotscheck, S. Chin. _A fourth-order real-space algorithm for solving local Schrödinger equations_. J. Chem. Phys. 115 (2001) 6841–6846.
  * [3] W. Auzinger, H. Hofstätter, O. Koch. _Non-existence of generalized splitting methods with positive coefficients of order higher than four_. Appl. Math. Letters 97 (2019) 48–52.
  * [4] 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.
  * [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] 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.
  * [8] 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.
  * [9] S. Blanes, F. Casas, A. Murua. _An efficient algorithm based on splitting for the time integration of the Schrödinger equation_. J. Comput. Phys. 303 (2015) 396–412.
  * [10] 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_. To appear in Commun. Comput. Phys. (2023).
  * [11] S. Blanes, F. Casas, C. González, M. Thalhammer. _Symmetric-conjugate splitting methods for evolution equations of parabolic type_. Submitted (2023).
  * [12] M. Caliari, S. Zuccher. _A fast time splitting finite difference approach to Gross–Pitaevskii equations_. Commun. Comput. Phys. 29 (2021) 1336–1364.
  * [13] F. Castella, P. Chartier, S. Decombes, G. Vilmart. _Splitting methods with complex times for parabolic equations_. BIT 49 (2009) 487–508.
  * [14] S. Chin. _Symplectic integrators from composite operator factorizations_. Phys. Lett. A 226 (1997) 344–348.
  * [15] S. Chin. _Structure of positive decomposition of exponential operators_. Phys. Rev. E 71 (2005) 016703.
  * [16] I. Danaila, B. Protas. _Computation of ground states of the Gross–Pitaevskii functional via Riemannian optimization_. SIAM J. Sci. Comput. 39 (2017) 1102–1129.
  * [17] E. Gross. _Structure of a quantized vortex in boson systems_. Nuovo Cimento 20 (1961) 454–477.
  * [18] E. Hairer, Ch. Lubich, G. Wanner. _Geometric Numerical Integration_. Springer, Second Edition, 2006.
  * [19] E. Hansen, A. Ostermann. _High order splitting methods for analytic semigroups exist_. BIT 49 (2009) 527–542.
  * [20] E. Hille, R. Phillips. _Functional Analysis and Semi-Groups_. American Mathematical Society, 1974.
  * [21] A. Iserles. _A First Course in the Numerical Analysis of Differential Equations_. Cambridge University Press, Second Edition, 2008.
  * [22] E. Kieri. _Stiff convergence of force-gradient operator splitting methods_. Appl. Numer. Math. 94 (2015) 33–45.
  * [23] R. McLachlan, G. Quispel. _Splitting methods_. Acta Numerica 11 (2002) 341–434.
  * [24] A. Messiah. _Quantum Mechanics_. Dover, 1999.
  * [25] L. Pitaevskii. _Vortex lines in an imperfect Bose gas_. Sov. Phys. JETP 13 (1961) 451–454.
  * [26] 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.
  * [27] 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.
  * [28] R. Ruth. _A canonical integration technique_. IEEE Trans. Nucl. Sci. 30 (1983) 26–69.
  * [29] J. Sanz-Serna, M. Calvo. _Numerical Hamiltonian Problems_. Dover, 2018.
  * [30] Q. Sheng. _Solving linear partial differential equations by exponential splitting_. IMA J. Numer. Anal. 9 (1989) 199–212.
  * [31] M. Suzuki. _General theory of fractal path integrals with applications to many-body theories and statistical physics_. J. Math. Phys. 32 (1991) 400–407.
  * [32] M. Suzuki. _Hybrid exponential product formulas for unbounded operators with possible applications to Monte Carlo simulations_. Phys. Lett. A 201 (1995) 425–428.
  * [33] M. Thalhammer. _High-order exponential operator splitting methods for time-dependent Schrödinger equations_. SIAM J. Numer. Anal. 46 (2008) 2022–2038.
  * [34] M. Thalhammer, J. Abhau. _A numerical study of adaptive space and time discretisations for Gross–Pitaevskii equations_. J. Comput. Phys. 231 (2012) 6665–6681.
  * [35] V. Varadarajan. _Lie Groups, Lie Algebras, and Their Representations_. Springer, 1984.
  * [36] H. Yoshida. _Construction of higher order symplectic integrators_. Phys. Lett. A 150 (1990) 262–268.

Figure 1: Time integration of the Gross–Pitaevskii equation (2) involving a
quadratic potential by standard splitting methods and the novel modified
operator splitting method. Global errors versus time stepsizes in space
dimensions $d\in\\{1,2,3\\}$. Nonlinear ($\vartheta=1$) versus simplified
linear ($\vartheta=0$) case. Due to the validity of the invariance principle,
the application of an explicit Runge–Kutta method is not needed (RK0).

Figure 2: Time integration of the parabolic problem (3) involving a quadratic
potential by standard splitting methods and the novel modified operator
splitting method. Global errors versus time stepsizes in space dimensions
$d\in\\{1,2,3\\}$. Nonlinear ($\vartheta=1$) versus simplified linear
($\vartheta=0$) case. In order to resolve the nonlinear subproblem, a fourth-
order explicit Runge–Kutta method is applied (RK4). Depending on the stiffness
of the equation, stability is ensured for sufficiently small time stepsizes.
For a naive implementation of the Yoshida splitting method with complex
coefficients, an order reduction is observed, see Appendix B.

Figure 3: Corresponding results for the time-dependent Gross–Pitaevskii
equation (2) involving a fourth-order polynomial potential.

Figure 4: Corresponding results for the parabolic problem (3) involving a
fourth-order polynomial potential.

Figure 5: Time integration of the one-dimensional parabolic equation (3)
involving a quadratic potential (left) or a fourth-order potential (right),
respectively, by the modified operator splitting method (11). Global errors
versus time stepsizes. The original approach is based on the application of an
explicit fourth-order Runge–Kutta method for the numerical solution of the
nonlinear subproblem involving the double commutator (cf.
$\mathcal{E}_{\tau,\frac{2}{3}F_{2}-\frac{1}{72}\tau^{2}G_{2}}$). Alternative
approaches are based on the Strang splitting method (cf.
$\mathcal{E}_{\frac{1}{2}\tau,\frac{2}{3}F_{2}}\circ\mathcal{E}_{\tau,-\frac{1}{72}\tau^{2}G_{2}}\circ\mathcal{E}_{\frac{1}{2}\tau,\frac{2}{3}F_{2}}$).
Here, a reduced number of (inverse) fast Fourier transforms is required and an
improved accuracy is observed. Furthermore, the knowledge of the exact
solution to a component (cf. $\mathcal{E}_{\frac{1}{2}\tau,\frac{2}{3}F_{2}}$)
enhances the stability behaviour of the resulting time integration method for
larger time increments. Figure 6: Long-term integration of the one-dimensional
Gross–Pitaevskii equation (2) by standard and modified operator splitting
methods. Computation of numerical approximations to the values of the energy
at time grid points $t_{n}=n\,\tau$ for $\tau=10^{-3}$ and
$n\in\\{0,1,\dots,10^{5}\\}$ as well as corresponding deviations with respect
to the minimal values, see also (19). The obtained results confirm the
favourable geometric properties of the modified operator splitting method
(11).

##  Appendix A Matlab code

A Matlab code that illustrates the practical implementation of modified
operator splitting methods as outlined in Section 6 and reproduces the
numerical results displayed in Figures 1–4 is available at

doi.org/10.5281/zenodo.7945624.

By default, the space dimenson is set equal to one such that the overall
computation time amounts to a few minutes. For dimensions two and three,
respectively, the computational effort will increase accordingly to the
complexity of the problem.

Structure of the code. The underlying process carries out the time integration
of the Gross–Pitaevskii equation (2) and the related parabolic equation (3) by
the standard splitting methods with coefficiens specified in (8) and the novel
modified operator splitting methods (11) with iterated commutators given in
(14) and (15), respectively. The space discretisation is based on the Fourier
spectral method and realised by fast Fourier transforms and their inverses. In
loops over the selected space dimensions and the different test cases
(parabolic / Schrödinger equation, polynomial potential of degree two / four,
linear case $\vartheta=0$ / nonlinear case with $\vartheta=1$), the global
errors are computed for certain sequences of time stepsizes. In order to make
similarities apparent and contrast differences, several auxiliary functions
are defined.

  * •

_Core_. Selection of time integration methods. Definition of problem data
(initial and final time, initial state, potential and its derivatives). Choice
of sequences of time stepsizes. Performance of time integration. Computation
of global errors.

  * •

_PrecomputationFourier_. Computation of underlying space grid and eigenvalues
associated with Laplacian.

  * •

_FourierReal2Spectral, FourierSpectral2Real_. Fast Fourier transform and its
inverse.

  * •

_PartA_. Numerical solution of linear subproblems associated with Laplacian
based on Fourier transforms.

  * •

_BWithoutU, B, TimeStepRKM124B, PartB_. Numerical solution of nonlinear
subproblems in the context of standard operator splitting methods.

  * •

_DoubleCommutator, BModifiedWithoutU, BModified, TimeStepRKM124BModified,
PartBModified_. Numerical solution of nonlinear subproblems in the context of
modified operator splitting methods.

  * •

_TimeIntegration_. Time integration by standard and modified operator
splitting methods

  * •

_MyTestCases_. Definition of decisive quantities characterising different test
cases.

  * •

_MyPlot_. Visualisation of obtained results.

Specialisation and improvements. The main purpose of the elementary structured
Matlab code is the systematic comparison of the procedures required for the
novel schemes with those for standard operator splitting methods, on the one
hand for Schrödinger equations and on the other hand for parabolic equations.
We point out that a significant improvement of the performance will be
achieved by separating the different types of evolution equations as well as
time integrators and reconciling auxiliary functions. Further enhancements in
particular in connection with the more costly fast (inverse) Fourier
transforms concern the distinction of real and complex arithmetics as well as
the avoidance of redundancies, e.g., by reordering the eigenvalues associated
with the Laplacian instead of the solution values and omitting additional
scaling constants.

##  Appendix B Order reduction of complex splitting methods

In order to explain the observed order reduction for the complex splitting
method (8d), it suffices to study a nonlinear ordinary differential equation
of the form (1) with $F_{1}=0$ and $F_{2}(u)=|u|^{2}\,u$ on a single
subinterval of length $(0,\tau)$. More precisely, we consider the nonlinear
subproblems

$\begin{gathered}\begin{cases}\tfrac{\text{d}}{\text{d}t}\,u(t)=F_{2}^{(1)}\big{(}u(t)\big{)}=\big{(}u(t)\big{)}^{3}\,,\\\
u(0)=u_{0}\in\mathbb{R}\,,\quad
t\in(0,\tau)\,,\end{cases}\begin{cases}\tfrac{\text{d}}{\text{d}t}\,u(t)=F_{2}^{(2)}\big{(}u(t)\big{)}=|u(t)|^{2}\,u(t)\,,\\\
u(0)=u_{0}\in\mathbb{R}\,,\quad t\in(0,\tau)\,,\end{cases}\end{gathered}$

with coinciding real-valued solutions

$\mathcal{E}_{\tau,F_{2}^{(1)}}(u_{0})=u(\tau)=\mathcal{E}_{\tau,F_{2}^{(2)}}(u_{0})\,.$

On the one hand, using the Taylor series expansion

$u(\tau)=u(0)+\tau\,u^{\prime}(0)+\tfrac{1}{2}\,\tau^{2}\,u^{\prime\prime}(0)+\tfrac{1}{6}\,\tau^{3}\,u^{\prime\prime\prime}(0)+\tfrac{1}{24}\,\tau^{4}\,u^{\prime\prime\prime\prime}(0)+\mathcal{O}\big{(}\tau^{5}\big{)}$

in combination with the differential equation and derivatives thereof implies

$u(\tau)=\Big{(}1+\tau\,u_{0}^{2}+\tfrac{3}{2}\,\tau^{2}\,u_{0}^{4}+\tfrac{5}{2}\,\tau^{3}\,u_{0}^{6}+\tfrac{35}{8}\,\tau^{4}\,u_{0}^{8}\Big{)}\,u_{0}+\mathcal{O}\big{(}\tau^{5}\big{)}\,.$

As a consequence, for any real number $b\in\mathbb{R}$, the relation

$b\in\mathbb{R}:\quad
u(b\,\tau)=\Big{(}1+b\,\tau\,u_{0}^{2}+\tfrac{3}{2}\,b^{2}\,\tau^{2}\,u_{0}^{4}+\tfrac{5}{2}\,b^{3}\,\tau^{3}\,u_{0}^{6}+\tfrac{35}{8}\,b^{4}\,\tau^{4}\,u_{0}^{8}\Big{)}\,u_{0}+\mathcal{O}\big{(}\tau^{5}\big{)}$

is valid. For complex numbers, however, the more general expansion

$\begin{split}b,u_{0}\in\mathbb{C}:\quad\mathcal{E}_{\tau,bF_{2}^{(2)}}(u_{0})&=\Big{(}1+b\,\tau\,|u_{0}|^{2}+\big{(}b^{2}+\tfrac{1}{2}\,|b|^{2}\big{)}\,\tau^{2}\,|u_{0}|^{4}\\\
&\qquad+\big{(}b^{3}+(\tfrac{7}{6}\,b+\tfrac{1}{3}\,\overline{b})\,|b|^{2}\big{)}\,\tau^{3}\,|u_{0}|^{6}\\\
&\qquad+\big{(}b^{4}+\tfrac{1}{24}\,(46\,b^{2}+6\,\overline{b}^{2}+29\,|b|^{2})\,|b|^{2}\big{)}\,\tau^{4}\,|u_{0}|^{8}\Big{)}\,u_{0}\\\
&\qquad+\mathcal{O}\big{(}\tau^{5}\big{)}\end{split}$

is obtained by decomposing the solution and accordingly the defining function
into real and imaginary parts. In order to reproduce the approximation that
corresponds to the Yoshida splitting, we impose the basic symmetry and order
conditions

$b_{4}=b_{1}\,,\quad b_{3}=b_{2}\,,\quad b_{1}+b_{2}+b_{3}+b_{4}=1\,,$

and then perform the fourfold composition

$u_{\text{Splitting}}(\tau)=\big{(}\mathcal{E}_{\tau,b_{1}F_{2}^{(2)}}\circ\mathcal{E}_{\tau,b_{2}F_{2}^{(2)}}\circ\mathcal{E}_{\tau,b_{2}F_{2}^{(2)}}\circ\mathcal{E}_{\tau,b_{1}F_{2}^{(2)}}\big{)}(u_{0})\,.$

Requiring this expansion to coincide with the expansion of the exact solution

$u_{\text{Splitting}}(\tau)-u(\tau)=C(u_{0})\,\Im(b_{1})\,\big{(}1-\Re(b_{1})\big{)}\,\tau^{3}+\mathcal{O}\big{(}\tau^{4}\big{)}$

leads to a condition that obviously contradicts the order conditions for
splitting methods and explains the observation of local order three and global
order two. Similar arguments apply to evolution equations of Schrödinger type
such as the Gross–Pitaevskii equation. But, in this context, the invariance
principle permits to avoid the application of a fourth-order Runge–Kutta
method to the arising nonlinear subproblem

$\tfrac{\text{d}}{\text{d}t}\,u(t)=F_{2}^{(2)}\big{(}u(t)\big{)}=-\,\text{i}\,|u(t)|^{2}\,u(t)\,,\quad
t\in(0,\tau)\,.$


# Transmission properties of one dimensional metal and left-handed gratings

Gergely Kajtár<EMAIL_ADDRESS>Peter Markoš  Dept. of Physics, INPE
FEI STU, Ilkovičova 3, 812 19 Bratislava, Slovak Republic

###### Abstract

We provide rigorous numerical calculations (using the rigorous coupled-wave
analysis) of transmission spectra of various metallic and metamaterial one-
dimensional gratings and identify plasmonic resonances responsible for
enhanced transmission. We argue that the most important mechanism which
influences the resonant transmission is the coupling of incident
electromagnetic wave with two plasmonic waves: lengthwise plasmons, which
propagates along the grating, and crosswise plasmons excited in air gaps.

plasmons, metal and metamaterial grating, enhanced transmission, RCWA

###### pacs:

73.20.Mf, 78.66.Bz, 81.05.Xj, 41.20.Jb

††preprint: APS/123-QED

##  I INTRODUCTION

Surface plasmon economou ; maier  represents the basic eigenmode of a metallic
grating. Thanks to the spatial periodicity of the grating, surface plasmons
could be excited by incident electromagnetic (EM) wave. The interaction of
these two waves (surface plasmon and incident wave) strongly influences the
transmission properties of the periodic grating. Enhanced transmission of EM
waves through periodic metallic gratings has been experimentally observed in
Ref. ebbesen . Contrary, for very thin metallic gratings, transmission is
smaller due to the absorption of plasmon energy braun . Negative role of
surface plasmons is also discussed in Ref. lalanne .

In this paper, we calculate the transmission of EM waves through various one-
dimensional (1D) gratings and show that observed resonant frequency behavior
of the transmission can be explained as an interaction of incident wave with
excited plasmonic resonances. We distinguish two kinds of plasmons porto ; ros
. The first one, lengthwise (LW) plasmon, which propagates along the grating
surface and can be excited only when Bragg’s relation between the wave vector
of surface plasmon and grating period is fulfilled maier ; zayats . These LW
plasmons are also responsible for the extraordinary transmission through a
single aperture, if the incident surface is periodically corrugated lezec . LW
plasmon resonance appears as a sharp Fano-type peak in the spectrum. The
second one, crosswise (CW) plasmon, is excited inside the gaps takakura ; lee
. CW plasmon resonance appears as a broad Fabry-Perot type maximum in the
spectrum. We provide quantitative estimation of resonant frequencies for both
plasmons and show that they are in good agreement with rigorously calculated
transmission resonances.

The Paper is organized as follows: In Section II the transmittance of the EM
wave through a metallic grating, calculated with the use of the rigorous
coupled wave analysis (RCWA) rcwa , is presented. In order to reveal the
effect of grating geometry (grating thickness, grating period, width of gap),
transmittance is calculated for different grating parameters. In Section III,
observed resonances in the transmission spectra are interpreted as a coupling
effect of incident EM wave with the surface plasmons excited on the surface of
the structure or inside it. Section IV deals with the application of the model
to left-handed material gratings. Discussion of obtained results is given in
Section V.

##  II Transmission spectra

Figure 1:  Cross section of the 1D metallic grating. Incident EM wave has a
wave vector $\mathbf{k}$. Figure 2:  (color online) Reference transmission
spectrum of all diffracted orders versus the normalized frequency. Arrows
indicate enhanced transmissions of Fabry-Perot type (continuous) and Fano type
(dashed arrow) resonances. Red line shows the zeroth order spectrum, blue
lines represent higher order terms. Parameters of the metallic gratings are:
$p=1000$ nm, $d=300$ nm, $h=10$ nm (refer to Fig. 1). Note the change of the
transmission spectrum for $p/\lambda\sim 5.06$ and 7.16. At these frequencies
the real part of $\varepsilon_{m}$ becomes -1 and 0, respectively.

The reference 1D metallic grating is shown in Fig. 1. The grating is uniformly
extended in the $y$-direction. Its thickness in the $z$-direction is $d$ and
its spatial period in the $x$-direction is $p$. The width of air gap is $h$.
The metallic relative permittivity is $\varepsilon_{m}$. The embedding medium
is air with relative permittivity $\varepsilon_{d}=1$. The relative
permittivity of the metal is given by the Drude’s formula

$\varepsilon_{m}\left(f\right)=1-\frac{f_{p}^{2}}{f^{2}+\mathrm{i}f\gamma}$
(1)

where $f_{p}$ is the plasmon frequency and $\gamma$ is the absorption
coefficient. For our calculations $f_{p}=2147$ THz and $\gamma=5$ THz, typical
for Ag and Au metals .

Numerical simulations of the transmission were realized with our own numerical
program based on the RCWA. In the RCWA, the structure is divided into multiple
sandwiched layers along the direction of propagation. The permittivity is a
periodic function of $x$ in each sublayer. The numerical calculation involves
spatial Fourier expansion of the EM field and of the dielectric function in
each sublayer of the structure. The EM field determined by the RCWA satisfies
Maxwell’s equation within each sublayer as well as the boundary conditions
between adjacent layers. The numerical accuracy is only limited by the number
of orders used in the Fourier expansion. The typical number of Fourier modes
we used is 106. If not stated otherwise, the incident electromagnetic plane
wave is normal to the grating plane and the vector of magnetic field is
parallel to the $z$-direction (TM polarization).

Figure 2 shows the calculated transmission coefficient against the normalized
frequency $p/\lambda$. There are two types of enhanced transmission. The first
one is a Fabry-Perot resonance (five lowest resonances at $p/\lambda=0.65$,
1.30, 1.82, 2.32) are marked by arrows in Fig. 2. The second type, Fano
resonances fan , observed at $p/\lambda=0.99$, 1.92, 2.73, etc. are marked by
dashed arrows. Both resonances are present in all of the diffracted orders’
transmission.

To identify the physical origin of these resonances, we gradually change one
of the parameters of the structure and investigate how the transmission
spectrum changes compared to the reference one shown in Fig. 2.

Figure 3 shows how the transmission spectra changes when the angle of
incidence (in the $xz$ plane) increases from $0^{\circ}$ (this is the
reference curve from Fig. 2) to $2^{\circ}$, $5^{\circ}$ and $10^{\circ}$.
While the Fabry-Perot transmission maxima have not changed, Fano resonances
split to two separate frequencies which strongly depend on the incident angle.

Figure 3:  (color online) Transmission spectrum of the zeroth diffraction
order for various angles of incidence.

Transmission spectra for various grating thicknesses $d$ are shown in Fig. 4.
The position of the first Fabry-Perot type maximum increases for decreasing
thickness $d$. On the other hand, the position of the first Fano resonance
hardly depends on $d$. It is interesting to see how the Fabry-Perot peak
interacts with the Fano resonance. The “travelling” Fabry-Perot peak actually
transforms to Fano type peak and back. After further decreasing grating
thickness a transmission minimum appears on interval $0<p/\lambda<1)$.

Figure 4:  Transmission spectra for different grating thicknesses $d=(10,320)$
nm. Dashed line is the transmission for reference grating thickness $d=300$ nm
shown in Fig. 2.

Similarly, a changing air gap width $h$ affects the position of Fabry-Perot
type maximum (Fig. 5) while sharp Fano peaks remain unaffected.

Figure 5:  Transmission spectra for different air gap widths $h=5,10,15,20$
nm. Spectrum for reference air gap width $h=10$ nm is dashed.

Finally, we probe transmission against the grating period $p$. As it is seen
in Fig. 6, position of Fabry-Perot maxima is left unchanged and the sharp Fano
resonances are shifted to the right for decreasing grating period $p$ and
become less sharper.

Figure 6:  Transmission spectra for different grating periods. Curves are
normalized to reference period $p=1000$ nm.

##  III Theoretical framework

In this section we present an interpretation of the origin of enhanced
transmission phenomenon. We show that Fano-type sharp peaks are caused by the
interaction of incident EM wave with surface plasmons propagating along the
surface of the metallic grating. The excitation of these lengthwise (LW)
plasmons is possible thanks to the spatial periodicity of the grating.
Similarly, Fabry-Perot-type resonances are due to crosswise (CW) plasmons
excited and guided along the air gap resonators (Fig. 7).

Figure 7:  Three types of surface plasmons, namely lengthwise symmetric (A
left) and antisymmetric (A right) plasmon and crosswise (B) symmetric plasmon,
play role in enhanced transmission.

###  III.1 Lengthwise plasmons

Consider a homogeneous metallic slab with thickness $d$ and negative relative
permittivity $\varepsilon_{m}$ given by the real part of Eq. 1. The metal slab
is surrounded by dielectric relative permittivity $\varepsilon_{d}$. The
dispersion relation for the TM polarized surface waves excited on both sides
of the metal-dielectrics interface is given economou ; markos ; maier  as

$\frac{\varepsilon_{m}}{\varepsilon_{d}}=-\frac{\kappa^{m}_{z}}{\kappa^{d}_{z}}\tanh{\frac{\kappa_{z}^{m}d}{2}}$
(2)

for the symmetric plasmon, and

$\frac{\varepsilon_{m}}{\varepsilon_{d}}=-\frac{\kappa^{m}_{z}}{\kappa^{d}_{z}d}\coth{\frac{\kappa_{z}^{m}}{2}}$
(3)

for the antisymmetric one. Parameters $\kappa_{z}^{m}$ and $\kappa_{z}^{d}$
determine exponentially decreasing intensity of the EM field by the distance
from the metal-air interface:

$\left(\kappa_{z}^{m}\right)^{2}=k_{x}^{2}-\varepsilon_{m}k^{2}$ (4)
$\left(\kappa_{z}^{d}\right)^{2}=k_{x}^{2}-\varepsilon_{d}k^{2}$ (5)

Here, $k=\omega/c$ is the wave vector of the incident wave in vacuum and
$k_{x}$ is the $x$-component of the plasmon. Since $\varepsilon_{d}$ is
positive, Eq. 5 would be satisfied only if $k_{x}>k$ for $\varepsilon_{d}=1$.
This is not achievable solely by the incident $k_{x}=k\sin{\varphi}$
($\varphi$ is the angle of incidence), therefore the grating periodicity is
required to modulate $k_{x}$:

$k_{x}=k\sin{\varphi}+m\frac{2\pi}{p},~{}~{}~{}~{}m=\pm 1,\pm 2,\dots$ (6)

Using Eqs. (2-6) we numerically obtain the positions of lengthwise plasmon
resonances for different integers $m$ versus the slab thickness $d$. Results
are plotted in Fig. 8 for normal incidence.

As it follows from Eq. 6, the position of the LW resonance strongly depends on
the angle of incidence and period $p$. On the other hand, the dependence on
the thickness is weak, at least for sufficiently large $d$. However, air gap
width has some effect on this resonance since it changes the effective
relative permittivity of the metallic grating (Fig. 5).

Considering these observations of numerical data from transmission spectra
(Figs. 3, 6 and 4) we conclude that Fano-type sharp resonances are generated
by lengthwise plasmons propagating on interfaces along the grating. Evidently,
lengthwise plasmon resonance creates negative transmission effect in
conjunction with Ref. lalanne . In contrast to Ref. lalanne  we can not
confirm that these resonances are closely related to diffraction orders.
Although the first two sharp resonances are close to the positions of
diffracted orders (because the dispersion curve of the plasmon lies very close
to the light cone $\omega=ck$), higher order resonances are not connected to
higher diffracted orders (see Table 1).

###  III.2 Crosswise plasmons

Crosswise plasmon resonance propagates across the grating along the interfaces
of air gap (see Fig. 7). Corresponding equations are similar to that for LW
plasmons:

$\frac{\varepsilon_{d}}{\varepsilon_{m}}=-\frac{\kappa^{d}_{x}}{\kappa^{m}_{x}}\tanh{\frac{\kappa_{x}^{d}h}{2}}$
(7)

for the symmetric mode and

$\frac{\varepsilon_{d}}{\varepsilon_{m}}=-\frac{\kappa^{d}_{x}}{\kappa^{m}_{x}}\coth{\frac{\kappa_{x}^{d}h}{2}}$
(8)

for the antisymmetric mode. Again, $\kappa_{x}^{m}$ and $\kappa_{x}^{d}$
determine exponential decrease of the EM field by the distance from the air-
metal interface

$\left(\kappa_{x}^{m}\right)^{2}=k_{z}^{2}-\varepsilon_{m}k^{2}$ (9)
$\left(\kappa_{x}^{d}\right)^{2}=k_{z}^{2}-\varepsilon_{d}k^{2}$ (10)

The wave vector in the $z$ direction, $k_{z}$, is given by a Fabry-Perot
resonance condition

$k_{z}=\frac{\pi}{d}m,~{}~{}~{}~{}m=1,2,\dots$ (11)

Figure 8:  (color online) Plasmon resonances for different grating thicknesses
$d$. There are five orders ($m=1,2,3,4,5$) showed for CW plasmons and three
orders ($m=1,2,3$) for LW plasmons, respectively. Symbols represent
numerically obtained position of the first Fabry-Perot resonance (Fig. 4)
which is very close to the resonant frequency of the $m=1$ CW plasmon.

Accordingly, the wave is guided through the air gap if $k_{z}$ is supported by
the gap-resonator’s length $d$. Solving Eq. 7-11 resonant frequencies are
obtained versus the grating thickness $d$ (Fig. 8). Apparently, there is no
resonance for the antisymmetric mode for the displayed interval
$0<p/\lambda<3$.

Crosswise plasmon resonance does not depend on the incident angle $\varphi$
and grating period $p$ but it depends on the thickness $d$ and gap width $h$
(refer to Figures in Section II and Fig. 9). This supports the idea that
Fabry-Perot type resonances are due to plasmons propagating inside air gaps.

###  III.3 Identifying plasmon resonances

We use the presented model for the identification of numerically calculated
resonances. Calculated resonant frequencies, summarized in Tables 1 and 2
confirm that there is a very good agreement between rigorous calculations and
the model.

Figure 9:  Symmetric crosswise plasmon resonance versus air gap width $h$ for the first three orders ($m=1,2,3$) starting from the bottom curve. Table 1:  Comparison of calculated values using RCWA simulation from Fig. 2 and the theoretical model from Fig. 8.  m | lengthwise | RCWA | crosswise | RCWA
---|---|---|---|---
1 | 0.99 | 0.99 | 0.70 | 0.65
2 | 1.92 | 1.92 | 1.36 | 1.29
3 | 2.73 | 2.73 | 1.93 |  1.82111Can not be clearly determined because it overlaps a lengthwise plasmon resonance.
4 | 3.33 | 3.38 | 2.34 | 2.32
5 | 3.86 | 3.85 | 2.79 |  2.701
Table 2:  Comparison of resonant frequencies calculated with the RCWA simulations (Fig. 3) to the prediction of the model. $\varphi$ | m | RCWA | lengthwise plasmon222CW plasmon resonance is not dependent on $\varphi$.
---|---|---|---
0 | 1 | 0.99 | 0.99
$2^{\circ}$ | 1 | 1.03 | 1.02
$2^{\circ}$ | $-1$ | 0.96 | 0.96
$5^{\circ}$ | 1 | 1.08 | 1.08
$5^{\circ}$ | $-1$ | 0.91 | 0.91
$10^{\circ}$ | 1 | 1.19 | 1.19
$10^{\circ}$ | $-1$ | 0.85 | 0.85

For grating thickness $d=300$ nm symmetric and antisymmetric LW plasmon
resonances are barely distinguishable. To test the theory, we calculate the
transmission spectra for thickness $d=39$ nm (Fig. 10). Now, the first ($m=1$)
symmetric and antisymmetric LW plasmon should appear at 0.98 and 1.00,
respectively. The second and the third antisymmetric and symmetric modes are
clearly distinguishable at 1.86 and 1.96 ($m=2$) and 2.58 and 2.84 ($m=3$).
This agrees with observed positions of antisymmetric LW resonances (1.00, 1.90
and 2.62) and symmetric LW resonances (1.96 and 2.80). Apparently, the first
symmetric resonance is suppressed by the first antisymmetric resonance.

Apart from plasmonic resonances, we observe a transmission dip at wavelength
$p/\lambda\sim 0.68$. The position of this dip depends on the width of the air
gap and on the thickness of the grating. The incoming EM wave stops to
interact with the metallic grating for frequencies below this dip, thus the
transmission increases to unity for the limit $\lambda\to\infty$.

Figure 10:  Transmission spectrum for grating thickness $d=39$ nm. Symmetric
(dashed arrow) and antisymmetric (continuous arrow) LW plasmon resonances are
visible. Dashed curve represents the transmission through homogeneous metallic
slab.

##  IV Left-handed material gratings

In this Section we apply the previous theory to the transmission of EM wave
through periodic grating made of left-handed metamaterial (LHM). We
characterize the metamaterial by the frequency dependent permittivity and
magnetic permeability using the model described in markos . The permittivity
is given by Eq. (1) but with plasma frequency $f_{p}=10$ GHz and absorption
$\gamma=0.01$ GHz. Frequency dependent permeability is given by the relation

$\mu(f)=1-\frac{Ff^{2}}{f^{2}-f_{0}^{2}}$ (12)

used for the description of the effective permeability of periodic array of
split-ring resonators srr . The resonant magnetic frequency $f_{0}=4$ GHz and
filling factor $F=0.56$. The grating possesses a transmission band within the
frequency interval

$f_{0}~{}\leq~{}f~{}\leq~{}\displaystyle{\frac{f_{0}}{\sqrt{1-F}}}\approx
6~{}\textrm{GHz}.$ (13)

Bellow and above this transmission band, the metamaterial behaves as a metal
with frequency dependent positive magnetic permeability.

Figure 11 shows the dispersion curves of three surface plasmons propagating
along the interface LHM - air ruppin ; markos . Since the permeability of the
LHM is negative on frequency interval $(0.4,0.6)$, a TE polarized surface wave
might also be excited.

Figure 11:  (color online) Thick solid lines show dispersion curves of surface
plasmons propagating along the air - LHM interface. Note the logarithmic scale
on the horizontal axis. For a thin air gap ($h=0.3$ mm), each plasmon splits
to a symmetric (dashed lines) and an antisymmetric (dot-dashed lines) plasmon.
Note that some dispersion curves lie very close to the LHM light cone.  Figure
12:  (color online) Transmission of EM wave through the LHM grating. Black
solid line: TM polarization. Red dashed line: TE polarization. Parameters of
the structure are $p=3$ cm, $d=3$ cm, air gap width $h=p/100$. Insets show
details of the transmission in the vicinity of the lengthwise TM plasmon
resonances.

Figure 12 presents the frequency dependent transmission through the LHM
grating. As expected, transmission is close to unity within the LHM
transmission band (Eq. 13). Nevertheless, we observed strong Fano resonances
here (not clearly visible in the Fig. 12) due to the coupling of incident EM
wave with LW guided waves shadrivov  propagating along the air-LHM-air
structure in the $x$-direction and $z$-direction. For the TM modes, the
dispersion curves are

$\mu(f)=\frac{k_{z}^{\rm
LHM}}{\kappa_{z}^{d}}\times\left\\{\begin{array}[]{l}\tan{k_{z}^{\rm
LHM}p}/{2}\\\\[8.53581pt] \cot{k_{z}^{\rm LHM}p}/{2}\end{array}\right.$ (14)

for the symmetric and antisymmetric modes, respectively. Here,

$(k_{z}^{\rm LHM})^{2}=\varepsilon(f)\mu(f)k^{2}-k_{x}^{2}$ (15)

and

$(\kappa_{z}^{d})^{2}=k_{x}^{2}-\varepsilon k^{2}$ (16)

for the LHM and air, respectively. The value of $k_{x}$ is given by Eq. 6.

On the frequency region $f<f_{0}$ we observe a series of transmission maxima
which could be identified as Fabry-Perot resonances based on the excitation of
CW plasmons inside air gaps. Since the thickness of the grating is rather
large ($d=p$), we observe large number of CW plasmons (Fig. 13). These
plasmonic resonances exist also in the transmission band (Eq. 13). Above the
transmission band, CW plasmons cause weak oscillations of the transmission
(not visible in Fig. 12).

For the TM polarization, we identify two LW plasmon resonances at wavelengths
$p/\lambda=0.396$ and $0.65$. The latter is less pronounced since the
exponential decay of plasmon is very slow $d\kappa_{\rm LHM}\sim 1$. We have
not succeeded to identify the TE LW plasmon resonance, expected at
$p/\lambda\approx 0.485$, because the transmission of the TE wave also
fluctuates in the LHM band (Eq. 13) due to excitations of CW plasmons.

Figure 13:  (color online) Resonant frequencies for crosswise plasmon resonances for LHM grating with period $p=3$ cm and air gap width $h=p/100$ as a function of the grating thickness $d$. The wave vector of LW plasmon is $k_{z}=(\pi/d)m$. CW plasmons exist also in the LHM transmission band (Eq. 13) and are responsible for fluctuations of the transmission curve.  Table 3:  Comparison of calculated values using RCWA from Fig. 12 and the theory from Fig. 13 for LHM grating.  m | crosswise plasmon resonance | RCWA
---|---|---
1 | 0.08 | 0.08
2 | 0.14 | 0.14
3 | 0.19 | 0.19
4 | 0.23 | 0.22
5 | 0.26 | 0.25

##  V Conclusions

We have shown that enhanced transmission phenomena for periodic metallic
gratings are caused by two clearly distinguishable plasmon resonances, namely
lengthwise and crosswise plasmon resonances. Lengthwise plasmons (symmetric or
antisymmetric) propagate along the two interfaces induced by the grating
period. Crosswise plasmons propagate inside the air gaps. Only symmetric CW
plasmons were observed. In the metallic gratings, the frequency of
antisymmetric CW plasmons is much higher for a narrow air gap, and in the LHM
gratings, antisymmetric CW plasmons lie very close to the LHM transmission
band (Fig. 11).

For completeness, we note that another excitations exist in the air gap,
namely the guiding modes with both components of the wave vector, $k_{x}$ and
$k_{z}$, real inside the gap. However, similarly to antisymmetric CW plasmons
in metallic gratings, the frequency of these excitations is not accessible
neither for metallic nor for the LHM gratings.

Our analysis does not consider the role of the absorption, which is given by
the small absorption parameter $\gamma$ in the frequency-dependent
permittivity. Good agreement of transmission data with the prediction of the
model indicates that for these parameters the transmission does not play
significant role in the propagation of EM wave.

We have shown that plasmon resonance does not rely on diffracted orders but
solely on the dispersion relation of propagating plasmons. Our rigorous
calculations were used to prove the validity of a simple model of propagating
plasmons. This simple model can be used to identify the type of enhanced
transmission phenomenon (lengthwise or crosswise plasmon resonance) and to
design various metal filters with subwavelength features which support
enhanced transmission.

This work was supported by the Slovak Research and Development Agency under
the contract No. APVV-0108-11 and VEGA Agency, project No. 1/0372/13.

## References

  * (1) E. N. Economou, Phys. Rev. 182, 569 (1969);
  * (2) S. A. Maier, Plasmonics: Fundamentals and Applications Springer Science (2002).
  * (3) T. W. Ebbesen, H. J. Lezec, H.F. G haemi, T. Thio, P. A. Wolff, Nature 391, 667 (1998); H. F. Ghaemi, T. Thio, D. E. Grupp, T. W. Ebbesen and H. J. Lezec, Phys. Rev. B 58, 6779 (1998).
  * (4) J. Braun, B. Grompf, G. Kobiela, M. Dressel, Phys. Rev. Lett. 103, 203901 (2009).
  * (5) Q. Cao and Ph. Lalanne, Phys. Rev. Lett. 88, 057403 (2002).
  * (6) J. A. Porto, J. García-Vidal, J. B. Pendry, Phys. Rev. Lett. 83, 2845 (1999).
  * (7) A. Roskiewicz and W. Nasalski, J. Phys. B: At. Mol. Opt. Phys.  46, 025401 (2013).
  * (8) A. V. Zayats, I. I. Smolyaninov and A. A. Maradudin, Phys. Rep. 408, 131 (2005).
  * (9) H. J. Lezec et al., Science 297, 820 (2002); L. Martin-Moreno, F. J. García-Vidal, H. J. Lezec, A. Degiron, and T. W. Ebbesen, Phys. Rev. Lett. 90, 167401 (2003).
  * (10) Y. Takakura, Phys. Rev. Lett. 86, 5601 (2001).
  * (11) K. G. Lee and Q-Han Park, Phys. Rev. Lett. 98, 103902 (2005).
  * (12) J. J. Hench and Z. Strakoš, Electron. Trans. Num. Anal. 31, 331 (2008).
  * (13) I. El-Kady et al., Phys Rev. B 62, 15299 (2000).
  * (14) Fan, Shanhui, Wonjuoo Suh, and J. D. Joannopoulos, J. Opt. Soc. Am. A 20, 569–572 (2003).
  * (15) P. Markoš and C. M. Soukoulis, Wave Propagation (Princeton University Press, Princeton, 2008).
  * (16) J. B. Pendry et al., IEEE Trans. Microwave Theory Tech. 47, 2075 (1999).
  * (17) N. Engheta and R. W. Ziolkowski (Eds.) Metamaterials Willey-Interscience (2006).
  * (18) R. Ruppin, Phys. Lett. A 277, 63 (2000); J. Phys.: Condens. Matt. 13, 1811 (2001).
  * (19) I. V. Shadrivov, A. A. Sukhorukov, and Yu. S. Kivshar, Phys. Rev. E  67, 057602 (2003).