# An efficient spectral method for the fractional Schrödinger equation on the
real line

Mengxia Shen111School of Mathematics and Statistics, Huazhong University of
Science and Technology, Wuhan 430074, P. R. China.
<EMAIL_ADDRESS>and Haiyong Wang11footnotemark: 1 222Hubei Key
Laboratory of Engineering Modeling and Scientific Computing, Huazhong
University of Science and Technology, Wuhan 430074, P. R. China

###### Abstract

The fractional Schrödinger equation (FSE) on the real line arises in a broad
range of physical settings and their numerical simulation is challenging due
to the nonlocal nature and the power law decay of the solution at infinity. In
this paper, we propose a new spectral discretization scheme for the FSE in
space based upon Malmquist-Takenaka functions. We show that this new
discretization scheme achieves much better performance than existing
discretization schemes in the case where the underlying FSE involves the
square root of the Laplacian, while in other cases it also exhibits comparable
or even better performance. Numerical experiments are provided to illustrate
the effectiveness of the proposed method.

Keywords: Fractional Laplacian, fractional Schrödinger equation, Malmquist-
Takenaka functions, spectral Galerkin method

AMS classifications: 65M70, 41A20

##  1 Introduction

Nonlocal and fractional models have become increasingly important because of
their connection with many real-world phenomena that appear in physics,
biology and materials science. In contrast to local models, nonlocal and
fractional models are more suitable to model complex systems exhibiting
singularities and anomalies as well as involving nonlocal interactions (see,
e.g., [10]). The fractional Laplacian operator, which can be seen as the
infinitesimal generator of a symmetric $\alpha$-stable Lévy process in
probability theory, is one of the fundamental nonlocal operators and it arises
in a number of applications such as anomalous diffusion, image denoising,
finance. The fractional Laplacian operator of a function
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is defined by

$(-\Delta)^{\alpha/2}f(x):=\alpha\frac{{2}^{\alpha-1}\Gamma(\frac{\alpha+d}{2})}{\pi^{d/2}\Gamma{(\frac{2-\alpha}{2}})}~{}\mathrm{p.v.}\int_{\mathbb{R}^{d}}\frac{f(x)-f(y)}{|x-y|^{d+\alpha}}\mathrm{d}y,\quad
x\in{\mathbb{R}^{d}},$ (1.1)

where $\alpha\in{(0,2)}$ and $\mathrm{p.v.}$ stands for the Cauchy principle
value. Equivalently, it can also be defined as a pseudo-differential operator
via the Fourier transform

$(-\Delta)^{\alpha/2}f(x):=\mathcal{F}^{-1}[|\xi|^{\alpha}\mathcal{F}[f](\xi)](x),$
(1.2)

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the Fourier and inverse
Fourier transforms, respectively. It is well known that the fractional
Laplacian operator reduces to the identity operator whenever
$\alpha\rightarrow 0$ and to the negative Laplacian operator whenever
$\alpha\rightarrow 2$. When discretizing nonlocal models involving the
fractional Laplacian operators, the main difficulty stems from the nonlocal
and singular nature of the fractional Laplacian operators and the slow decay
of the underlying solution at infinity. In the past decade, numerical methods
for such nonlocal models have attracted a lot of attention and many
significant advances have been made, such as the finite element methods [1,
5], the finite difference methods [12, 17, 28] and spectral methods [2, 6, 27,
30, 31, 34]. Among these methods, the finite element and finite difference
methods are typically studied for nonlocal models on bounded domains. While
for nonlocal models on unbounded domains, spectral methods are particularly
attractive due to their global character.

In this paper, we are interested in the fractional Schrödinger equation (FSE)
on the real line

$\displaystyle\begin{cases}\mathrm{i}\partial_{t}\psi(x,t)=\gamma(-\Delta)^{\alpha/2}\psi(x,t)+\mathcal{T}\psi(x,t),\quad
x\in\mathbb{R},\quad t>0,\\\\[5.0pt] \psi(x,0)=\psi_{0}(x),\quad
x\in\mathbb{R},\\\\[3.0pt]
\lim_{|x|\rightarrow{\infty}}\psi(x,t)=0,\end{cases}$ (1.3)

where $\mathrm{i}=\sqrt{-1}$, $\gamma\in\mathbb{R}$ and $\gamma\neq 0$,
$\mathcal{T}$ is a linear or nonlinear operator (e.g.,
$\mathcal{T}\psi(x,t)=V(x)\psi(x,t)$ or
$\mathcal{T}\psi(x,t)=\pm|\psi(x,t)|^{2}\psi(x,t)$) and $\psi(x,t)$ is a
complex-valued wave function. Equation (1.3), which was introduced by Laskin
in [25], is a natural generalization of the standard Schrödinger equation that
arises in the context of the well-known Feynman path integrals approach to
quantum mechanics when the Brownian trajectories are replaced by Lévy flights.
More recently, experimental realizations and rigorous derivations of FSE have
been widely investigated in different branches of physics, such as the
continuum limit of certain discrete physical systems with long-range
interactions [22], beam propagation [16, 26] and the Lévy crystal in a
condensed-matter environment [32]. Around the same time, numerical methods
that combine spectral methods for the spatial discretization and time-stepping
methods for the temporal discretization for the solution of (1.3) have
attracted considerable attention (see, e.g., [11, 23, 27, 30, 31]). However,
when the solution decays slowly with a power law at infinity, exponential
convergence of these existing spectral discretization schemes using either
Fourier, Hermite or mapped Chebyshev functions in space has so far not been
observed. Here, we propose a novel spectral discretization scheme using
Malmquist-Takenaka functions which have excellent approximation properties for
functions with poles in the complex plane. We show that this new
discretization method achieves exponential convergence in the particular case
of $\alpha=1$, regardless of the underlying FSE is linear or nonlinear, and
exhibits a comparable or even better performance than state-of-the-art
discretization schemes in other cases.

The rest of this paper is organized as follows. In section 2, we briefly
review some properties of Malmquist-Takenaka functions that will be used for
the construction of our spectral discretization method. In section 3, we
present a numerical method for FSE by combining a spectral Galerkin method
using Malmquist-Takenaka functions for the space discretization coupled with
time-stepping schemes for the temporal discretization. We perform two
numerical examples to illustrate the performance of the proposed method in
section 4 and conclude the paper with some remarks in section 5.

##  2 Malmquist-Takenaka functions

The Malmquist-Takenaka functions***Up to some constant and scaling factors,
the Malmquist-Takenaka functions are also known as the Christov functions in
some literature. Here we adopt the name used in [18]. (MTFs) are defined by

$\varphi_{n}(x)=\mathrm{i}^{n}\sqrt{\frac{2}{\pi}}\frac{(1+2\mathrm{i}x)^{n}}{(1-2\mathrm{i}x)^{n+1}},\quad
n\in{\mathbb{Z}}.$ (2.1)

Let $L^{2}(\mathbb{R})$ denote the space of square integrable functions and
let $(\cdot,\cdot)$ denote the inner product defined by
$(f,g)=\int_{\mathbb{R}}f(x)\overline{g(x)}\mathrm{d}x$, where
$\overline{f(x)}$ denotes the complex conjugate of $f(x)$. It is well known
that the system $\\{\varphi_{n}\\}_{n\in\mathbb{Z}}$ forms a complete and
orthonormal basis in $L^{2}(\mathbb{R})$, i.e.,

$(\varphi_{n},\varphi_{m})=\delta_{n,m},$ (2.2)

where $\delta_{n,m}$ is the Kronecker delta. Theoretical aspects of MTFs as
well as their applications in designing algorithms for Fourier and Hilbert
transforms have been investigated during the past few decades (see, e.g., [8,
15, 18, 19, 35, 36]). We list below some of theoretical properties of MTFs
that will be used for the construction of spectral method.

  * (i)

They satisfy the following differential recurrence relation

$\frac{\mathrm{d}}{\mathrm{d}x}\varphi_{n}(x)=-n\varphi_{n-1}(x)+\mathrm{i}(2n+1)\varphi_{n}(x)+(n+1)\varphi_{n+1}(x),$
(2.3)

from which we can deduce that the differentiation matrix of the Malmquist-
Takenaka system is skew-Hermitian and tridiagonal. Moreover, higher order
derivatives of $\varphi_{n}(x)$ can be derived from repeated application of
(2.3). Moreover, they also satisfy

$x\frac{\mathrm{d}}{\mathrm{d}x}\varphi_{n}(x)=-\frac{n}{2}\mathrm{i}\varphi_{n-1}(x)-\frac{1}{2}\varphi_{n}(x)-\frac{n+1}{2}\mathrm{i}\varphi_{n+1}(x).$
(2.4)

  * (ii)

They are uniformly bounded on the real line and

$|\varphi_{n}(x)|=\sqrt{\frac{2}{\pi}}\frac{1}{\sqrt{1+4x^{2}}}\leq\sqrt{\frac{2}{\pi}},\quad
x\in\mathbb{R}.$ (2.5)

  * (iii)

For $N\in\mathbb{N}$, we define the space
$\mathbb{V}_{N}(\mathbb{R})=\mathrm{span}\\{\varphi_{k}(x)\\}_{k=-N}^{N-1}$
and denote by $\Pi_{N}:L^{2}(\mathbb{R})\rightarrow\mathbb{V}_{N}(\mathbb{R})$
the orthogonal projection operator, i.e.,

$(\Pi_{N}f)(x)=\sum_{k=-N}^{N-1}a_{k}\varphi_{k}(x),\quad
a_{k}=(f,\varphi_{k}).$ (2.6)

With the change of variable $x=\tan(\theta/2)/2$, we obtain

$\displaystyle a_{k}$
$\displaystyle=\frac{(-\mathrm{i})^{k}}{2\sqrt{2\pi}}\int_{-\pi}^{\pi}f\left(\frac{1}{2}\tan\frac{\theta}{2}\right)\left(1-\mathrm{i}\tan\frac{\theta}{2}\right)e^{-\mathrm{i}k\theta}\mathrm{d}\theta$
$\displaystyle\approx\frac{\mathrm{(-i)}^{k}}{2N}\sqrt{\frac{\pi}{2}}\sum_{j=0}^{2N-1}f\left(\frac{1}{2}\tan\frac{\theta_{j}}{2}\right)\left(1-\mathrm{i}\tan\frac{\theta_{j}}{2}\right)e^{-\mathrm{i}k\theta_{j}},$
(2.7)

where $\theta_{j}=-\pi+{\pi}j/N$ and we have used the fact that integrals of
periodic functions can be computed efficiently by using the composite
trapezoidal rule. Hence, the computation of $\\{a_{k}\\}_{k=-N}^{N-1}$ can be
performed rapidly with a single FFT in $\mathcal{O}(N\log N)$ operations.

  * (iv)

They are eigenfunctions of the Hilbert transform

$\mathcal{H}[\varphi_{n}](x)=(-\mathrm{i})\mathrm{sgn}(n)\varphi_{n}(x),\quad\mathcal{H}[f](x)=\frac{1}{\pi}\mathrm{p.v.}\int_{\mathbb{R}}\frac{f(z)}{x-z}\mathrm{d}z,$
(2.8)

where $\mathrm{sgn}(n)=1$ for $n=0,1,\ldots$ and $\mathrm{sgn}(n)=-1$ for
$n=-1,-2,\ldots$. Moreover, the Hilbert transform is related to the square
root of the Laplacian by $(-\Delta)^{1/2}f(x)=\mathcal{H}[f{{}^{\prime}}](x)$.

Figure 1: The magnitude of the Malmquist-Takenaka coefficients of
$\\{a_{k}\\}_{k\in\mathbb{Z}}$. Left: $f(x)=\exp(-x^{2})$ and
$\mathrm{sech}(x)$. Right: $f(x)=1/(1+x^{2})$ and $1/(1+x^{4})$.

Having introduced basic properties of MTFs, we now move on to their spectral
approximation properties. To gain insight, we plot in Figure 1 the absolute
value of the Malmquist-Takenaka coefficients of
$f(x)=\exp(-x^{2}),\mathrm{sech}(x),1/(1+x^{2}),1/(1+x^{4})$. Clearly, we
observe that the coefficients of the first two functions decay at
subexponential rates and the coefficients of the last two functions decay at
exponential rates. Indeed, Weideman in [36] had studied the asymptotic
estimate of the coefficients $\\{a_{k}\\}_{k\in\mathbb{Z}}$ of several
representative functions, including these four test functions, and his
analysis provides an important insight to understand the approximation power
of MTFs. Below we present a simplified version of Weideman’s result on the
exponential convergence of Malmquist-Takenaka approximation and we provide a
short proof for the purpose of being self-contained.

######  Theorem 2.1.

Let $\Pi_{N}f$ denote the spectral approximation defined in (2.6) and let
$\mathcal{A}_{\rho}$ denote the annulus defined by
$\mathcal{A}_{\rho}=\\{z\in\mathbb{C}:\rho^{-1}<|z|<\rho\\}$. Moreover, let

$\hat{f}(z)=\frac{2}{1+z}f\left(\frac{z-1}{2\mathrm{i}(z+1)}\right).$ (2.9)

If $\hat{f}(z)$ is analytic in $\mathcal{A}_{\rho}$ for some $\rho>1$. Then
for $x\in\mathbb{R}$,

$\displaystyle\|f-\Pi_{N}f\|_{L^{\infty}(\mathbb{R})}=\mathcal{O}(\rho^{-N}).$
(2.10)

###### Proof.

By using the map $z=(1+2\mathrm{i}x)/(1-2\mathrm{i}x)$, we have

$\displaystyle a_{k}$
$\displaystyle=(f,\varphi_{k})=\frac{(-\mathrm{i})^{k+1}}{\sqrt{2\pi}}\oint_{|z|=1}\hat{f}(z)z^{-k}\mathrm{d}z,$

which implies that $a_{k}$ is also the Laurent coefficient of $\hat{f}(z)$.
Furthermore, since $\hat{f}(z)$ is analytic in the annulus
$\mathcal{A}_{\rho}$, it follows that $a_{k}=\mathcal{O}(\rho^{-|k|})$ for all
$k\in\mathbb{Z}$. Consequently,

$\displaystyle\|f-\Pi_{N}f\|_{L^{\infty}(\mathbb{R})}\leq\sum_{k=N}^{\infty}|a_{k}|\|\varphi_{k}\|_{L^{\infty}(\mathbb{R})}+\sum_{k=-\infty}^{-N-1}|a_{k}|\|\varphi_{k}\|_{L^{\infty}(\mathbb{R})},$

and the desired result (2.10) follows by combining the above inequality with
(2.5). ∎

A few remarks on the approximation power of MTFs are in order.

######  Remark 2.2.

If $f(x)$ has a partial fraction decomposition of the form

$f(x)=\sum_{k=1}^{m}\sum_{j=1}^{\chi_{k}}\frac{\eta_{k,j}}{(x-s_{k})^{j}},$
(2.11)

where $\\{s_{k}\\}_{k=1}^{m}$ is a set of poles in the complex plane but not
on $\mathbb{R}$ and $\\{\chi_{k}\\}_{k=1}^{m}$ is a set of orders associated
with those poles, then the above theorem indicates that the Malmquist-Takenaka
projection $\Pi_{N}f$ converges at an exponential rate.

######  Remark 2.3.

Recently, Iserles, Luong and Webb in [19] compared the approximation power of
the Malmquist-Takenaka, Hermite and stretched Fourier functions for Gaussian
wave packet functions of the form
$f(x)=\exp(-\beta(x-x_{0})^{2})\cos(\omega{x})$, where $\beta>0$ and
$x_{0},\omega\in\mathbb{R}$. After some lengthy algebra, they derived the
decay rates of the coefficients with respect to these three orthogonal
systems, respectively, and concluded that the MTFs are superior to the other
two functions. Note that all those three functions have banded skew-Hermitian
differentiation matrices.

######  Remark 2.4.

Mapped Chebyshev functions (MCFs) were recently used to develop spectral
methods for PDEs with fractional Laplacian on the unbounded domain [30].
Specifically, the MCF is defined by

$\mathbb{T}_{k}(x)=\frac{1}{\sqrt{c_{k}\pi/2}}\frac{1}{\sqrt{1+x^{2}}}T_{k}\left(\frac{x}{\sqrt{1+x^{2}}}\right),\quad
k=0,1,\ldots.$

where $c_{0}=2$ and $c_{k}=1$ for $k\geq 1$ and $T_{k}(x)$ is the Chebyshev
polynomial of the first kind of degree $k$. It is easy to verify that
$\\{\mathbb{T}_{k}(x)\\}_{k=0}^{\infty}$ forms an orthonormal system on
$\mathbb{R}$. Let $S_{N}^{\mathrm{MC}}(x)$ denote the MCF approximation of the
form

$S_{N}^{\mathrm{MC}}(x)=\sum_{k=0}^{2N-1}a_{k}^{\mathrm{MC}}\mathbb{T}_{k}(x),\quad
a_{k}^{\mathrm{MC}}=\int_{\mathbb{R}}f(x)\mathbb{T}_{k}(x)\mathrm{d}x.$ (2.12)

Concerning MTF and MCF approximations, it is natural to ask which one is
better. Note that both $\Pi_{N}f$ and $S_{N}^{\mathrm{MC}}(x)$ have the same
number of terms. Figure 2 displays the maximum errors of $\Pi_{N}(f)$ and
$S_{N}^{\mathrm{MC}}(x)$ for $f(x)=\exp(-x^{2})/(1+\mathrm{i}x)$,
$\exp(-x^{2})$, $1/(x^{2}+4)$, $\mathrm{sech}(x)$, $1/(x^{4}+1)$ and
$1/(x^{4}+1)^{1.2}$. Clearly, we can see that both approximations have their
own advantages. A systematic study on the comparison of both approximations is
beyond the scope of the present paper and will be discussed in another
publication.

Figure 2: Maximum errors of MTF and MCF approximations for
$f(x)=\exp(-x^{2})/(1+\mathrm{i}x)$, $1/(x^{2}+4)$, $1/(x^{4}+1)$ (left
column) and $f(x)=\exp(-x^{2})$, $\mathrm{sech}(x)$, $1/(x^{4}+1)^{1.2}$
(right column).

######  Remark 2.5.

In practical calculations, it is beneficial to introduce a scaling parameter
in MTFs. More specifically, the scaled MTFs are defined by

$\varphi_{n}^{\mathrm{S}}(x)=\nu^{-1/2}\varphi_{n}\left(\frac{x}{\nu}\right),\quad
n\in\mathbb{Z},$ (2.13)

where $\nu>0$ is a scaling parameter. It is easy to check that
$\\{\varphi_{n}^{\mathrm{S}}\\}_{n\in\mathbb{Z}}$ forms an orthonormal
sequence on $\mathbb{R}$. In the rest of this paper, we will use the scaled
Malmquist-Takenaka approximations when we mention the scaling parameter $\nu$
explicitly.

The generalized Laguerre polynomials (GLPs), denoted by $L_{n}^{(\alpha)}(x)$
with $\alpha>-1$, are orthogonal with respect to the weight function
$\omega_{\alpha}(x)=x^{\alpha}e^{-x}$ on the half line
$\mathbb{R}_{+}:=(0,+\infty)$ and

$\int_{\mathbb{R}_{+}}L_{n}^{(\alpha)}(x)L_{m}^{(\alpha)}(x)\omega_{\alpha}(x)\mathrm{d}x=\gamma_{n}^{(\alpha)}\delta_{n,m},$
(2.14)

where $\gamma_{n}^{(\alpha)}=\Gamma(n+\alpha+1)/\Gamma(n+1)$. In particular,
we will drop the superscript in $L_{n}^{(\alpha)}(x)$ whenever $\alpha=0$,
i.e., $L_{n}^{(0)}(x)=L_{n}(x)$. With the Laguerre polynomials we introduce a
sequence of functions on the real line for all $n\in\mathbb{Z}$ (see [21, 35])

$\Psi_{n}(x)=\left\\{\begin{array}[]{ll}e^{-x/2}L_{n}(x)H(x),&n\geq
0,\\\\[8.0pt] -e^{x/2}L_{-n-1}(-x)H(-x),&n<0,\end{array}\right.$ (2.15)

where $H(x)$ is the Heaviside step function. It is easy to check that
$\Psi_{n}(x)=-\Psi_{-n-1}(-x)$ for all $n\in\mathbb{Z}$ and
$\\{\Psi_{n}\\}_{n\in\mathbb{Z}}$ forms an orthonormal function sequence on
the real line.

In the following result, we state the connection between
$\\{\Psi_{n}\\}_{n\in\mathbb{Z}}$ and $\\{\varphi_{n}\\}_{n\in\mathbb{Z}}$ and
present an explicit formula for the fractional Laplacian of
$\\{\varphi_{n}\\}_{n\in\mathbb{Z}}$.

######  Lemma 2.6.

The Fourier transform of $\varphi_{n}(x)$ is

$\mathcal{F}[\varphi_{n}](\xi)=(-\mathrm{i})^{n}\Psi_{n}(\xi),$ (2.16)

and the fractional Laplacian of $\varphi_{n}(x)$ is

$(-\Delta)^{\alpha/2}\varphi_{n}(x)=\frac{\Gamma(\alpha+1)}{\sqrt{2\pi}}\left\\{\begin{array}[]{ll}{\displaystyle\frac{(-\mathrm{i})^{n}}{(\frac{1}{2}-\mathrm{i}x)^{\alpha+1}}~{}_{2}F_{1}\left(\begin{gathered}-n,\alpha+1\\\
1\end{gathered}\middle|\,\frac{1}{\frac{1}{2}-\mathrm{i}x}\right)},&n\geq
0,\\\\[14.0pt]
{\displaystyle\frac{-(-\mathrm{i})^{n}}{(\frac{1}{2}+\mathrm{i}x)^{\alpha+1}}~{}_{2}F_{1}\left(\begin{gathered}n+1,\alpha+1\\\
1\end{gathered}\middle|\,\frac{1}{\frac{1}{2}+\mathrm{i}x}\right)},&n<0,\end{array}\right.$
(2.17)

where ${}_{2}F_{1}(\cdot)$ is the Gauss hypergeometric function (see, e.g.,
[29, Chapter 15]).

###### Proof.

We only consider the case of $n\geq 0$ since the case of $n\leq-1$ can be
proved in a similar way. Taking the inverse Fourier transform of
$\Psi_{n}(x)$, we obtain

$\displaystyle\mathcal{F}^{-1}[\Psi_{n}](x)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}\Psi_{n}(\xi)e^{\mathrm{i}x\xi}\mathrm{d}\xi=\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}_{+}}e^{-\xi/2}L_{n}(\xi)e^{\mathrm{i}x\xi}\mathrm{d}\xi=\mathrm{i}^{n}\varphi_{n}(x),$

where we have used the formula [29, Equation (18.17.34)] in the last step. The
desired result (2.16) follows immediately by taking the Fourier transform of
the above equality. As for (2.17), using (1.2) and (2.16), we have

$\displaystyle(-\Delta)^{\alpha/2}\varphi_{n}(x)=(-\mathrm{i})^{n}\mathcal{F}^{-1}[|\xi|^{\alpha}\Psi_{n}(\xi)]=\frac{(-\mathrm{i})^{n}}{\sqrt{2\pi}}\int_{\mathbb{R}}|\xi|^{\alpha}\Psi_{n}(\xi)e^{\mathrm{i}\xi
x}\mathrm{d}\xi.$

The desired result (2.17) follows from applying [14, Equation (7.414.7)] to
the last equation. This ends the proof. ∎

######  Remark 2.7.

Let $\mu_{n}(x)=\mathrm{i}^{n}(\pi/2)^{1/2}\varphi_{n}(-x/2)$ for
$n\in\mathbb{Z}$. The fractional Laplacian of $\mu_{n}(x)$ was recently
derived in [7, Proposition 3.1] based on the techniques of complex analysis.
Here we provide an alternative approach for its derivation.

######  Remark 2.8.

From Lemma 2.6 we obtain immediately that
$\mathcal{F}[\varphi_{n}^{\mathrm{S}}](\xi)=(-\mathrm{i})^{n}\sqrt{\nu}\Psi_{n}(\nu\xi)$.

##  3 Malmquist-Takenaka spectral method

In this section we present a novel spectral discretization using MTFs in space
combined with time-stepping schemes for the temporal discretization for
solving the equation (1.3).

###  3.1 Spatial discretization

Let $\psi_{N}(x,t)=\sum_{k=-N}^{N-1}\zeta_{k}(t)\varphi_{k}(x)$ be the
spectral approximation to the solution of (1.3) and let $U(t)$ denote the
coefficient vector of $\psi_{N}(x,t)$, i.e.,
$U(t)=(\zeta_{-N}(t),\ldots,\zeta_{N-1}(t))^{T}$. From (2.5) it is easy to see
that $\psi_{N}(x,t)$ automatically satisfies the boundary condition
$\lim_{|x|\rightarrow{\infty}}\psi_{N}(x,t)=0$. Our spectral Galerkin method
is to find $\psi_{N}\in{\mathbb{V}_{N}}(\mathbb{R})$ such that

$\displaystyle\mathrm{i}(\partial_{t}{\psi_{N}},\phi)$
$\displaystyle=\gamma((-\Delta)^{\alpha/2}\psi_{N},\phi)+(\mathcal{T}\psi_{N},\phi),\quad\forall\phi\in{\mathbb{V}_{N}(\mathbb{R})},$
(3.1)

and $\psi_{N}(x,0)=\Pi_{N}\psi_{0}(x)$. Setting $\phi(x)=\varphi_{j}(x)$ in
(3.1) with $j=-N,\ldots,N-1$ and using the orthogonality property of MTFs, we
obtain that

$\displaystyle U^{\prime}(t)=-\mathrm{i}\gamma
AU(t)-\mathrm{i}\mathcal{N}(U,t),$ (3.2)

where $A\in\mathbb{C}^{2N\times{2N}}$, $\mathcal{N}(U,t)\in\mathbb{C}^{2N}$
are defined by

$\displaystyle
A=\big{(}(-\Delta)^{\alpha/2}\varphi_{k},\varphi_{j}\big{)}_{j,k=-N}^{N-1},\quad\mathcal{N}(U,t)=\big{(}\mathcal{T}\psi_{N},\varphi_{j}\big{)}_{j=-N}^{N-1}.$
(3.3)

Furthermore, recalling the Parseval’s equality, i.e.,

$\int_{\mathbb{R}}f(x)\overline{g(x)}\mathrm{d}x=\int_{\mathbb{R}}\mathcal{F}[f](\xi)\overline{\mathcal{F}[g](\xi)}\mathrm{d}\xi,\quad\forall
f,g\in L^{2}(\mathbb{R}),$

it follows that
$((-\Delta)^{\alpha/2}\varphi_{k},\varphi_{j})=(|\xi|^{\alpha}\mathcal{F}[\varphi_{k}],\mathcal{F}[\varphi_{j}])$.
Hence, the matrix $A$ in (3.3) can also be written as

$\displaystyle
A=\big{(}|\xi|^{\alpha}\mathcal{F}[\varphi_{k}],\mathcal{F}[\varphi_{j}]\big{)}_{j,k=-N}^{N-1}.$
(3.4)

Now, we consider the elements of the matrix $A$.

######  Lemma 3.1.

Let $A$ be the matrix defined in (3.3) or (3.4). Then $A$ is a Hermitian
matrix and can be written as a block two-by-two diagonal matrix of the form

$A=\left[\begin{array}[]{cc}PCP&\\\ &C\\\ \end{array}\right],$ (3.5)

where $P\in\mathbb{R}^{N\times N}$ is the permutation matrix which reverses
the order of a vector, i.e.,
$P(x_{1},\ldots,x_{N})^{T}=(x_{N},\ldots,x_{1})^{T}$, and
$C\in\mathbb{C}^{N\times N}$ is a Hermitian matrix whose elements are given by

$C_{j,k}=\mathrm{i}^{j-k}\sum_{\ell=0}^{\min\\{j,k\\}}\frac{(\alpha+1)_{\ell}(-\alpha)_{k-\ell}(-\alpha)_{j-\ell}}{\ell!(k-\ell)!(j-\ell)!},\quad
j,k=0,\ldots,N-1,$ (3.6)

and $(z)_{n}$ is the Pochhammer symbol defined by $(z)_{n}=(z)_{n-1}(z+n-1)$
for $n\geq 1$ and $(z)_{0}=1$. In the particular case of $\alpha=1$, then $C$
reduces to a tridiagonal and Hermitian matrix whose elements are given
explicitly by

$C_{j,k}=\begin{cases}j(-\mathrm{i}),\quad&k=j-1,\\\ (2j+1),\quad&k=j,\\\
(j+1)\mathrm{i},\quad&k=j+1.\end{cases}\quad j,k=0,\ldots,N-1.$ (3.7)

###### Proof.

Combining (3.4) with Lemma 2.6 we have

$\displaystyle
A_{j,k}=\left(|\xi|^{\alpha}\mathcal{F}[\varphi_{k}],\mathcal{F}[\varphi_{j}]\right)$
$\displaystyle=\mathrm{i}^{j-k}\int_{\mathbb{R}}|\xi|^{\alpha}\Psi_{k}(\xi)\Psi_{j}(\xi)\mathrm{d}\xi$
$\displaystyle=\mathrm{i}^{j-k}\begin{cases}{\displaystyle\int_{\mathbb{R}_{+}}\xi^{\alpha}e^{-\xi}L_{k}(\xi)L_{j}(\xi)\mathrm{d}\xi},\quad&k,j\geq
0,\\\\[12.0pt]
{\displaystyle\int_{\mathbb{R}_{+}}\xi^{\alpha}e^{-\xi}L_{-k-1}(\xi)L_{-j-1}(\xi)\mathrm{d}\xi},\quad&k,j<0.\end{cases}$
(3.8)

It is easily seen that $A$ is a block two-by-two diagonal matrix and the first
block can be derived by reversing the order of rows and columns of the second
block. Therefore, we restrict our attention to the second block, which is
denoted by $C$. Recalling the connection formula of Laguerre polynomials (see
[29, Equation (18.18.18)]), we have

$L_{n}(\xi)=\sum_{j=0}^{n}\frac{(-\alpha)_{n-j}}{(n-j)!}L^{(\alpha)}_{j}(\xi).$
(3.9)

The desired result (3.6) follows by combining (3.1), (3.9) and the
orthogonality property of the Laguerre polynomials
$\\{L_{k}^{(\alpha)}(\xi)\\}$. In the particular case of $\alpha=1$, recalling
the three term recurrence relation of Laguerre polynomials, i.e.,

$(k+1)L_{k+1}(\xi)=(2k+1-\xi)L_{k}(\xi)-kL_{k-1}(\xi),$ (3.10)

The desired result (3.7) follows immediately by combining (3.1) and the
orthogonality of Laguerre polynomials $\\{L_{k}(\xi)\\}$. This ends the proof.
∎

######  Remark 3.2.

We define

$\beta_{\ell,n}^{(\alpha)}=\frac{(-\alpha)_{n-\ell}}{(n-\ell)!}\sqrt{\frac{(\alpha+1)_{\ell}}{\ell!}},\quad\ell=0,\ldots,n,$
(3.11)

where $n=0,1,\ldots$. Note that $\beta_{\ell,n}^{(\alpha)}=0$ for $\alpha=1$
and $n-\ell>1$. Combining (3.1), (3.9) and the orthogonality property of the
Laguerre polynomials $\\{L_{k}^{(\alpha)}(\xi)\\}$, we obtain that

$\displaystyle C_{j,k}$
$\displaystyle=\mathrm{i}^{j-k}\int_{\mathbb{R}_{+}}\xi^{\alpha}e^{-\xi}L_{k}(\xi)L_{j}(\xi)\mathrm{d}\xi=\mathrm{i}^{j-k}\sum_{\ell=0}^{\min\\{j,k\\}}\beta_{\ell,j}^{(\alpha)}\beta_{\ell,k}^{(\alpha)},$
(3.12)

and thus

$\displaystyle
C=\left(\begin{array}[]{ccc}\mathrm{i}^{0}&\cdots&\mathrm{i}^{1-N}\\\
\vdots&\ddots&\vdots\\\ \mathrm{i}^{N-1}&\cdots&\mathrm{i}^{0}\\\
\end{array}\right)\circ$
$\displaystyle\left[\left(\begin{array}[]{ccc}\beta_{0,0}^{(\alpha)}&&\\\
\vdots&\ddots&\\\
\beta_{0,N-1}^{(\alpha)}&\cdots&\beta_{N-1,N-1}^{(\alpha)}\\\
\end{array}\right)\right.$ (3.19)
$\displaystyle~{}~{}~{}~{}~{}\times\left.\left(\begin{array}[]{ccc}\beta_{0,0}^{(\alpha)}&&\\\
\vdots&\ddots&\\\
\beta_{0,N-1}^{(\alpha)}&\cdots&\beta_{N-1,N-1}^{(\alpha)}\\\
\end{array}\right)^{T}\right],$ (3.23)

where $\circ$ denotes the Hadamard product.

###  3.2 Temporal discretization

In this section we consider time-stepping schemes for the temporal
discretization of (3.2). We divide our discussion into two cases according to
$\mathcal{T}$ is a linear or nonlinear operator.

####  3.2.1 The linear case

We restrict our attention to the case $\mathcal{T}\psi=V(x)\psi(x,t)$, where
$V(x)$ is a smooth potential. From the definition of $\mathcal{N}(U,t)$ in
(3.3) we obtain $\mathcal{N}(U,t)=MU(t)$, where $M\in\mathbb{C}^{2N\times 2N}$
is defined by $M_{j,k}=(\mathcal{T}\varphi_{k},\varphi_{j})$ with
$j,k=-N,\ldots,N-1$. For $k\in\mathbb{Z}$, we define

$\mu_{k}=\frac{\mathrm{i}^{-k}}{2\pi}\int_{-\pi}^{\pi}V\left(\frac{1}{2}\tan{\frac{\theta}{2}}\right)e^{-\mathrm{i}k\theta}\mathrm{d}\theta.$
(3.24)

It can be verified by direct calculation that $M_{j,k}=\mu_{j-k}$, and thus

$M=\left(\begin{array}[]{cccc}\mu_{0}&\mu_{-1}&\cdots&\mu_{1-2N}\\\
\mu_{1}&\mu_{0}&\cdots&\mu_{2-2N}\\\ \vdots&\vdots&\ddots&\vdots\\\
\mu_{2N-1}&\mu_{2N-2}&\cdots&\mu_{0}\end{array}\right).$ (3.25)

Clearly, we see that $M$ is a Toeplitz and Hermitian matrix. Furthermore,
observe that the integrand on the right-hand of (3.24) is periodic in $\theta$
with period $2\pi$, the elements of $M$ (i.e.,
$\\{\mu_{k}\\}_{k=1-2N}^{2N-1}$) can be computed rapidly with the FFT in
$\mathcal{O}(N\log N)$ operations.

Let us now turn to the numerical solution of (3.2). In view of
$\mathcal{N}(U,t)=MU(t)$, we obtain the following ODE system

$U^{\prime}(t)=-\mathrm{i}(\gamma{A}+M)U(t).$ (3.26)

The exact solution of (3.26) is $U(t)=\exp(-\mathrm{i}(\gamma{A}+M)t)U(0)$ and
$U(0)$ can be computed from the Malmquist-Takenaka coefficients of
$\psi_{0}(x)$ by the FFT. Note that the exact solution requires the
computation of the matrix exponential $\exp(-\mathrm{i}(\gamma{A}+M)t)$, which
is generally not a good idea to compute it directly. Indeed, as we will see in
Figures 3 and 4, the decay rates of the magnitudes of the elements of $A$ and
$M$ (or, equivalently, $C$ and $M$) are quite different. To avoid this issue,
we consider the use of splitting method to solve (3.26). Specifically, let
$t_{k}=k\tau$ denote the time grid points, where $\tau>0$ is the time step
size, and let $U_{k}$ denote the approximation to the exact value $U(t_{k})$
and $U_{0}=U(0)$. From $t_{n-1}$ to $t_{n}$, the splitting method reads

$U_{n}=\underbrace{\left[\prod_{j=1}^{m}\exp(-\mathrm{i}a_{j}\gamma\tau
A)\exp(-\mathrm{i}b_{j}\tau M)\right]}_{:=~{}S(\tau)}U_{n-1},\quad n\geq 1,$
(3.27)

where $a_{j}$ and $b_{j}$ are some suitably chosen coefficients to ensure that
the method achieves some order $p$, i.e.,
$S(\tau)=\exp(-\mathrm{i}(\gamma{A}+M)\tau)+\mathcal{O}(\tau^{p+1})$. We list
below three symmetric splitting methods which achieve orders two, four and
six, respectively:

  * •

SM1:

$a_{1}=a_{2}=\frac{1}{2},~{}~{}~{}b_{1}=1,~{}b_{2}=0.$ (3.28)

This method is known as the Strang splitting [33].

  * •

SM2:

$a_{1}=a_{4}=\frac{\kappa_{1}}{2},~{}a_{2}=a_{3}=\frac{\kappa_{0}+\kappa_{1}}{2},~{}~{}~{}b_{1}=b_{3}=\kappa_{1},~{}b_{2}=\kappa_{0},~{}b_{4}=0,$
(3.29)

where $\kappa_{0}=-2^{1/3}/(2-2^{1/3})$ and $\kappa_{1}=1/(2-2^{1/3})$.

  * •

SM3:

$\displaystyle a_{1}$
$\displaystyle=a_{8}=\frac{w_{3}}{2},~{}a_{2}=a_{7}=\frac{w_{2}+w_{3}}{2},~{}a_{3}=a_{6}=\frac{w_{1}+w_{2}}{2},~{}a_{4}=a_{5}=\frac{w_{0}+w_{1}}{2},$
$\displaystyle b_{1}$
$\displaystyle=b_{7}=w_{3},~{}b_{2}=b_{6}=w_{2},~{}b_{3}=b_{5}=w_{1},~{}b_{4}=w_{0},~{}b_{8}=0,$
(3.30)

where $w_{1}=-1.17767998417887$, $w_{2}=0.235573213359$,
$w_{3}=0.784513610477$ and $w_{0}=1-2(w_{1}+w_{2}+w_{3})$.

For the derivation of the coefficients of these splitting methods, we refer to
[37] for more details.

Figure 3: The logarithm of the magnitude of the elements of $C$ for
$\alpha=0.6$ (left) and $\alpha=1.4$ (right).

######  Remark 3.3.

Since both $A$ and $M$ are Hermitian matrices, those matrix exponentials in
(3.27) are all unitary and thus unitary evolution and unconditional stability
of the method are guaranteed.

When implementing splitting methods, it is necessary to evaluate matrix
exponential of the forms $\exp(-\mathrm{i}\lambda{A})$ and
$\exp(-\mathrm{i}\lambda M)$, where $\lambda>0$ is a constant. Moreover, from
(3.5) we obtain that

$\exp\left(-\mathrm{i}\lambda{A}\right)=\left[\begin{array}[]{cc}P\exp\left(-\mathrm{i}\lambda{C}\right)P&\\\
&\exp\left(-\mathrm{i}\lambda{C}\right)\\\ \end{array}\right],$ (3.31)

and thus the computation of $\exp(-\mathrm{i}\lambda{A})$ can be reduced to
the computation of $\exp(-\mathrm{i}\lambda{C})$. We now consider the
structure of the matrices $C$ and $M$. In the case $\alpha\neq 1$, from (3.11)
we deduce that $\beta_{\ell,n}^{(\alpha)}=\mathcal{O}(n^{-\alpha-1})$ for
fixed $\ell$ and $n\gg 1$ and therefore the magnitude of the elements of $C$
decays algebraically away from the diagonal; see Figure 3. In the case
$\alpha=1$, from Lemma 3.1 we know that $C$ is a Hermitian and tridiagonal
matrix. Moreover, from (3.25) we know that $M$ is also a Hermitian matrix and
from (3.24) it is easily seen that ${\mu_{k}}$ decays rapidly whenever
$V\in{H^{s}(\mathbb{R})}$ for sufficiently large $s\in\mathbb{N}$, and thus we
can expect that $M$ is near a banded matrix; see Figure 4. Now we turn to the
computation of $\exp(-\mathrm{i}\lambda{C})$ and
$\exp(-\mathrm{i}\lambda{M})$. For this, we use the recently developed
algorithm in [3], which is based on Chebyshev approximation of complex
exponentials. We only consider the computation of
$\exp(-\mathrm{i}\lambda{C})$ since the computation of
$\exp(-\mathrm{i}\lambda{M})$ is similar. Specifically, if the eigenvalues of
$C$ are contained in the interval $[\zeta,\eta]\subset\mathbb{R}$ and
$\lambda(\eta-\zeta)/2\leq 2.212$, the algorithm in [3] reads

$\exp(-\mathrm{i}\lambda{C})\approx\exp\left(-\mathrm{i}\lambda\frac{\eta+\zeta}{2}\right)\left[c_{0}I_{N}+2\sum_{k=1}^{m}c_{k}T_{k}\left(\frac{2}{\eta-\zeta}\left(C-\frac{\zeta+\eta}{2}I_{N}\right)\right)\right],$
(3.32)

where $I_{N}$ is the identity matrix of order $N$ and
$c_{k}=(-\mathrm{i})^{k}J_{k}(\lambda(\eta-\zeta)/2)$ and $J_{k}(x)$ is the
Bessel function of the first kind of order $k$. When choosing $m=18$, the
approximation error in (3.32) will be less than the machine precision (i.e.,
$2^{-53}\approx 1.11\times 10^{-16}$) and the calculation of the right-hand
side of (3.32) can be achieved with only five matrix-matrix products.
Otherwise, if $\lambda(\eta-\zeta)/2>2.212$, then the scaling and squaring
technique, i.e.,
$\exp(-\mathrm{i}\lambda{C})=(\exp(-\mathrm{i}\lambda{C}/2^{s}))^{2^{s}}$ for
some $s\in\mathbb{N}$, should be used such that the exponential
$\exp(-\mathrm{i}\lambda{C}/2^{s})$ can be evaluated by (3.32). Note that the
values of $\zeta$ and $\eta$ have to be specified before embarking on the
algorithm. If they cannot be specified in advance, then the algorithm will
simply take $\eta=-\zeta=\|C\|_{1}$.

Figure 4: The logarithm of the magnitude of the elements of $M$ for
$V(x)=1/(1+x^{2})$ (left) and $V(x)=\exp(-x^{2})$ (right).

######  Remark 3.4.

It is possible to improve the efficiency of the algorithm (3.32) by taking the
structure of $C$ and $M$ into account. For example, we can simply set the
elements of the matrices involved to zero whenever their magnitude is less
than the machine precision.

###  3.3 The nonlinear case

In the case where $\mathcal{T}$ is a nonlinear operator, using the variation-
of-constant formula to (3.2), we have

$U(t_{n+1})=e^{-\mathrm{i}\gamma\tau{A}}U(t_{n})-\mathrm{i}\tau\int_{0}^{1}e^{-\mathrm{i}\gamma\tau{A}(1-y)}\mathcal{N}(U(t_{n}+y\tau),t_{n}+y\tau)\mathrm{d}y.$
(3.33)

To approximate (3.33), the fourth-order exponential time differencing Runge-
Kutta (ETDRK4) method and its various modifications are preferable (see, e.g.,
[4, 9, 20, 24]). Here we utilize the Krogstad-P22 scheme developed in [4].
More specifically, let $L=\mathrm{i}\gamma{A}$ and let
$R_{2,2}(\tau{L})=(12I-6\tau{L}+\tau^{2}{L}^{2})(12{I}+6\tau{L}+\tau^{2}{L}^{2})^{-1}$,
where we have omitted the subscript $2N$ on the identity matrix $I$ for
notational simplicity. Moreover, we define

$\displaystyle P_{1}(\tau{L})$
$\displaystyle=12\tau(12{I}+6\tau{L}+\tau^{2}{L}^{2})^{-1},$ $\displaystyle
P_{2}(\tau{L})$
$\displaystyle=\tau(6{I}+\tau{L})(12{I}+6\tau{L}+\tau^{2}{L}^{2})^{-1},$
$\displaystyle P_{3}(\tau{L})$
$\displaystyle=2\tau(4{I}+\tau{L})(12{I}+6\tau{L}+\tau^{2}{L}^{2})^{-1}.$

Then, the Krogstad-P22 scheme reads

$\displaystyle U_{n+1}$
$\displaystyle=R_{2,2}(\tau{L})U_{n}-\mathrm{i}{P}_{1}(\tau{L})\mathcal{N}(U_{n},t_{n})-\mathrm{i}{P}_{2}(\tau{L})\bigg{[}-3\mathcal{N}(U_{n},t_{n})+2\mathcal{N}\left(a_{n},t_{n}+\frac{\tau}{2}\right)$
$\displaystyle~{}~{}~{}+2\mathcal{N}\left(b_{n},t_{n}+\frac{\tau}{2}\right)-\mathcal{N}(c_{n},t_{n}+\tau)\bigg{]}-\mathrm{i}P_{3}(\tau{L})\bigg{[}\mathcal{N}(U_{n},t_{n})$
(3.34)
$\displaystyle~{}~{}~{}-\mathcal{N}\left(a_{n},t_{n}+\frac{\tau}{2}\right)-\mathcal{N}\left(b_{n},t_{n}+\frac{\tau}{2}\right)+\mathcal{N}(c_{n},t_{n}+\tau)\bigg{]},$

where

$\displaystyle a_{n}$
$\displaystyle=\tilde{R}_{2,2}(\tau{L})U_{n}-\mathrm{i}\tilde{P}_{1}(\tau{L})\mathcal{N}(U_{n},t_{n}),$
$\displaystyle b_{n}$
$\displaystyle=\tilde{R}_{2,2}(\tau{L})U_{n}-\mathrm{i}\tilde{P}_{1}(\tau{L})\mathcal{N}(U_{n},t_{n})-\mathrm{i}\tilde{P}_{2}(\tau{L})\left[\mathcal{N}\left(a_{n},t_{n}+\frac{\tau}{2}\right)-\mathcal{N}(U_{n},t_{n})\right],$
$\displaystyle c_{n}$
$\displaystyle=R_{2,2}(\tau{L})U_{n}-\mathrm{i}P_{1}(\tau{L})\mathcal{N}(U_{n},t_{n})-2\mathrm{i}P_{2}(\tau{L})\left[\mathcal{N}\left(b_{n},t_{n}+\frac{\tau}{2}\right)-\mathcal{N}(U_{n},t_{n})\right],$

and

$\displaystyle\tilde{R}_{2,2}(\tau{L})$
$\displaystyle=(48{I}-12\tau{L}+\tau^{2}{L}^{2})(48{I}+12\tau{L}+\tau^{2}{L}^{2})^{-1},$
$\displaystyle\tilde{P}_{1}(\tau{L})$
$\displaystyle=24\tau(48{I}+12\tau{L}+\tau^{2}{L}^{2})^{-1},$
$\displaystyle\tilde{P}_{2}(\tau{L})$
$\displaystyle=2\tau(12{I}+\tau{L})(48{I}+12\tau{L}+\tau^{2}{L}^{2})^{-1}.$

######  Remark 3.5.

The main drawback of the ETDRK4 and ETDRK4-B schemes [9, 24] is the
computation of the following matrix functions

$\varphi_{0}(\tau{L})=\exp(-\tau{L}),\quad\varphi_{k}(\tau{L})=(-\tau
L)^{-k}\left(\varphi_{0}(\tau{L})-\sum_{j=0}^{k-1}\frac{(-\tau{L})^{j}}{j!}\right),\quad
k=1,2,3,$

which suffers from cancellation errors, especially when the eigenvalues of $L$
are very close to zero. Kassam and Trefethen in [20] proposed to evaluate
these functions by using complex contour integrals. However, the choice of the
contour is problem dependent. The Krogstad-P22 scheme [4] is developed by
using Padé approximations to the above functions, which avoids direct
computation of the higher powers of matrix inverse. Moreover, as observed in
[4], the factors ${L}^{-1}$ and ${L}^{-3}$ that appear in ETDRK4 and ETDRK4-B
schemes cancel out in the Krogstad-P22 scheme.

The stability of the Krogstad-P22 scheme can be analyzed by using a similar
argument as for the ETDRK4 scheme in [9]. For the nonlinear and autonomous ODE
of the form $U_{t}=-\lambda{U}-\mathrm{i}\mathcal{N}(U)$. Suppose that there
exists a fixed point $U_{0}$ such that
$\lambda{U_{0}}+\mathrm{i}\mathcal{N}(U_{0})=0$. Linearizing about this fixed
point yields

$U_{t}=-\lambda{U}-\mathrm{i}cU,$ (3.35)

where $U$ is now the perturbation to $U_{0}$ and
$c=\mathcal{N}^{\prime}(U_{0})$. Applying the Krogstad-P22 scheme (3.3) to the
linearized equation (3.35) and setting $r=U_{n+1}/U_{n}$, $x=-\mathrm{i}c\tau$
and $y=-\lambda\tau$, we then obtain

$r(x,y)=c_{0}+c_{1}x+c_{2}x^{2}+c_{3}x^{3}+c_{4}x^{4},$ (3.36)

where

$\displaystyle c_{0}$ $\displaystyle=\frac{12+6y+y^{2}}{12-6y+y^{2}},$
$\displaystyle c_{1}$
$\displaystyle=\frac{144y^{3}-432y^{2}-1728y+6912}{{(12-6y+y^{2})}^{2}{(48-12y+y^{2})}},$
$\displaystyle c_{2}$
$\displaystyle=\frac{36y^{5}-648y^{4}+4032y^{3}+3456y^{2}-82944y+165888}{{(12-6y+y^{2})}^{2}{(48-12y+y^{2})}^{2}},$
$\displaystyle c_{3}$
$\displaystyle=\frac{-48y^{4}+768y^{3}+576y^{2}-27648y+55296}{{(12-6y+y^{2})}^{2}{(48-12y+y^{2})}^{2}},$
$\displaystyle c_{4}$
$\displaystyle=\frac{-96y^{3}+1920y^{2}-10368y+13824}{{(12-6y+y^{2})}^{2}{(48-12y+y^{2})}^{2}}.$

Finally, the stability region of the Krogstad-P22 scheme can be obtained by
requiring $|r(x,y)|\leq 1$. Since all eigenvalues of ${L}$ are pure imaginary,
it is therefore enough to consider the case where $\lambda$ is a pure
imaginary number. In Figure 5 we plot the stability regions for several pure
imaginary numbers of $y$ for (3.2). We observe that each of the stability
region includes an interval of the imaginary axis and its length increases as
$|y|$ increases. This gives an indication of stability of the Krogstad-P22
scheme.

Figure 5: Stability regions of the Krogstad-P22 scheme for several values of
$y$.

##  4 Numerical examples

In this section, we present two examples to show the performance of the
proposed method.

Example 1. Consider the following linear fractional Schrödinger equation

$\mathrm{i}\partial_{t}\psi(x,t)=\gamma(-\Delta)^{\alpha/2}\psi(x,t)+V(x)\psi(x,t),$
(4.1)

where we take $\gamma=1/2$ and $V(x)=1/(1+x^{2})$. We first compare the
performance of spatial discretizations using MTFs and MCFs, respectively. To
avoid the influence of the error due to temporal discretization, we consider a
pure version of spectral methods, that is, we evaluate $U(t)$ by the exact
formula $U(t)=\exp(-\mathrm{i}(\gamma{A}+M)t)U(0)$ and compute the involved
matrix exponential by the expm function in Matlab. Since the exact solution of
(4.1) is not known, we define a reference solution which is computed by the
MTF spectral Galerkin method with $N=500$ (note that the number of terms of
this spectral method is $2N$). Moreover, we take the scaling parameter $\nu=4$
for both methods.

In the first row of Figure 6 we plot the maximum error of both spectral
methods at time $t=1$ for the initial data $\psi_{0}(x)=\mathrm{sech}(x)$. We
see that MTF spectral discretization converges much faster than its MCF
counterpart in the case of $\alpha=1$ and both spectral discretizations
converge almost at the same rate otherwise. Indeed, in the case of $\alpha=1$,
MTF spectral discretization scheme converges at exponential rates, while MCF
spectral discretization scheme converges only at algebraic rates. In the
second row of Figure 6 we plot the maximum error at time $t=1$ for the initial
data $\psi_{0}(x)=(\mathrm{i}x+10)/(x^{2}+4)$. We observe that MTF spectral
discretization converges always faster than its MCF counterpart. Moreover,
similar to the previous case, our MTF spectral discretization converges at
exponential rates in the particular case of $\alpha=1$ while its MCF
counterpart converges only at algebraic rates.

Figure 6: Maximum errors of MTF and MCF spectral discretizations for
$\alpha=1$ (left), $\alpha=0.6,1.4$ (right). Here
$\psi_{0}(x)=\mathrm{sech}(x)$ (top row) and
$\psi_{0}(x)=(\mathrm{i}x+10)/(x^{2}+4)$ (bottom row).

Figure 7: Temporal orders of splitting methods SM1 (left), SM2 (middle) and
SM3 (right) coupled with MTF spectral Galerkin method for (4.1).

Next, we consider the temporal order of convergence of the MTF spectral
Galerkin method coupled with the splitting methods. In this case, the
reference solution is computed with the MTF spectral Galerkin method coupled
with the splitting method (• ‣ 3.2.1) with $\tau=2^{-11}$ and $N=500$. We
consider the maximum error of the numerical solution at time $t=1$ for the
initial data $\psi_{0}(x)=1/(1+x+x^{2})$. In Figure 7 we plot the maximum
error of MTF spectral Galerkin method coupled with the splitting methods
(3.28), (3.29) and (• ‣ 3.2.1) as a function of the time step size $\tau$.
Clearly, we see that the classical order of each splitting scheme is retained.

Finally, we consider the evolution of $|\psi(x,t)|$ in (4.1) with a double-
barrier potential of the form $V(x)=100(\exp(-(x-10)^{2})+\exp(-(x+10)^{2}))$.
The initial data is taken as $\psi_{0}(x)=\exp(-x^{2})\exp(-\mathrm{i}\kappa
x)$, where $\kappa\geq 0$. This equation was studied in [16] in the context of
beam propagation and the variables $x$ and $t$ denote the normalized
transverse and longitudinal coordinates, respectively. In our simulations, we
use the MTF spectral Galerkin method with $N=500$ and the scaling parameter
$\nu=4$ in space, combined with the splitting method SM3 with $\tau=0.001$ in
time. The top row of Figure 8 shows the evolution of $|\psi(x,t)|$ for
$\kappa=0$. We see that $|\psi(x,t)|$ splits and is diffraction-free in the
case of $\alpha=1$ and splits and diffracts in the case of $\alpha\in(1,2)$.
Moreover, we also see that the diffraction of the beam becomes stronger as
$\alpha$ increases. The bottom row of Figure 8 shows the evolution of
$|\psi(x,t)|$ for $\kappa=10$, which indicates that the initial data can be
viewed as an oblique Gaussian beam. We see that $|\psi(x,t)|$ exhibits
diffraction-free propagation in the case of $\alpha=1$ and exhibits
diffraction in the case of $\alpha\in(1,2)$. Moreover, similar to the previous
case, the diffraction of the beam becomes stronger as $\alpha$ increases. Our
results are consistent with the observations in [16].

Figure 8: Time evolution of $|\psi(x,t)|$ for the initial data
$\psi_{0}(x)=\exp(-x^{2})\exp(-\mathrm{i}\kappa x)$. Here $\kappa=0$ (top row)
and $\kappa=10$ (bottom row).

Example 2. Let us consider the following focusing nonlinear FSE (see, e.g.,
[7, 11, 23])

$\mathrm{i}\partial_{t}\psi(x,t)=\gamma(-\Delta)^{\alpha/2}\psi(x,t)-|\psi(x,t)|^{2}\psi(x,t),$
(4.2)

where $\gamma=1/2$. It is known that the solution satisfies the mass
conservation, i.e.,

$M(t)=\int_{\mathbb{R}}|\psi(x,t)|^{2}\mathrm{d}x=M(0).$ (4.3)

We discretize the equation (4.2) using the MTF spectral Galerkin method in
space coupled with the Krogstad-P22 scheme (3.3) in time. The nonlinear terms
in (3.3) are computed as follows: We only consider the term
$\mathcal{N}(U_{n},t_{n})$ since the other terms
$\mathcal{N}(a_{n},t_{n}+\tau/2)$, $\mathcal{N}(b_{n},t_{n}+\tau/2)$ and
$\mathcal{N}(c_{n},t_{n}+\tau)$ can be computed similarly. Firstly, let
$\\{\varrho_{k}\\}_{k\in\mathbb{Z}}$ be the sequence defined by

$\varrho_{k}=\frac{\mathrm{i}^{-k}}{\pi^{2}}\int_{-\pi}^{\pi}\bigg{|}\cos\left(\frac{\theta}{2}\right)\sum_{j=-N}^{N-1}\zeta_{j}(t_{n})\mathrm{i}^{j}e^{\mathrm{i}j\theta}\bigg{|}^{2}e^{-\mathrm{i}k\theta}\mathrm{d}\theta.$
(4.4)

From (3.3) one can verify that
$\mathcal{N}(U_{n},t_{n})=\mathcal{B}(U_{n})U_{n}$, where
$\mathcal{B}(U_{n})\in\mathbb{C}^{2N\times 2N}$ is defined by

$\mathcal{B}(U_{n})=\left(\begin{array}[]{cccc}\varrho_{0}&\varrho_{-1}&\cdots&\varrho_{1-2N}\\\
\varrho_{1}&\varrho_{0}&\cdots&\varrho_{2-2N}\\\
\vdots&\vdots&\ddots&\vdots\\\
\varrho_{2N-1}&\varrho_{2N-2}&\cdots&\varrho_{0}\end{array}\right).$ (4.5)

It is clear that $\mathcal{B}(U_{n})$ is a Toeplitz and Hermitian matrix.
Furthermore, note that the integrand on the right-hand of (4.4) is periodic in
$\theta$ with period $2\pi$, hence the elements of $\mathcal{B}(U_{n})$ (i.e.,
$\\{\varrho_{k}\\}_{k=1-2N}^{2N-1}$) can be computed rapidly with the FFT. On
the other hand, since $\mathcal{B}(U_{n})$ is a Toeplitz matrix,
$\mathcal{N}(U_{n},t_{n})=\mathcal{B}(U_{n})U_{n}$ can also be computed by the
FFT [13].

We first consider the initial data $\psi_{0}(x)=\mathrm{sech}(x)$. The top row
of Figure 9 shows the time evolution of the solution $|\psi(x,t)|^{2}$ for the
initial data $\psi_{0}(x)=\mathrm{sech}(x)$. Note that the focusing nonlinear
equation (4.2) will display the finite time blow-up phenomenon whenever
$\alpha\in(1/2,1]$ (see [23]). In our simulations, we take
$N=300,\nu=10,\tau=0.001$. Clearly, the finite time blow-up phenomenon is
observed for $\alpha=0.8,1$. We then consider the initial data
$\psi_{0}(x)=\mathrm{exp}(-x^{2})$. In this case, the mass can be calculated
as

$M(0)=\int_{\mathbb{R}}\mathrm{exp}(-2x^{2})\mathrm{d}x=\sqrt{\frac{\pi}{2}}.$
(4.6)

The bottom row of Figure 9 shows the evolution of the solution
$|\psi(x,t)|^{2}$ for different values of $\alpha$ and we take
$N=300,~{}\nu=4,~{}\tau=0.001$ in our simulations. Again, the finite time
blow-up phenomenon is observed for $\alpha=0.6,0.8$. Now, we turn our
attention to the property of mass conservation. Using the change of variable
$x=\tan(\theta/2)/2$, the mass error can be written as

$\displaystyle\mathrm{err}(t)=\left|\int_{\mathbb{R}}|\psi_{N}(x,t)|^{2}\mathrm{d}x-M(0)\right|=\bigg{|}\frac{1}{2\pi}\int_{-\pi}^{\pi}\bigg{|}\sum_{k=-N}^{N-1}\zeta_{k}(t)\mathrm{i}^{k}e^{\mathrm{i}{k}\theta}\bigg{|}^{2}\mathrm{d}\theta-\sqrt{\frac{\pi}{2}}\bigg{|},$

and the integral in the last equation can be evaluated by the inverse FFT. In
Figure 10 we plot the mass error of the MTF spectral method coupled with the
Krogstad-P22 scheme for $\alpha=1.99$ and we choose $N=150,\nu=4,\tau=0.001$.
We see that the mass errors are of the order $10^{-13}$ for $t\in[0,50]$,
which indicates that our method is mass conserved. Moreover, compared with the
mass error in [7, Figure 12], which achieves the order $10^{-7}$ with the same
number of spectral discretizations, our method is clearly much better.

Figure 9: Time evolution of $|\psi(x,t)|^{2}$ in (4.2) with the initial data
$\psi_{0}(x)=\mathrm{sech}(x)$ (top row) and
$\psi_{0}(x)=\mathrm{exp}(-x^{2})$ (bottom row).  Figure 10: Mass error of the
MTF spectral Galerkin method coupled with the Krogstad-P22 scheme to (4.2) for
$\alpha=1.99$.

Finally, we illustrate the accuracy of the MTF spectral discretization in
space and the temporal order of convergence of the Krogstad-P22 scheme (3.3)
in time. The reference solution is computed by using the MTF spectral Galerkin
method coupled with (3.3) with $N=300$ and $\tau=10^{-4}$. In Figure 11 we
plot the maximum error of MTF spectral method with the scaling parameter
$\nu=4$ coupled with the Krogstad-P22 scheme (3.3) at time $t=1$. We can see
that, similar to the linear case, the MTF spectral discretization converges at
an exponential rate in the case $\alpha=1$ and at an algebraic rate in the
case $\alpha\neq{1}$. In the left graph of Figure 12, we plot the maximum
error of the MTF spectral method with $\nu=4$ coupled with (3.3) at time $t=1$
as a function of the time step size $\tau$. We see that the temporal order of
convergence is four for all choices of $\alpha$. In the right graph of Figure
12, we plot the asymptotic behavior of the computed solutions at time $t=1$
for three values of $\alpha$. We see that all computed solutions decay at the
rate $\mathcal{O}(|x|^{-\alpha-1})$, which are in agreement with the numerical
observation in [34].

Figure 11: Maximum errors of the MTF spectral method coupled with the
Krogstad-P22 scheme at time $t=1$ as a function of $N$.

Figure 12: Left: Maximum errors of MTF spectral method coupled with the
Krogstad-P22 scheme (3.3) at time $t=1$ as a function of $\tau$. Right:
Asymptotic behavior of the solution $\psi(x,t=1)$ for three values of
$\alpha$.

##  5 Conclusions

In this paper, we have studied the numerical solution of the FSE on the real
line. We proposed a new spectral discretization using MTFs in space combined
with some time-stepping methods for the discretization in time. This new
spectral discretization can achieve exponential convergence in space in the
case of $\alpha=1$, regardless of the underlying FSE is linear or nonlinear,
and exhibits a comparable or even better performance than state-of-the-art
spectral discretization schemes in other cases. We conclude that spectral
methods using MTFs are competitive for solving PDEs whose solution has slow
decay behavior at infinity.

Before closing this paper, we list several problems for future research:

  * •

Our study in this work is restricted to the one-dimensional FSE problem.
However, an extension of the current work to the high dimensional FSE problem
is straightforward. An issue arised in the process of extension is a
multivariate counterpart of the integral in (3.1), which might be difficult to
evaluate due to the singular and nonseparable factor $|\xi|^{\alpha}$.

  * •

Both MTFs and MCFs are orthogonal systems on the real line and spectral
approximations using them can be achieved rapidly by using the FFT. An
interesting problem is to compare their approximation powers for functions
with exponential or algebraic decay behavior at infinity.

  * •

Conservation laws, such as mass and energy, are important in quantum
mechanics. In the case of the nonlinear FSE, Figure 10 implies that the MTF
spectral Galerkin method coupled with the Krogstad-P22 scheme is mass
conserved. However, a rigorous analysis of this observation is still lacking.

We will address these issues in the future.

## Acknowledgements

This work was supported by the National Natural Science Foundation of China
under grant 11671160. The authors are grateful to Prof. Chengming Huang for
his valuable suggestions in improving this paper.

## References

  * [1] G. Acosta, F. M. Bersetche and J. P. Borthagaray, A short FE implementation for a 2d homogeneous Dirichlet problem of a fractional Laplacian, Comput. Math. Appl., 74(4):784-816, 2017.
  * [2] G. Acosta, J. P. Borthagaray, O. Bruno and M. Maas, Regularity theory and high order numerical methods for the (1D)-fractional Laplacian, Math. Comp., 87(312):1821-1857, 2018.
  * [3] P. Bader, S. Blanes, F. Casas and M. Seydaoğlu, An efficient algorithm to compute the exponential of skew-Hermitian matrices for the time integration of the Schrödinger equation, Math. Comput. Simulation, 194:383-400, 2022.
  * [4] H. P. Bhatt and A. Q. M. Khaliq, Fourth-order compact schemes for the numerical simulation of coupled Burgers’ equation, Comput. Phys. Commun., 200:117-138, 2016.
  * [5] A. Bonito, W.-Y. Lei and J. E. Pasciak, Numerical approximation of the integral fractional Laplacian, Numer. Math., 142(2):235-278, 2019.
  * [6] J. Cayama, C. M. Cuesta and F. de la Hoz, A pseudospectral method for the one-dimensional fractional Laplacian on $\mathbb{R}$, Appl. Math. Comput., 389:125577, 2021.
  * [7] J. Cayama, C. M. Cuesta and F. de la Hoz, Numerical approximation of the fractional Laplacian on $\mathbb{R}$ using orthogonal families, Appl. Numer. Math., 158:164-193, 2020.
  * [8] C. Christov, A complete orthonormal system of functions in $L^{2}(-\infty,\infty)$ space, SIAM J. Appl. Math., 42(6):1337-1344, 1982.
  * [9] S. M. Cox and P. C. Matthews, Exponential time differencing for stiff systems, J. Comput. Phys., 176(2):430-455, 2002.
  * [10] Q. Du, Nonlocal Modeling, Analysis, and Computation, SIAM, Philadelphia, 2019.
  * [11] S.-W. Duo and Y.-Z. Zhang, Mass-conservative Fourier spectral methods for solving the fractional nonlinear Schrödinger equation, Comput. Math. Appl., 71(11):2257-2271, 2016.
  * [12] S.-W. Duo, H. W. van Wyk and Y.-Z. Zhang, A novel and accurate finite difference method for the fractional Laplacian and the fractional Poisson problem, J. Comput. Phys., 355:233-252, 2018.
  * [13] G. H. Golub and C. F. van Loan, Matrix Computations, Fourth Edition, John Hopkins University Press, Baltimore, 2013.
  * [14] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products, Seventh Edition, Academic Press, 2007.
  * [15] J. R. Higgins, Completeness and Basis Properties of Sets of Special Functions, Cambridge University Press, Cambridge, 1977.
  * [16] C.-M. Huang and L.-W. Dong, Beam propagation management in a fractional Schrödinger equation, Sci. Rep., 7:5442, 2017.
  * [17] Y.-H. Huang and A. Oberman, Numerical methods for the fractional Laplacian: a finite difference-quadrature approach, SIAM J. Numer. Anal., 52(6):3056-3084, 2014.
  * [18] A. Iserles and M. Webb, A family of orthogonal rational functions and other orthogonal systems with a skew-Hermitian differentiation matrix, J. Fourier Anal. Appl., 26(1):No. 19, 2020.
  * [19] A. Iserles, K. Luong and M. Webb, Approximation of wave packets on the real line, arXiv:2101.02566, 2021.
  * [20] A. K. Kassam and L. N. Trefethen, Fourth-order time stepping for stiff PDEs, SIAM J. Sci. Comput., 26(4):1214-1233, 2005.
  * [21] J. Keilson, W. Nunn and U. Sumita, The bilateral Laguerre transform, Appl. Math. Comput., 8(2):137-174, 1981.
  * [22] K. Kirkpatrick, E. Lenzmann, G. Staffilani, On the continuum limit for discrete NLS with long-range lattice interactions, Comm. Math. Phys., 317(3):563-591, 2013.
  * [23] C. Klein, C. Sparber and P. Markowich, Numerical study of fractional nonlinear Schrödinger equations, Proc. Ser. A Math. Phys. Eng. Sci., 470(2172):20140364, 2014.
  * [24] S. Krogstad, Generalized integrating factor methods for stiff PDEs, J. Comput. Phys., 203(1):72-88, 2005.
  * [25] N. Laskin, Fractional Schrödinger equation, Phys. Rev. E, 66(5):056108, 2002.
  * [26] S. Longhi, Fractional Schrödinger equation in optics, Opt. Lett., 40(6):1117-1120, 2015.
  * [27] Z.-P. Mao and J. Shen, Hermite spectral methods for fractional PDEs in unbounded domains, SIAM J. Sci. Comput., 39(5):A1928-A1950, 2017.
  * [28] V. Minden and L.-X. Ying, A simple solver for the fractional Laplacian in multiple dimensions, SIAM J. Sci. Comput., 42(2):A878-A900, 2020.
  * [29] F. W. J. Olver, D. W. Lozier, R. F. Boisvert and C. W. Clark, NIST Handbook of Mathematical Functions, Cambridge University Press, Cambridge, UK, 2010.
  * [30] C.-T. Sheng, J. Shen, T. Tang, L.-L. Wang and H.-F. Yuan, Fast Fourier-like mapped Chebyshev spectral-Galerkin methods for PDEs with integral fractional Laplacian in unbounded domains, SIAM J. Numer. Anal., 58(5):2435-2464, 2020.
  * [31] C.-T. Sheng, S.-N. Ma, H.-Y. Li, L.-L. Wang and L.-L. Jia, Nontensorial generalised Hermite spectral methods for PDEs with fractional Laplacian and Schrödinger operators, ESIAM: M2NA, 55(5):2141-2168, 2021.
  * [32] B. A. Stickler, Potential condensed-matter realization of space-fractional quantum mechanics: The one-dimensional Lévy crystal, Phys. Rev. E, 88(1):012120, 2013.
  * [33] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal., 5(3):506-517, 1968.
  * [34] T. Tang, L.-L. Wang, H.-F. Yuan and T. Zhou, Rational spectral methods for PDEs involving fractional Laplacian in unbounded domains, SIAM J. Sci. Comput., 42(2):A585-A611, 2020.
  * [35] H. Weber, Numerical computation of the Fourier transform using Laguerre functions and the fast Fourier transform, Numer. Math., 36(2):197-209, 1980.
  * [36] J. A. C. Weideman, Computing the Hilbert transform on the real line, Math. Comp., 64(210):745-762, 1995.
  * [37] H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A, 150(5-7):262-268, 1990.