# Pseudospectral methods with PML for nonlinear Klein-Gordon equations in
classical and non-relativistic regimes

Xavier Antoine  X. Antoine: Université de Lorraine, CNRS, Inria, IECL, F-54000
Nancy, France<EMAIL_ADDRESS>and  Xiaofei Zhao  X. Zhao:
School of Mathematics and Statistics & Computational Sciences Hubei Key
Laboratory, Wuhan University, 430072 Wuhan, China<EMAIL_ADDRESS>
###### Abstract.

Two different Perfectly Matched Layer (PML) formulations with efficient
pseudo-spectral numerical schemes are derived for the standard and non-
relativistic nonlinear Klein-Gordon equations (NKGE). A pseudo-spectral
explicit exponential integrator scheme for a first-order formulation and a
linearly implicit preconditioned finite-difference scheme for a second-order
formulation are proposed and analyzed. To obtain a high spatial accuracy, new
regularized Bermúdez type absorption profiles are introduced for the PML. It
is shown that the two schemes are efficient, but the linearly implicit scheme
should be preferred for accuracy purpose when used within the framework of
pseudo-spectral methods combined with the regularized Bermúdez type functions.
In addition, in the non-relativistic regime, numerical examples lead to the
conclusion that the error related to regularized Bermúdez type absorption
functions is insensitive to the small parameter $\varepsilon$ involved in the
NKGE. The paper ends by a two-dimensional example showing that the strategy
extends to the rotating NKGE where the vortex dynamics is very well-
reproduced.

Keywords: nonlinear Klein-Gordon equation, perfectly matched layer, absorption
function, spectral accuracy, non-relativistic limit, pseudo-spectral method,
exponential integrator

AMS Subject Classification: 81Q05, 80M22, 37M15

##  1\. Introduction

The Klein-Gordon equation was originally proposed to describe the dynamics of
spinless particles, and is known as the relativistic version of the
Schrödinger equation. It has been widely applied in the studies of quantum
field theory, cosmology and plasma physics [16, 26, 41, 42]. In this work, we
consider the following $d$-dimensional ($d=1,2,3$) nonlinear Klein-Gordon
equation (NKGE) [15, 26, 34]:

$\left\\{\begin{split}&\frac{1}{c^{2}}\partial_{tt}u(\mathbf{x},t)-\Delta
u(\mathbf{x},t)+\frac{m^{2}c^{2}}{\hbar^{2}}u(\mathbf{x},t)+\lambda|u(\mathbf{x},t)|^{2}u(\mathbf{x},t)=0,\quad
t>0,\ \mathbf{x}\in{\mathbb{R}}^{d},\\\
&u(\mathbf{x},0)=u_{0}(\mathbf{x}),\quad\partial_{t}u(\mathbf{x},0)=v_{0}(\mathbf{x}),\quad\mathbf{x}\in{\mathbb{R}}^{d},\end{split}\right.$
(1.1)

where $c>0$ is the speed of light, $\hbar>0$ denotes the Planck constant,
$m>0$ is the particle mass and $\lambda>0$ is a given constant describing the
strength of the defocusing interaction. In addition, the unknown function is
$u=u(\mathbf{x},t):{\mathbb{R}}^{d}\times[0,\infty)\to{\mathbb{C}}$, where
$u_{0},v_{0}:{\mathbb{R}}^{d}\to{\mathbb{C}}$ are two given initial data.
Here, the global well-posedness of the defocusing NKGE (1.1) is guaranteed by
[28], while for the focusing case $\lambda<0$, there could be finite time
blow-up [22]. The following energy or Hamiltonian of (1.1) is also conserved

$H(t):=\int_{{\mathbb{R}}^{d}}\left[\frac{1}{c^{2}}|\partial_{t}u(\mathbf{x},t)|^{2}+|\nabla
u(\mathbf{x},t)|^{2}+\frac{m^{2}c^{2}}{\hbar^{2}}|u(\mathbf{x},t)|^{2}+\frac{\lambda}{2}|u(\mathbf{x},t)|^{4}\right]d\mathbf{x}\equiv
H(0),\quad t\geq 0.$ (1.2)

Since the initial-value problem (1.1) is set in ${\mathbb{R}}^{d}$, the
spatial domain has to be truncated to use a standard numerical discretization
scheme, e.g. finite-difference, finite-element or pseudo-spectral scheme [15].
In the present paper, our aim is to apply the Fourier pseudo-spectral
discretization scheme [38] which leads to highly accurate numerical solutions
of PDEs with smooth coefficients and is widely used in computations of quantum
mechanics [4, 5, 6, 7, 8, 11, 12, 13, 14, 15, 16, 27, 35]. As a consequence,
one needs to impose periodic boundary conditions on the fictitious boundary
delimiting the (rectangular) computational domain. Because of this constraint,
any truncation technique based on non-reflecting/artificial/absorbing boundary
condition for nonlinear PDEs [3, 9] cannot be applied. Since its introduction
in the seminal paper by Bérenger [18], the Perfectly Matched Layer (PML)
approach provides an alternative powerful tool to simulate the numerical
solution of PDEs in unbounded domains [1, 2, 3, 9, 10, 17, 19, 20, 21, 24, 25,
32, 43]. Concerning nonlinear PDEs, the PML technique has been applied for
example to the nonlinear Schrödinger equations [3, 9, 10, 43], Euler and
Navier-Stokes equations [31] or two-fluid plasma equations [40]. Concerning
the nonlinear wave equations, Appelö and Kreiss [2] proposed and studied a
first-order PML formulation with Dirichlet/Neumann boundary conditions related
to Hagstrom techniques [1, 29]. In addition, some first numerical experiments
were presented to analyze the potentiality of the approach which was also
tested for the one-dimensional linear Klein-Gordon equation in [9] again using
finite-differences. The aim of the present paper is to address the application
and assessment of the PML approach to solve complex NKGE, both in the
classical and non-relativistic regimes, considering Fourier pseudo-spectral
approximation schemes.

In Section 2, we introduce the first-order PML formulation in the classical
scaling, inspired by the works by Appelö and Kreiss in [2]. To obtain an
explicit pseudo-spectral scheme, we use an exponential-wave integrator
discretization. In addition, we also derive a second-order PML formulation of
the NKGE, related to the standard developments for time-harmonic wave and
Schrödinger equations. The discretization is now based on a Crank-Nicolson
scheme combined with the pseudo-spectral method and an efficient
preconditioned Krylov (GMRES) subspace iterative solver. This yields a
linearly implicit scheme with a cost similar to the one for the explicit
first-order PML formulation. Since we use a pseudo-spectral approach, the PDE
with PML that has to be discretized requires some smooth coefficients but also
some stability of the PML for the tuning parameters. To this end, we introduce
some generalized singular Bermúdez-type absorption functions which are locally
smoothed at the inner PML boundary. This allows us to achieve simultaneously a
high (near-spectral) accuracy of the scheme and stable PML layers that are
less sensitive to the tuning parameters, for the second-order formulation. The
first-order PML formulation is shown to be less accurate and is not
recommended in the pseudo-spectral framework. The conclusions are supported by
a thorough numerical study. In Section 3, we extend and evaluate the numerical
methods for the non-relativistic scaling. This shows that the second-order
formulation combined with the pseudo-spectral scheme is again an efficient and
accurate method for solving the NKGE, and most particularly thanks to the
small $\varepsilon$ parameter appearing in the non-relativistic regime. In
Section 4, we extend the method to a two-dimensional rotating NKGE that models
the dynamics of the cosmic superfluid set in a rotating frame [35, 42]. We
show in particular that the numerical method in the bounded domain is able to
simulate very accurately the dynamics of the vortices. Finally, we conclude in
Section 5.

##  2\. PML in classical scaling

To simplify the presentation, we first consider the one-dimensional real-
valued case of the NKGE, i.e. $d=1$, $\mathbf{x}=x$ and $u\in{\mathbb{R}}$ in
(1.1), and we begin by looking at its classical dimensionless form [28, 36]:

$\left\\{\begin{split}&\partial_{tt}u(x,t)-\partial_{xx}u(x,t)+u(x,t)+\lambda
u(x,t)^{3}=0,\quad t>0,\ x\in{\mathbb{R}},\\\
&u(x,0)=u_{0}(x),\quad\partial_{t}u(x,0)=v_{0}(x),\quad
x\in{\mathbb{R}}.\end{split}\right.$ (2.1)

Such classical scaling describes the physical system with wave speed at the
same order of the speed of light. We shall consider two types of PML
formulations for (2.1) in the sequel. One is a first-order formulation
proposed by Appelö and Kreiss in [2] for general nonlinear wave equations, and
the other one is a second-order analogy of the PML widely applied e.g. for
nonlinear Schrödinger equations [9, 10, 43].

###  2.1. PML-I

Let us begin by the modal ansatz construction with postulation

$\displaystyle
u(x,t)=\mathrm{e}^{kx+st}\mathrm{e}^{\frac{k}{s+\alpha}\int_{0}^{x}\sigma(\rho)d\rho},\quad\partial_{x}\to\frac{s+\alpha}{s+\alpha+\sigma}\partial_{x},$
(2.2)

where $\sigma=\sigma(x)\geq 0$ is some chosen function known as the
_absorption function_ , $s$ is interpreted as the variable in the Laplace-
transform domain, and $\alpha\geq 0$ is a chosen parameter. By assuming that
${\mathrm{Re}}(s)\geq 0$ and $k\in{\mathbb{C}}$ in (2.2) satisfy the
dispersion relation for the linear version of the model, e.g.
$\partial_{tt}u-\partial_{xx}u+u=0$, a formulation of the PML equations as a
first-order system has been proposed for a class of general nonlinear wave
equations in [2]. Here, we can directly apply such PML formulation (called
PML-I) to the NKGE (2.1), leading to

$\left\\{\begin{split}&\partial_{tt}u-\partial_{xx}u+u+\lambda
u^{3}=\sigma\left[\eta_{2}-\partial_{t}u+\alpha
u\right]+\partial_{x}(\sigma\eta_{1}),\quad t>0,\ x\in I^{*},\\\
&\partial_{t}\eta_{1}+(\sigma+\alpha)\eta_{1}+\partial_{x}u=0,\\\
&\partial_{t}\eta_{2}+\alpha\eta_{2}+(\alpha^{2}+1)u+\lambda u^{3}=0,\\\
&u(x,0)=u_{0}(x),\quad\partial_{t}u(x,0)=v_{0}(x),\quad\eta_{1}(x,0)=\eta_{2}(x,0)=0,\quad
x\in I^{*},\\\
&u(-L^{*},t)=u(L^{*},t),\quad\eta_{1}(-L^{*},t)=\eta_{1}(L^{*},t),\quad\eta_{2}(-L^{*},t)=\eta_{2}(L^{*},t),\quad
t\geq 0,\end{split}\right.$ (2.3)

where

$I^{*}=(-L^{*},L^{*})\quad\mbox{with}\quad L^{*}=L+\delta$

is a bounded interval for some $\delta>0$, and $\eta_{1}=\eta_{1}(x,t),\
\eta_{2}=\eta_{2}(x,t)$ are two auxiliary functions introduced in order to
localize the layer equations in time. The initial conditions for the two
auxiliary functions in (2.3) are suggested to be zero in [43].

The absorption function $\sigma(x)$ is required to satisfy

$\sigma(x)=0,\quad-L\leq x\leq L,\quad\mbox{and}\quad\sigma(x)>0,\quad
L<|x|\leq L^{*},$

such that inside the physical domain $[-L,L]$ the NKGE and its solution remain
the same. The damping effect takes place only inside the layer, i.e.
$L<|x|\leq L^{*}$. The smoothness of $\sigma(x)$ in the layer region
$L\leq|x|\leq L^{*}$, particularly at the interface $|x|=L$, determines the
regularity of the PML solution $u(x,t)$. When the NKGE (2.1) is linear
($\lambda=0$), the absorbing layer (2.3) is perfectly matching [2]. Similar
PMLs of the first-order formulation have also been derived for other
hyperbolic systems [1, 32]. The boundary conditions for (2.3) at $|x|=L^{*}$
can be for example the homogeneous Dirichlet or Neumann boundary condition
which were considered in [2]. In either case, the reflected waves at the
boundary enter again the layer region and so are damped twice. Here we adopt
periodic boundary conditions as suggested in [9, 10], so that the Fourier
pseudo-spectral method [38] can be easily applied for the spatial
discretization of the PML equations. Note that the periodicity of the boundary
ensures that the waves hitting the outer boundary enter the layer from the
other side of the domain which is the same as the ‘round-way’ damping.

Let us discuss briefly the damping effect brought by (2.2) for the Klein-
Gordon equation and the choice of $\alpha$. Firstly, we remark that with
$\alpha=0$, the damping effect

$\frac{k}{s}\int_{0}^{x}\sigma(\rho)d\rho$

in (2.2) is a typical choice for dispersionless materials, where the phase
velocity $k/s$ is constant so that all waves attenuate at the same rate [32].
For the PML-I formulation for the NKGE (2.3), we consider the linear case,
i.e. $\lambda=0$ in (2.1) or (2.3) to illustrate the damping effect. By taking
the modal solution $u=\mathrm{e}^{kx+st}$, we have the dispersion relation for
the original linear Klein-Gordon equation:

$s^{2}=k^{2}-1\ \Rightarrow\ k=-\sqrt{s^{2}+1},$ (2.4)

and so (2.2) reads

$u(x,t)=\mathrm{e}^{kx+st}\mathrm{e}^{-\frac{\sqrt{s^{2}+1}}{s+\alpha}\int_{0}^{x}\sigma(\rho)d\rho}.$

This means that the modal solution $\mathrm{e}^{kx+st}$ for PML-I (2.3) is
damped inside the layer $L<|x|\leq L^{*}$ by the absorption function $\sigma$
modulated by the factor $g(s)$:

$g(s)=-\frac{\sqrt{s^{2}+1}}{s+\alpha},\quad\mbox{for}\quad{\mathrm{Re}}(s)>0.$
(2.5)

We show in Figure 1 the mapping of the right half plane of the complex domain
${\mathrm{Re}}(s)>0$ under the function $g(s)$ for $\alpha=0$ and
$\alpha=0.5$. We can see that for ${\mathrm{Re}}(s)>0$ and $\alpha\geq 0$, we
always have ${\mathrm{Re}}(g)<0$. However, it is difficult to say if there
exists or not an optimal choice of the shifting parameter $\alpha$ even for
such a linear equation. The reason is that there could always be some $g(s)$
getting arbitrarily close to the imaginary axis for some $s$, with
${\mathrm{Re}}(s)$ arbitrarily small.

$\begin{array}[]{c}\psfig{height=113.81102pt,width=170.71652pt}\psfig{height=113.81102pt,width=170.71652pt}\\\
\psfig{height=113.81102pt,width=170.71652pt}\end{array}$ Figure 1. The right
half complex plane ${\mathrm{Re}}(s)>0$ (top left) under mapping
$g=-\sqrt{s^{2}+1}/(s+\alpha)$: $\alpha=0$ (top right) and $\alpha=0.5$
(bottom).

To numerically solve the PML-I equations (2.3), instead of the finite
difference time-domain discretization in [2], our periodic boundary setup in
(2.3) is natural to apply the Fourier spectral discretization which is
directly diagonalized in space, and we integrate in time under the framework
of the exponential wave integrator [14, 27, 30]. Denoting $\tau=\Delta t>0$ as
the time step and $t_{n}=n\tau$ for $n\in{\mathbb{N}}$ as the discrete times,
we shall derive in the following a fully explicit exponential wave integrator
Fourier pseudo-spectral method for (2.3). For simplicity of notations, we will
omit in the following the spatial variable $x$ in a space-time dependent
function, e.g. $u(t)=u(x,t)$. Let

$v(t):=\partial_{t}u(t),\quad t\geq 0.$

Then, by applying the Duhamel’s formula to (2.3), we have

$\displaystyle u(t_{n+1})=$
$\displaystyle\cos(\langle\partial_{x}\rangle\tau)u(t_{n})+\frac{\sin(\langle\partial_{x}\rangle\tau)}{\langle\partial_{x}\rangle}v(t_{n})+\int_{0}^{\tau}\frac{\sin(\langle\partial_{x}\rangle(\tau-\rho))}{\langle\partial_{x}\rangle}f(t_{n}+\rho)d\rho,$
$\displaystyle v(t_{n+1})=$
$\displaystyle-\langle\partial_{x}\rangle\sin(\langle\partial_{x}\rangle\tau)u(t_{n})+\cos(\langle\partial_{x}\rangle\tau)v(t_{n})+\int_{0}^{\tau}\cos(\langle\partial_{x}\rangle(\tau-\rho))f(t_{n}+\rho)d\rho,$
$\displaystyle\eta_{1}(t_{n+1})=$
$\displaystyle\mathrm{e}^{-(\sigma+\alpha)\tau}\eta_{1}(t_{n})-\int_{0}^{\tau}\mathrm{e}^{(\sigma+\alpha)(\rho-\tau)}\partial_{x}u(t_{n}+\rho)d\rho,$
$\displaystyle\eta_{2}(t_{n+1})=$
$\displaystyle\mathrm{e}^{-\alpha\tau}\eta_{2}(t_{n})-\int_{0}^{\tau}\mathrm{e}^{\alpha(\rho-\tau)}\left[(\alpha^{2}+1)u(t_{n}+\rho)+\lambda
u(t_{n}+\rho)^{3}\right]d\rho,\quad n\geq 0,\ x\in I^{*},$

where we set $\langle\partial_{x}\rangle=\sqrt{1-\partial_{xx}}$ and

$\displaystyle f(t_{n}+\rho):=\sigma(\eta_{2}(t_{n}+\rho)-v(t_{n}+\rho)+\alpha
u(t_{n}+\rho))+\partial_{x}(\sigma\eta_{1}(t_{n}+\rho))-\lambda
u(t_{n}+\rho)^{3}.$

Applying the trapezoidal rule to the above integrals, we obtain the following
exponential-wave integrator Fourier pseudo-spectral (EWI-FP) scheme: by
introducing $u^{n}\approx u(t_{n}),v^{n}\approx
v(t_{n}),\eta_{1}^{n}\approx\eta_{1}(t_{n}),\eta_{2}^{n}\approx\eta_{2}(t_{n})$,
for $n\geq 0$, we have

$\displaystyle u^{n+1}=$
$\displaystyle\cos(\langle\partial_{x}\rangle\tau)u^{n}+\frac{\sin(\langle\partial_{x}\rangle\tau)}{\langle\partial_{x}\rangle}v^{n}+\frac{\tau\sin(\langle\partial_{x}\rangle\tau)}{2\langle\partial_{x}\rangle}f^{n},$
(2.6a) $\displaystyle\eta_{1}^{n+1}=$
$\displaystyle\mathrm{e}^{-(\sigma+\alpha)\tau}\eta_{1}^{n}-\frac{\tau}{2}\left[\mathrm{e}^{-(\sigma+\alpha)\tau}\partial_{x}u^{n}+\partial_{x}u^{n+1}\right],$
(2.6b) $\displaystyle\eta_{2}^{n+1}=$
$\displaystyle\mathrm{e}^{-\alpha\tau}\eta_{2}^{n}-\frac{\tau}{2}\left[\mathrm{e}^{-\alpha\tau}\left(\left(\alpha^{2}+1\right)u^{n}+\lambda(u^{n})^{3}\right)+\left(\alpha^{2}+1\right)u^{n+1}+\lambda\left(u^{n+1}\right)^{3}\right],$
(2.6c) $\displaystyle v^{n+1}=$
$\displaystyle-\langle\partial_{x}\rangle\sin(\langle\partial_{x}\rangle\tau)u^{n}+\cos(\langle\partial_{x}\rangle\tau)v^{n}+\frac{\tau}{2}\left[\cos(\langle\partial_{x}\rangle\tau)f^{n}+f^{n+1}\right],$
(2.6d)

with

$f^{n}=\sigma\left(\eta_{2}^{n}-v^{n}+\alpha
u^{n}\right)+\partial_{x}(\sigma\eta_{1}^{n})-\lambda(u^{n})^{3}.$

Here, the scheme (2.6) and the numerical solution
$(u^{n},v^{n},\eta_{1}^{n},\eta_{2}^{n})$ are defined on the spatial grids:
$x_{j}=-L^{*}+jh$ for $j=0,1,\ldots,N-1$ with some even integer $N>0$ and mesh
size $h=2L^{*}/N$.

The above EWI-FP scheme is fully explicit in time by computing (2.6a) to
(2.6d). The involved spatial differentiation operators can all be implemented
efficiently under Fourier pseudo-spectral method [38] by fast Fourier
transform (FFT), resulting in a computational cost $O(N\log N)$. Concerning
the stability of the EWI-FP scheme, noting the term $\partial_{x}u$ in (2.6b)
and the term $\partial_{x}(\sigma\eta_{1}^{n})$ appearing in (2.6a) and
(2.6d), we have some unbalanced energy norms on both sides of (2.6).
Therefore, we are expecting a CFL stability constraint:

$\tau<Ch,\quad\mbox{for some}\quad C>0.$

For the accuracy, noticing the quadrature error from the trapezoidal rule and
the interpolation error from the Fourier pseudo-spectral method, it can be
analyzed similarly as in [14, 27] that the global error of the above EWI-FP
method for (2.3) up to some finite time reads

$O(\tau^{2})+O(h^{m_{0}}),$

for some $m_{0}>0$ which depends on the smoothness of the PML solution of
(2.3). We shall show numerically that by using a very smooth absorption
function $\sigma(x)$ particularly with high-order continuous derivatives at
the interface point $|x|=L$, the EWI-FP could reach a spectral-like accuracy
for spatial discretization. Thanks to the spectral accuracy, we can use a
rather large $h>0$ in practice so that the above CFL condition is not too
restrictive.

###  2.2. PML-II

For dispersive models like the nonlinear Schrödinger equations (NSE), the
phase velocity is naturally non-constant with respect to the wavelength.
Therefore, instead of (2.2), an average damping effect has been widely
considered in the postulation of a modal solution:

$u(x,t)=\mathrm{e}^{kx+st}\mathrm{e}^{kR\int_{0}^{x}\sigma(\rho)d\rho},$ (2.7)

for some chosen $R\in{\mathbb{C}}$ as for the NSE. Using (2.7) for the NKGE
(2.1), by directly changing

$\partial_{x}\rightarrow\frac{1}{1+R\sigma}\partial_{x},$

we consider the following second PML formulation (PML-II):

$\left\\{\begin{split}&\partial_{tt}u-\frac{1}{1+R\sigma}\partial_{x}\left(\frac{1}{1+R\sigma}\partial_{x}u\right)+u+\lambda
u^{3}=0,\quad t>0,\ x\in I^{*},\\\
&u(x,0)=u_{0}(x),\quad\partial_{t}u(x,0)=v_{0}(x),\quad x\in I^{*},\\\
&u(-L^{*},t)=u(L^{*},t),\quad t\geq 0.\end{split}\right.$ (2.8)

Here, we also impose the periodic boundary condition in (2.8) for the
convenience of spatial discretization.

An obvious advantage of PML-II (2.8) over PML-I (2.3) is that PML-II does not
require any auxiliary variables. We only need to solve one equation in (2.8),
instead of three equations in (2.3). The PML-II can be considered as a second-
order formulation, which is an analogy of the PML e.g. for NSE [9, 10, 43].
Note that for the NSE, the shifting parameter $R$ is usually chosen as
$R=\mathrm{e}^{i\pi/4}$ [43]. However, for the NKGE (2.8), this choice
generates some instabilities, and we need to consider $R>0$. Indeed, assuming
again $\lambda=0$ in (2.1) with the modal solution $u=\mathrm{e}^{kx+st}$, by
the dispersion relation (2.4), we have in (2.7):

$u(x,t)=\mathrm{e}^{kx+st}\mathrm{e}^{-R\sqrt{s^{2}+1}\int_{0}^{x}\sigma(\rho)d\rho}.$

Then in order to make the solution decay inside the layer, we ask for

${\mathrm{Re}}\left(-R\sqrt{s^{2}+1}\right)<0,\quad\forall\
{\mathrm{Re}}(s)>0.$

We illustrate in Figure 2 the mapping $g(s)=-\sqrt{s^{2}+1}$ of the right half
complex plane. It can be clearly seen that the choice of $R$ for the PML (2.8)
is restricted to a positive real number, otherwise any non-zero imaginary part
of $R$ leads to some ${\mathrm{Re}}\left(-R\sqrt{s^{2}+1}\right)>0$ which
triggers instability. This will be verified by numerical tests later. Since
$R$ must be positive, it can be viewed in (2.8) only as a modification of
$\sigma_{0}>0$, namely the strength of the absorption function $\sigma(x)$
(see e.g. (2.13) or (2.15)). Therefore, we will only consider the impact from
the choice of $\sigma_{0}>0$ on the performance of the PML-II with $R>0$ fixed
in this section.

$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 2. The right half complex plane ${\mathrm{Re}}(s)>0$ (left) under the
mapping $g=-\sqrt{s^{2}+1}$ (right).

For the numerical discretization of the PML-II formulation (2.8), we use the
following direct semi-implicit finite-difference Fourier pseudo-spectral (FD-
FP) method:

$\displaystyle\frac{u^{n+1}-2u^{n}+u^{n-1}}{\tau^{2}}+\frac{A}{2}(u^{n+1}+u^{n-1})+\frac{1}{2}(u^{n+1}+u^{n-1})+\lambda(u^{n})^{3}=0,\quad
n\geq 1,$ (2.9a) $\displaystyle u^{1}=u_{0}+\tau
v_{0}-\frac{\tau^{2}}{2}\left[Au_{0}+u_{0}+\lambda u_{0}^{3}\right],$ (2.9b)

where $u^{n}\in{\mathbb{R}}^{N}$ and the matrix

$A=-d_{0}D_{1}d_{0}D_{1},\quad\mbox{with}\quad
d_{0}=\mathrm{diag}(1/(1+R\sigma)),$

and $D_{1}\in{\mathbb{C}}^{N\times N}$ is the matrix representation [38] of
the approximate first-order differentiation operator $\partial_{x}$ in the
periodic domain $I^{*}$. The time averaging in (2.9a) indeed avoids any CFL
type stability condition. The starting value (2.9b) is obtained directly from
the Taylor expansion for the three-level scheme (2.9a). Note that (2.9a) reads
explicitly

$\displaystyle
u^{n+1}=-u^{n-1}+G^{-1}\left[\frac{2}{\tau^{2}}u^{n}-\lambda(u^{n})^{3}\right],\quad
n\geq 1,\quad
G:=\left[\left(\frac{1}{\tau^{2}}+\frac{1}{2}\right)\mathbb{I}+\frac{1}{2}A\right],$
(2.10)

where $\mathbb{I}$ is the identity matrix in ${\mathbb{R}}^{N\times N}$. The
matrix $G^{-1}$ is dense here, but it can be pre-computed with the chosen
$\tau$ and $N$ once for all. The computational cost of FD-FP at each time
level is therefore $O(N^{2})$ if one performs the previous matrix-vector
product. However, by borrowing the idea from [10], the scheme (2.10) can be
implemented efficiently by incorporating the GMRES solver with FFT. Here we
briefly describe such a strategy for evaluating (2.10). The equation (2.10)
equivalently reads

$\left\\{\begin{array}[]{l}\displaystyle u^{n+1}=-u^{n-1}+w^{n},\quad n\geq
1,\\\ \displaystyle
Gw^{n}=\frac{2}{\tau^{2}}u^{n}-\lambda(u^{n})^{3}.\end{array}\right.$ (2.11)

The second equation of (2.11) can be treated as a linear system with unknown
$w^{n}$, which can be solved iteratively e.g. by GMRES [37] since it is a
matrix-free solver. Within the GMRES iteration, the involved spatial
derivatives in $G$ can be straightforwardly implemented by applying the FFT to
discretize the operators in the iteration process. The GMRES convergence can
be strongly improved by introducing a Fourier diagonalizable preconditioner
for $\sigma=0$, similarly as in [5, 6, 7], i.e.

$\mathcal{P}Gw^{n}=\mathcal{P}\left[\frac{2}{\tau^{2}}u^{n}-\lambda(u^{n})^{3}\right],$

with

$\mathcal{P}=\left(\frac{1}{\tau^{2}}+\frac{1}{2}-\frac{\partial_{xx}}{2}\right)^{-1}.$
(2.12)

Let us remark that the preconditioning operator $\mathcal{P}$ can also be
directly implemented at low cost by FFT, since it has constant coefficients
and is therefore diagonalizable. The computational cost of the FD-FP method is
then nearly $O(N\log N)$. The global error of FD-FP is similar to EWI-FP,
which is $O(\tau^{2})+O(h^{m_{0}})$.

###  2.3. Absorption function

In this subsection, we discuss the choices for the absorption function. To
enhance the efficiency of the PMLs, we look for spectral accuracy in the
related spatial discretization. Let us recall that getting a smooth absorption
function is crucial to get the spectral accuracy, where the parameter
$m_{0}>0$ specifies the accuracy induced by the PML absorption function on the
error estimates. We shall consider two types of absorption functions
$\sigma(x)$ for the PMLs: a classical polynomial choice and a singular
function case.

The first possibility is to choose a high degree _polynomial type_ function as
proposed in [2] for wave equations

$\sigma_{P}(x):=\left\\{\begin{split}&\sigma_{0}\left[1-\left(\frac{x-L^{*}}{\delta}\right)^{2}\right]^{8},\quad
L\leq|x|\leq L^{*}=L+\delta,\\\ &0,\qquad\qquad\qquad\qquad\qquad\
\mbox{else}.\end{split}\right.$ (2.13)

It can be viewed as a smoothed version of the commonly used quadratic or cubic
absorbing functions for PMLs [10, 18, 24, 43]. Thanks to the large power,
(2.13) is very smooth at the interfaces $|x|=L$ so that one can maintain the
high order accuracy from the Fourier spectral discretization in space.

Nevertheless, it is well-known since the work by Bermúdez et al. [19, 20, 21]
that singular type functions are much more adapted than the polynomial
functions for PMLs applied to wave-like equations. More specifically, the
$\sigma$ functions must be such that

$\int_{L}^{L^{*}}\sigma(x)dx=+\infty.$ (2.14)

In particular, they are less subject to parameter tuning problems than in the
polynomial case (2.13) to fix the values of $\sigma_{0}$ and $\delta$, in
particular thanks to the speed $c$. This is crucial in our case when
prospecting the non-relativistic limit of the NKGE in the next section.

The concrete examples provided by Bermúdez et al. do not correspond however to
smooth functions, and therefore they lead to a limited accuracy for the
pseudo-spectral approach. Here, we propose some simple locally corrected
smoother _Bermúdez type_ functions as follows: for $k=-1,0,1,\ldots$, we
define the absorption functions

$\sigma_{B_{k}}(x):=\left\\{\begin{split}&\sigma_{0}\beta_{k}(|x|-L^{*}),\quad
L\leq|x|\leq L^{*}=L+\delta,\\\
&0,\qquad\qquad\qquad\quad\mbox{else},\end{split}\right.$ (2.15)

where we set

$\displaystyle\beta_{k}(|x|-L^{*})=\beta_{-1}(|x|-L^{*})-\sum_{j=0}^{k}\frac{1}{j!}\frac{d^{j}\beta_{-1}(z)}{dz^{j}}\bigg{|}_{z=-\delta}(|x|-L)^{j},\quad\beta_{-1}(z)=-\frac{1}{z},\quad
z\in{\mathbb{R}}.$

The function $\sigma_{B_{-1}}$ (i.e. $k=-1$ in (2.15)) is the original
absorption function proposed in [19, 20, 21]. It contains both discontinuity
and singularity over the domain. For $k\geq 0$ as above, we have introduced a
truncated Taylor series expansion of $\beta_{-1}$ in $\beta_{k}$ so that the
function $\sigma_{B_{k}}(x)$ has $k$-th order continuous derivatives on the
interval $|x|<L^{*}$, which is expected for the spectral scheme.

The polynomial type choice $\sigma(x)=\sigma_{P}(x)$ is the classical
absorption function which is bounded on the whole interval $|x|\leq L+\delta$.
Here, $\sigma_{0}=\|\sigma_{P}\|_{\infty}>0$ and $\delta>0$ are known as the
_strength_ and _thickness_ of the PML, respectively. The function $\sigma_{P}$
can be directly applied to both PML formulations, i.e. PML-I (2.3) and PML-II
(2.8). The Bermúdez type choice $\sigma(x)=\sigma_{B_{k}}(x)$ has
singularities at the two boundary points $|x|=L^{*}$. The strength of such
absorption function can be considered as infinite in some sense. This indeed
is a problem for the PML-I formulation (2.3), where in the equations we can
see that the infinite function value of $\sigma_{B_{k}}(x)$ at the boundary is
inconsistent with the periodic boundary condition mathematically. Moreover,
when the propagating waves are getting close to the boundary, the large values
of $\sigma_{B_{k}}(x)$ near the boundary could make the right-hand-side of the
first equation in (2.3) stiff which causes numerical instability problems in
schemes like EWI-FP. Thus, for PML-I (2.3), we shall only consider the
polynomial absorption function $\sigma_{P}$ (2.13). As for the PML-II
formulation (2.8), in contrast, $\sigma_{B_{k}}(x)$ can be directly applied in
the equation, where we have a well-defined function $S(x)$ on $I^{*}$ with
zero boundary values, i.e. in (2.8) with $\sigma=\sigma_{B_{k}}$ for any
$k\in{\mathbb{N}}$,

$S(x):=\frac{1}{1+R\sigma_{B_{k}}(x)},\quad|x|\leq
L^{*},\quad\mbox{where}\quad S(\pm L^{*})=0.$

Therefore, we shall consider both types of absorption functions for PML-II
(2.8).

###  2.4. Numerical results

In this subsection, we conduct some numerical experiments to illustrate the
performance of the two presented types of PML formulations as well as the
proposed numerical schemes. We shall denote in the following $u=u(x,t)$ as the
exact solution of the NKGE (2.1), $u_{\textrm{pml}}=u_{\textrm{pml}}(x,t)$ as
the exact solution of the PML-I (2.3) or PML-II (2.8), and
$u_{\textrm{pml}}^{n}\approx u_{\textrm{pml}}(x,t_{n})$ as the corresponding
numerical solution from EWI-FP (2.6) or FD-FP (2.9).

PML-I. We begin with the PML-I formulation (2.3). Within these tests, we fix
the absorption function as the high order polynomial type function (2.13),
i.e.

$\sigma(x)=\sigma_{P}(x),\quad|x|\leq L^{*}.$

First of all, we verify the accuracy of our numerical solver EWI-FP (2.6) for
approximating (2.3) with

$\lambda=1,\quad u_{0}=5\mathrm{e}^{-x^{2}},\quad
v_{0}=\frac{1}{2}\mathrm{sech}(x^{2}),$ (2.16)

and fix $\alpha=0$, $\sigma_{0}=8,\,\delta=0.5,\,L=4$ for the layer. Then we
compute the relative maximum error

$e^{n,\infty}_{\textrm{pml}}:=\frac{\|u_{\textrm{pml}}^{n}-u_{\textrm{pml}}\|_{L^{\infty}(I)}}{\|u_{\textrm{pml}}\|_{L^{\infty}(I)}}$
(2.17)

on the physical domain $I=(-L,L)$ at some time $t=t_{n}$, where the reference
solution $u_{\textrm{pml}}$ of (2.3) is obtained here numerically by using
EWI-FP under a fine mesh with $\tau=1\times 10^{-4}$ and $h=1/128$. The
temporal and spatial errors of the EWI-FP scheme (2.6) at $t=4$ are shown in
Figure 3, where we fix $h=1/32$ for the temporal error (left) and
$\tau=1\times 10^{-4}$ for the spatial error (right). From the numerical
results in Figure 3, we can clearly see that the proposed EWI-FP (2.6) for
solving (2.3) with (2.13) converges in time with a second-order accuracy and
converges in space with a near spectral accuracy.

$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 3. The temporal (left) and spatial (right)
$e^{n,\infty}_{\textrm{pml}}$-errors of EWI-FP for PML-I.

Next, we test the error of the PML-I formulation (2.3) for approximating the
solution to the NKGE (2.1). We use EWI-FP with a very fine mesh, e.g.
$\tau=1\times 10^{-4}$ and $h=1/128$, such that the numerical discretization
error is rather negligible. We compute the following relative PML error as a
function of time $t>0$

$e^{2}_{\textrm{pml}}(t):=\frac{\|u(\cdot,t)-u_{\textrm{pml}}(\cdot,t)\|_{L^{2}(I)}}{\|u(\cdot,t)\|_{L^{2}(I)}},\quad
I=(-L,L),$ (2.18)

under some different strength $\sigma_{0}$ or thickness $\delta$ parameters
for the absorption function (2.13). We use the same numerical example (2.16)
as above, and we measure the error (2.18) on the physical domain $I=(-4,4)$.
The exact solution $u(x,t)$ of (2.1) is obtained by using the EWI-FP scheme
(for $u$ only) to solve (2.1) accurately on a large enough interval, e.g.
$(-16,16)$, where the solution within the computational time is still away
from the boundary. The relative errors (2.18) as a function of time under
different parameters for (2.13) are plotted in Figure 4. To compare the
performance of the PML-I (2.3) in the nonlinear case of the Klein-Gordon
equation with its linear case, we also include in Figure 4 the error result
(2.18) under the same example (2.16) but with $\lambda=0$ in (2.1) and (2.3).
For (2.16), the profiles of the exact solution of (2.1) and the PML solution
of (2.3) at time $t=4$ and $t=6$ are reported in Figure 5 for the PML
parameters $\sigma_{0}=8$ and $\delta=7/8$.

Linear case: $\lambda=0$

$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$

Nonlinear case: $\lambda=1$

$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$

Figure 4. The PML-I $e^{2}_{\textrm{pml}}$-error as a function of time for
$\sigma=\sigma_{P}$ in the linear and nonlinear cases: $\delta=1/2$ and
varying strength $\sigma_{0}$ (left), and $\sigma_{0}=6$ and varying thickness
$\delta$ (right).
$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 5. The exact solution and the PML-I solution at $t=4$ (left) and $t=6$
(right).

From the numerical results in Figures 4 and 5, we can see that the PML-I (2.3)
is effectively approximating the solution of NKGE (2.1) within the physical
domain (see Figure 5). The PML-I formulation in the linear case, i.e.
$\lambda=0$ in (2.3), was proposed to be perfectly matched for the linear
Klein-Gordon equation in [2]. In the nonlinear case, i.e. $\lambda\neq 0$ in
NKGE (2.1), although PML-I (2.3) is not theoretically perfectly matched, its
performance is very close to the linear case (cf. the first row and the second
row in Figure 4). With increasing strength of the absorption function and
thickness of the fictitious layer, the accuracy of PML-I could be improved
(see Figure 4).

######  Remark 2.1.

With a usual low order polynomial absorption function $\sigma$, e.g. the
quadratic or cubic polynomial in [10, 18, 24, 43], the performance of the
PML-I (2.3) is similar to that in Figure 4. However, the spatial convergence
rate of the EWI-FP method would only be at the second- or third-order. These
numerical results are omitted here for brevity.

######  Remark 2.2.

The choice of $\alpha\geq 0$ barely impacts the performance of PML-I (2.3),
where the PML errors (2.18) when using different $\alpha$ are very close. This
is not surprising within the context of the analysis of the damping effect in
(2.5). Therefore, we only present the numerical results for $\alpha=0$.

PML-II. Under the same numerical example (2.16), we present the numerical
experiments for PML-II given by (2.8). The physical domain is fixed as
$I=(-4,4)$. For these tests, we will consider both the polynomial choice
(2.13) and the Bermúdez type function (2.15) as the absorption function for
PML-II.

Firstly, we take the polynomial choice (2.13) as the absorption function, i.e.
$\sigma(x)=\sigma_{P}(x)$ for PML-II. Similarly as before, with $\sigma_{0}=8$
and $\delta=0.5$ fixed for the absorption function, we first show the temporal
and spatial accuracy of the numerical method, i.e. FD-FP (2.9) at $t=4$ in
Figure 6. Next, by taking in FD-FP a very fine mesh ($\tau=10^{-4}$ and
$h=1/256$) so that the numerical discretization error is negligible, the
corresponding PML error (2.18) of the PML-II formulation (2.8) is shown in
Figure 7 as a function of time for different values of $\sigma_{0}$ and
$\delta$ (for $R=1$ in (2.8)).

$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 6. The temporal (left) and spatial (right)
$e^{n,\infty}_{\textrm{pml}}$-errors of FD-FP for PML-II with
$\sigma=\sigma_{P}$.
$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 7. The PML-II $e^{2}_{\textrm{pml}}$-error as a function of time under
$\sigma=\sigma_{P}$: for $\delta=1/2$ and varying strength $\sigma_{0}$
(left), and for $\sigma_{0}=6$ and varying thickness $\delta$ (right).

As explained in Section 2.2 by the dispersion relation in the linear case, the
factor $R$ must be positive to avoid any instability. To confirm the formal
analysis and illustrate the stability issue in the PML-II equation (2.8) with
respect to $R$, we show in Figure 8 the evolution of the maximum norm of the
PML solution $u_{\textrm{pml}}$ of (2.8) for $R>0$ or $R\in{\mathbb{C}}$.

$\begin{array}[]{c}\psfig{height=113.81102pt,width=156.49014pt}\psfig{height=113.81102pt,width=156.49014pt}\\\
\psfig{height=113.81102pt,width=156.49014pt}\end{array}$ Figure 8. Stability
of PML-II: $\|u_{\textrm{pml}}\|_{L^{\infty}}$ as a function of time for $R=1$
(top left), $R=\mathrm{e}^{i\pi/4}$ (top right), $R=\mathrm{e}^{i\pi/16}$
(bottom).

From the numerical results in Figures 6-8, we can conclude that:

  * 1)

The PML-II formulation (2.8) with $R>0$ is stable (see Figure 8), while, for
$R$ with any non-zero imaginary part, instability will occur in (2.8).

  * 2)

The FD-FP method (2.9) is also second-order accurate in time and nearly
spectrally accurate in space. Under the same mesh size and time step, it is
slightly less accurate than the EWI-FP method (cf. Figures 6 and 3).
Therefore, from the numerical discretization point of view, PML-II is a little
bit less efficient than the PML-I.

  * 3)

PML-II effectively approximates the exact solution of NKGE (2.1) (see Figure
7). Moreover, under the same absorption function $\sigma(x)$ with the same
strength and thickness, the error from the PML-II is smaller than PML-I
particularly as time increases (cf. Figures 7 and 4). Therefore, from the
modelling truncation point of view, PML-II is more effective than PML-I. We
remark that the results of using the usual low order polynomial absorption
function for PML-II is similar.

Next, we choose the Bermúdez type function (2.15) as the absorption function,
i.e. $\sigma(x)=\sigma_{B_{k}}(x)$ for PML-II (2.8). We use the same example
(2.16) with $\sigma_{0}=8,\delta=0.5$ in (2.15), and we first test the
numerical discretization error of the FD-FP method (2.9) for PML-II (2.8) with
$R=1$. The spatial error at $t=4$ is presented in Figure 9 under different
smoothing orders $k\geq 0$ in the Bermúdez function (2.15). The corresponding
temporal error is totally the same as in the previous example with the
polynomial choice (2.13), so it is omitted here for brevity. On the other
hand, to illustrate the efficiency brought by the preconditioner (2.12) for
solving (2.11), we fix in this example $\sigma=\sigma_{B_{2}},\ \tau=0.02$ and
consider the GMRES solver for the time step $n=1$ in (2.11). We show the
number of iterations needed by the GMRES solver (without restart) to converge
at a threshold $\epsilon>0$ with or without using the preconditioner (2.12)
under different values of $h$ in Table 1. This confirms that the
preconditioner leads to a convergence independent of the mesh refinement,
which is well-adapted to the pseudo-spectral method. The corresponding PML
error (2.18) with respect to time in this example is presented in Figure 10
for different values of the parameters $\sigma_{0}$ and $\delta$.

$\begin{array}[]{c}\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 9. Spatial $e^{n,\infty}_{\textrm{pml}}$-error of FD-FP for PML-II with
$\sigma=\sigma_{B_{k}}$, for different values of $k$.
$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 10. The PML-II $e^{2}_{\textrm{pml}}$-error as a function of time for
$\sigma=\sigma_{B_{2}}$: for $\delta=1/2$ and varying strength $\sigma_{0}$
(left), and for $\sigma_{0}=6$ and varying thickness $\delta$ (right). Table
1. Number of iterations needed by GMRES under threshold $\epsilon$ for (2.11)
at $n=1$: with (pre) or without (non) the preconditioner.

$\epsilon=10^{-10}$ $h=1/128$ $h=1/256$ $h=1/512$ $\epsilon=10^{-13}$
$h=1/128$ $h=1/256$ $h=1/512$ pre 2 2 2 pre 8 7 7 non 11 20 35 non 24 46 92

Based on the numerical results in Table 1, Figures 9 and 10, we can draw the
following observations:

  * 1)

For a high enough regularization parameter $k\geq 0$, the spatial error of FD-
FP method converges fast. For $k\geq 2$, the Bermúdez function
$\sigma_{B_{k}}$ is able to offer a near spectral accuracy for the Fourier
pseudo-spectral method, although the error is a little bit larger than that
for the polynomial choice (2.13) (see Figure 9). Practically, taking $k=2$ or
$3$ is enough, since increasing $k$ does not give any more improvements.

  * 2)

The proposed GMRES solver (2.11) for FD-FP method works very well with the
help of the preconditioner (2.12). The number of iterations needed to reach a
threshold has been significantly reduced to $O(1)$ for all mesh sizes. This
makes the practical efficiency of the FD-FP method comparable to the EWI-FP
method.

  * 3)

The PML error of PML-II with Bermúdez function (2.15) is much smaller than
that of the polynomial choice (2.13) under the same parameters $\sigma_{0}$
and $\delta$ (cf. Figures 10 and 7). Moreover, the accuracy of the PML-II with
Bermúdez function is much less sensitive to the choice of the parameters
$\sigma_{0}$ and $\delta$ (see Figure 10). Therefore, one does not need to
tune the parameters for the layer in a practical computation, which makes it
suitable for concrete applications and to later investigate the non-
relativistic regime of the NKGE.

Energy decay and comparison. To end, we test and compare the energy behaviour
of the two PML formulations. We define the following energy functional for
some function $w=w(x,t)$

$H_{I}(t;w):=\int_{I}\left[|\partial_{t}w(x,t)|^{2}+|\partial_{x}w(x,t)|^{2}+|w(x,t)|^{2}+\frac{\lambda}{2}|w(x,t)|^{4}\right]dx,\quad
t\geq 0.$

With the exact solution $w=u$ of (2.1), the above $H_{I}(t;u)$ denotes the
part of the total energy (1.2) inside the physical domain $I=(-L,L)$ for the
NKGE (2.1). Note that the quantity $H_{I}(t;u)$ is not conserved by (2.1), and
in fact it should be decaying with respect to time, since the waves in the
solution $u(x,t)$ keep propagating to the far field as time evolves.

We take the example (2.16) with $I=(-4,4)$. For PML-I given by (2.3), we
consider $\sigma=\sigma_{P}$ with $\sigma_{0}=8,\,\delta=6/8,\,\alpha=0$. For
PML-II as defined by (2.8), we take $\sigma=\sigma_{B_{2}}$ with
$\sigma_{0}=3,\,\delta=6/8,\,R=1$. We compute $H_{I}$ with the exact solution
$w=u$ of (2.1), and with the exact solution $w=u_{\textrm{pml}}$ of PML-I
(2.3) or PML-II (2.8). The “exact” solutions here are obtained numerically as
before. The numerical results are reported in Figure 11 until the final time
$t=22$. We see that both PML-I and PML-II follow the exact energy decay very
well in long time, and the approximation from PML-II is however more accurate
than PML-I for all times.

$\begin{array}[]{c}\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 11. Energy decay of the PMLs in physical domain: $H_{I}(t)$ obtained
from PML-I, PML-II and the NKGE.

Overall, based on the presented numerical experiments in this section for
approximating the NKGE (2.1) in the classical scaling, we suggest to rather
use the PML-II formulation (2.8) with the regularized Bermúdez absorption
function (2.15) for $k\geq 2$, not only from the point of view of the temporal
and spatial accuracy, but also for the stability to fix the tuning parameters
of the PML.

##  3\. PML for non-relativistic scaling

Physically when $c\to\infty$, the NKGE (1.1) is said to be in the non-
relativistic limit regime [33, 34, 42]. In such a case, it is often convenient
to introduce the scaling [13, 35]:

$x\to\frac{x}{x_{s}},\quad t\to\frac{t}{t_{s}},\quad\lambda\to\lambda
mx_{s}^{2},\quad\mbox{with}\quad
t_{s}=\frac{mx_{s}^{2}}{\hbar},\quad\varepsilon=\frac{\hbar}{mcx_{s}},$

where $t_{s}$ and $x_{s}$ are the time and length units, respectively, and
$\varepsilon>0$ is a dimensionless parameter. Then, the formulation of the
one-dimensional real-valued case of (1.1) in the non-relativistic scaling
reads [12, 15, 23]:

$\left\\{\begin{split}&\varepsilon^{2}\partial_{tt}u(x,t)-\partial_{xx}u(x,t)+\frac{1}{\varepsilon^{2}}u(x,t)+\lambda
u(x,t)^{3}=0,\quad t>0,\ x\in{\mathbb{R}},\\\
&u(x,0)=u_{0}(x),\quad\partial_{t}u(x,0)=v_{0}(x),\quad
x\in{\mathbb{R}}.\end{split}\right.$ (3.1)

In the following, we shall consider $\varepsilon\in(0,1]$. The NKGE (3.1) in
such scaling could describe the physical system where the wave speed is
smaller than the speed of light but with non-negligible relativistic effects
[15, 35, 42]. It also appears in the high-plasma-frequency limit regime for
plasma physics [14, 16]. For $\varepsilon=1,$ (3.1) leads to the classical
NKGE (2.1), and as $\varepsilon\to 0$, we are in the non-relativistic limit
regime.

###  3.1. PML-II

We adopt the same notations as introduced in the previous section. Now in the
same spirit as (2.8), the PML-II in the non-relativistic scaling reads:

$\left\\{\begin{split}&\varepsilon^{2}\partial_{tt}u-\frac{1}{1+R\sigma}\partial_{x}\left(\frac{1}{1+R\sigma}\partial_{x}u\right)+\frac{1}{\varepsilon^{2}}u+\lambda
u^{3}=0,\quad t>0,\ x\in I^{*},\\\
&u(x,0)=u_{0}(x),\quad\partial_{t}u(x,0)=v_{0}(x),\quad x\in I^{*},\\\
&u(-L^{*},t)=u(L^{*},t),\quad t\geq 0.\end{split}\right.$ (3.2)

For the linear case of (3.2), i.e. $\lambda=0$, let us now consider a plane
wave solution: $u(x,t)=\mathrm{e}^{i(kx-\omega t)}$ with
$k,\omega\in{\mathbb{R}}$ satisfying the dispersion relation:

$\omega=\pm\frac{\sqrt{1+\varepsilon^{2}k^{2}}}{\varepsilon^{2}},$

which provides the phase velocity

$v_{p}=\frac{\omega}{k}=\pm\frac{\sqrt{k^{-2}+\varepsilon^{2}}}{\varepsilon^{2}}.$

From the above, it is clear to see that waves with all wavelengths are
travelling to infinity at the speed $O(\varepsilon^{-2})$ as $\varepsilon\to
0$. Therefore, when $\varepsilon$ is small, the waves in (3.1) enter and pass
through the damping layer $L\leq|x|\leq L+\delta$ very quickly, and then get
reflected at the opposite outer boundary $|x|=L+\delta$ because of the
periodic boundary conditions. This fact may leave the layer not enough time to
effectively absorb all the waves. In fact, we will show later by numerical
tests, for a fixed PML setup under the classical polynomial absorption
function (2.13), i.e. fixed $\sigma,\,\delta,\,R$ as $\varepsilon$ decreases,
that the PML-II (3.2) will soon fail to approximate the exact solution of NKGE
in the physical domain.

In order to get a stable absorber with rather uniform damping effect for
$\varepsilon\in(0,1]$, we will verify later by numerical experiments that we
need to choose

$R=R^{\varepsilon}=O(1/\varepsilon^{2})\in{\mathbb{R}}^{+},\quad\mbox{if}\quad\sigma(x)=\sigma_{P}(x),$
(3.3)

for the PML-II (3.2) with fixed $\sigma_{P}(x)$. This is clearly equivalent to
say that with $R>0$ fixed as $\varepsilon\to 0$, we need to choose the
strength of the absorption (2.13) $\sigma_{0}=O(\varepsilon^{-2})$. Based on
the studies from the previous section, we know that the accuracy of the
Bermúdez absorption function (2.15) is not very sensitive to the choice of the
strength parameter $\sigma_{0}$. Therefore, it is hopeful and reasonable that
with the Bermúdez’s absorption function in the PML-II (3.2), the damping
effect is rather uniform for $\varepsilon\in(0,1]$ by using

$R=O(1)\in{\mathbb{R}}^{+},\quad\mbox{if}\quad\sigma(x)=\sigma_{B_{k}}(x).$

It is well-known that as $\varepsilon\to 0$, the solution of the NKGE (3.1)
contains rapid oscillations in time [15, 33, 34, 39]:

$u(x,t)=\mathrm{e}^{it/\varepsilon^{2}}z(x,t)+\mathrm{e}^{-it/\varepsilon^{2}}\overline{z}(x,t)+O(\varepsilon^{2}),\quad
t\geq 0,\ x\in{\mathbb{R}},$ (3.4)

where $z=z(x,t)$ solves a nonlinear Schrödinger equation independent of
$\varepsilon$. In space, the solution $u$ of NKGE (3.1) is rather smooth with
wavelength independent of $\varepsilon$. As an illustrative example, we take

$u_{0}=5\mathrm{e}^{-x^{2}},\quad
v_{0}=\frac{1}{2}\textrm{sech}(x^{2}),\quad\lambda=\frac{1}{2},$ (3.5)

and we show the profile of the exact solution $u(x,t)$ at $t=4$ in the
physical domain $x\in(-4.5,-4)$ under different $\varepsilon$ in Figure 12.
While in the PML-II formulation (3.2), due to the fast reflected waves (or the
incoming waves due to the periodic boundary condition), the
$\varepsilon$-dependent oscillations will be induced to the space variable $x$
of the solution. To see this, we consider in (3.2) the physical domain
$I=(-4,4)$ with the layer setup
$\sigma_{0}=3,\,\delta=0.5,\,R=\varepsilon^{-2}$ in $\sigma_{P}$ (2.13), and
we show in Figure 12 the PML solution $u_{\textrm{pml}}(x,t)$ of (3.2) at
$t=4$ inside the layer $-L-\delta\leq x\leq-L$. The PML solution
$u_{\textrm{pml}}(x,t)$ under the same parameters values but with $R=1$ and
$\sigma_{B_{2}}$ (2.15) is also given in Figure 12. It can be seen from Figure
12 that, as $\varepsilon$ decreases, the PML solution $u_{\textrm{pml}}(x,t)$
exhibits more and more spatial oscillations inside the layer. Such phenomenon
will cause numerical difficulties to the spatial discretization of the PML-II
formulation (3.2). As we shall report in the numerical tests, the accuracy of
the Fourier pseudo-spectral discretization for (3.2) depends on $\varepsilon.$

$\begin{array}[]{c}\psfig{height=113.81102pt,width=170.71652pt}\psfig{height=113.81102pt,width=170.71652pt}\\\
\psfig{height=113.81102pt,width=170.71652pt}\end{array}$ Figure 12. Profiles
of the exact solution $u(x,t)$ of NKGE (top left) and the PML-II solution
$u_{\textrm{pml}}(x,t)$ with $R=\varepsilon^{-2}$ and $\sigma=\sigma_{P}$ (top
right), or with $R=1$ and $\sigma=\sigma_{B_{2}}$ (bottom).

For the FD-FP method, the finite-difference discretization of (3.2) with
$R=R^{\varepsilon}$ directly reads:

$\varepsilon^{2}\frac{u^{n+1}-2u^{n}+u^{n-1}}{\tau^{2}}+\frac{A^{\varepsilon}}{2}(u^{n+1}+u^{n-1})+\frac{1}{2\varepsilon^{2}}(u^{n+1}+u^{n-1})+\lambda(u^{n})^{3}=0,$

which gives the FD-FP scheme

$\displaystyle
u^{n+1}=-u^{n-1}+(G^{\varepsilon})^{-1}\left[\frac{2\varepsilon^{2}}{\tau^{2}}u^{n}-\lambda(u^{n})^{3}\right],\quad
n\geq 1,$ (3.6)

with

$\displaystyle
G^{\varepsilon}=\left[\left(\frac{\varepsilon^{2}}{\tau^{2}}+\frac{1}{2\varepsilon^{2}}\right)\mathbb{I}+\frac{1}{2}A^{\varepsilon}\right],\quad
A^{\varepsilon}=-d_{0}D_{1}d_{0}D_{1},\quad
d_{0}=\mathrm{diag}(1/(1+R^{\varepsilon}\sigma)).$

For the above matrix $G^{\varepsilon}$ that depends on $\varepsilon$ in the
scheme, we will show later by numerical results that its condition number
remains small and stays uniformly bounded as $\varepsilon\to 0$. Similarly as
(2.10), it can also be implemented efficiently by the GMRES solver with a
preconditioner in the non-relativistic scaling:

$\left\\{\begin{array}[]{l}\displaystyle u^{n+1}=-u^{n-1}+w^{n},\quad n\geq
1,\\\
\displaystyle\mathcal{P}^{\varepsilon}G^{\varepsilon}w^{n}=\mathcal{P}^{\varepsilon}\left[\frac{2\varepsilon^{2}}{\tau^{2}}u^{n}-\lambda(u^{n})^{3}\right],\end{array}\right.\quad\mathcal{P}^{\varepsilon}=\left(\frac{\varepsilon^{2}}{\tau^{2}}+\frac{1}{2\varepsilon^{2}}-\frac{\partial_{xx}}{2}\right)^{-1}.$
(3.7)

The starting value for (3.6) if directly using the construction from the
Taylor expansion reads

$u^{1}=u_{0}+\tau
v_{0}-\frac{\tau^{2}}{2\varepsilon^{2}}\left[A^{\varepsilon}u_{0}+\frac{u_{0}}{\varepsilon^{2}}+\lambda
u_{0}^{3}\right],$

where the $O(1/\varepsilon^{2})$ and $O(1/\varepsilon^{4})$ terms in the above
will induce instability to the scheme, since $u^{1}=O(1/\varepsilon^{4})$
could be very large for small $\varepsilon$. To avoid this, we consider a
filtered data:

$u^{1}=u_{0}+\tau
v_{0}-\frac{\tau}{2}\sin(\tau/\varepsilon^{2})\left[A^{\varepsilon}u_{0}+\lambda
u_{0}^{3}\right]-\frac{\tau}{2}\sin(\tau/\varepsilon^{4})u_{0},$ (3.8)

so that $u^{1}=O(1)$.

Thanks to the filtered starting value (3.8) as well as the time averaging, the
FD-FP method (3.6) is stable for all $\varepsilon\in(0,1]$. This means that we
have the boundedness of the numerical solution of FD-FP up to a fixed time if
$\tau\leq C$ for some $C>0$ independent of $\varepsilon$ or the spatial mesh
size $h$. However, owning to the $\varepsilon$-dependent oscillation
frequencies in both time and space directions of the PML solution of (3.2), it
is expected that the FD-FP method (2.9) is not uniformly accurate for
approximating (3.2) as $\varepsilon\to 0$. One then needs to use mesh sizes
$\tau>0$ and $h>0$ that depend on $\varepsilon$ to reach the same accuracy
level at a fixed time.

######  Remark 3.1.

It is beyond the scope of this paper to investigate whether the asymptotic
expansion (3.4) also holds for the PML-II equation (3.2) or not. If it is
still valid, then formally as $\varepsilon\to 0$, the PML-II (3.2) is
consistent with the limit: a nonlinear Schrödinger equation with PML presented
in [9, 43]. We shall address this question in a future study.

###  3.2. PML-I

Next, we consider the PML-I formulation under the non-relativistic scaling. To
write the correct $\varepsilon$-dependent PML model for (3.1), we shall
proceed the derivation as in [2].

We begin by considering a linear case of (3.1):

$\varepsilon^{2}\partial_{tt}u-\partial_{xx}u+\frac{1}{\varepsilon^{2}}u+\lambda
u=0.$ (3.9)

By taking the Laplace transform on both sides of the above equation, we get

$\left(\varepsilon^{2}s^{2}-\partial_{xx}+\frac{1}{\varepsilon^{2}}+\lambda\right)\widehat{u}(x,s)=0,\quad\widehat{u}(x,s)=\int_{0}^{\infty}u(x,t)\mathrm{e}^{-st}dt,$

where $s$ is the Laplace variable. By (2.2), we find formally, for some
$\alpha\geq 0$ and absorption function $\sigma(x)$,

$\left(\varepsilon^{2}s^{2}+\frac{1}{\varepsilon^{2}}+\lambda\right)\widehat{u}=\frac{s+\alpha}{s+\alpha+\sigma}\partial_{x}\left(1-\frac{\sigma}{s+\alpha+\sigma}\right)\partial_{x}\widehat{u}.$
(3.10)

Let

$\widehat{\eta}_{1}=-\frac{1}{s+\alpha+\sigma}\partial_{x}\widehat{u},$

then from (3.10) we have

$\displaystyle\varepsilon^{2}\left(s^{2}+s\sigma\right)\widehat{u}-\varepsilon^{2}\sigma\alpha\widehat{u}+\sigma\frac{\varepsilon^{2}\alpha^{2}}{s+\alpha}\widehat{u}=\partial_{x}\left(\partial_{x}\widehat{u}+\sigma\widehat{\eta}_{1}\right)-\left(\frac{1}{\varepsilon^{2}}+\lambda\right)\widehat{u}-\sigma\frac{\varepsilon^{-2}+\lambda}{s+\alpha}\widehat{u}.$

By taking the inverse Laplace transform of the above equation and further
letting

$\widehat{\eta}_{2}=-\frac{\varepsilon^{2}\alpha^{2}+\varepsilon^{-2}+\lambda}{s+\alpha}\widehat{u},$

we obtain

$\displaystyle\varepsilon^{2}\partial_{tt}u-\partial_{xx}u+\frac{1}{\varepsilon^{2}}u+\lambda
u=\partial_{x}(\sigma\eta_{1})+\varepsilon^{2}\sigma\alpha
u-\varepsilon^{2}\sigma\partial_{t}u+\sigma\eta_{2},$
$\displaystyle\partial_{t}\eta_{1}+(\alpha+\sigma)\eta_{1}+\partial_{x}u=0,\quad\partial_{t}\eta_{2}+\alpha\eta_{2}+\left(\varepsilon^{2}\alpha^{2}+\frac{1}{\varepsilon^{2}}\right)u+\lambda
u=0,$

which is a PML system for the linear model problem (3.9). To proceed from the
linear case to the nonlinear case, as proposed in [2], we directly replace the
term $\lambda u$ by $\lambda u^{3}$ in the above expressions. Then together
with the domain truncation similarly as (2.3), the PML-I formulation for NKGE
(3.1) in the non-relativistic scaling reads

$\left\\{\begin{split}&\varepsilon^{2}\partial_{tt}u-\partial_{xx}u+\frac{1}{\varepsilon^{2}}u+\lambda
u^{3}=\sigma\left(\eta_{2}+\varepsilon^{2}\alpha
u-\varepsilon^{2}\partial_{t}u\right)+\partial_{x}(\sigma\eta_{1}),\quad t>0,\
x\in I^{*},\\\
&\partial_{t}\eta_{1}+(\alpha+\sigma)\eta_{1}+\partial_{x}u=0,\\\
&\partial_{t}\eta_{2}+\alpha\eta_{2}+\left(\varepsilon^{2}\alpha^{2}+\frac{1}{\varepsilon^{2}}\right)u+\lambda
u^{3}=0,\\\
&u(x,0)=u_{0}(x),\quad\partial_{t}u(x,0)=v_{0}(x),\quad\eta_{1}(x,0)=\eta_{2}(x,0)=0,\quad
x\in I^{*},\\\
&u(-L^{*},t)=u(L^{*},t),\quad\eta_{1}(-L^{*},t)=\eta_{1}(L^{*},t),\quad\eta_{2}(-L^{*},t)=\eta_{2}(L^{*},t),\quad
t\geq 0.\end{split}\right.$ (3.11)

Asymptotically, the limit of the PML-I system (3.11) as $\varepsilon\to 0$ is
not clear. It can be seen from (3.11) that there is a $O(\varepsilon^{-2})$
term in the equation of $\eta_{2}$, which acts as a very stiff source term as
$0<\varepsilon\ll 1$. Numerically, to propose a stable and accurate numerical
algorithm for solving (3.11) is a challenge for small $\varepsilon$. In this
paper, we do not intend to address this issue, and we shall focus on the
investigation of the two PML formulations and the two choices of absorption
functions. In the following, let us simply apply the EWI-FP discretization
from Section 2.1.

The EWI-FP method for (3.11) then reads: for $n\geq 0,$

$\displaystyle u^{n+1}=$
$\displaystyle\cos(\langle\partial_{x}\rangle_{\varepsilon}\tau)u^{n}+\frac{\sin(\langle\partial_{x}\rangle_{\varepsilon}\tau)}{\langle\partial_{x}\rangle_{\varepsilon}}v^{n}+\frac{\tau\sin(\langle\partial_{x}\rangle_{\varepsilon}\tau)}{2\varepsilon^{2}\langle\partial_{x}\rangle_{\varepsilon}}f^{n},$
(3.12a) $\displaystyle\eta_{1}^{n+1}=$
$\displaystyle\mathrm{e}^{-(\sigma+\alpha)\tau}\eta_{1}^{n}-\frac{\tau}{2}\left[\mathrm{e}^{-(\sigma+\alpha)\tau}\partial_{x}u^{n}+\partial_{x}u^{n+1}\right],$
(3.12b) $\displaystyle\eta_{2}^{n+1}=$
$\displaystyle\mathrm{e}^{-\alpha\tau}\eta_{2}^{n}-\frac{\tau}{2}\left[\mathrm{e}^{-\alpha\tau}\left((\varepsilon^{2}\alpha^{2}+\varepsilon^{-2})u^{n}+\lambda(u^{n})^{3}\right)+(\varepsilon^{2}\alpha^{2}+\varepsilon^{-2})u^{n+1}+\lambda(u^{n+1})^{3}\right],$
(3.12c) $\displaystyle v^{n+1}=$
$\displaystyle-\langle\partial_{x}\rangle_{\varepsilon}\sin(\langle\partial_{x}\rangle_{\varepsilon}\tau)u^{n}+\cos(\langle\partial_{x}\rangle_{\varepsilon}\tau)v^{n}+\frac{\tau}{2\varepsilon^{2}}\left[\cos(\langle\partial_{x}\rangle_{\varepsilon}\tau)f^{n}+f^{n+1}\right],$
(3.12d)

with

$\langle\partial_{x}\rangle_{\varepsilon}=\frac{\sqrt{1-\varepsilon^{2}\partial_{xx}}}{\varepsilon^{2}},\quad
f^{n}=\sigma\left(\eta_{2}^{n}-\varepsilon^{2}v^{n}+\varepsilon^{2}\alpha
u^{n}\right)+\partial_{x}(\sigma\eta_{1}^{n})-\lambda(u^{n})^{3}.$

As we discussed in the previous section, for each fixed $\varepsilon$, the
EWI-FP has a CFL as $\tau<Ch$. Now, the scheme is much more constrained since
for each fixed $h$ as $\varepsilon\to 0$, the EWI-FP scheme (3.12) has a
severe stability constraint: $\tau\lesssim\varepsilon^{2}$, due to the
stiffness that we mentioned above. Moreover, we are not able to directly apply
the Bermúdez absorption function in (3.12) due to the singularity in
$\sigma_{B_{k}}$ (2.15) and stability issue. So, we have to consider the
classical choice (2.13) for (3.11) or (3.12) and tune the parameters for every
$\varepsilon$. Concerning these drawbacks, we recommend to rather use the PML-
II formulation (3.2) than PML-I (3.11) for NKGE in the non-relativistic
scaling.

###  3.3. Numerical results

In this subsection, we present the numerical results for the above two PML
formulations for the NKGE in the non-relativistic scaling. Let us take the
given data in (3.5) for the NKGE (3.1) in the following tests. We begin with
the numerical investigations of PML-II, and we shall consider similarly as
before its performance under the two kinds of absorption functions, and
finally analyze shortly PML-I.

PML-II with polynomial type absorption function. We start with the polynomial
type absorption function $\sigma(x)=\sigma_{P}(x)$ as (2.13) for the PML-II
(3.2). First of all, we illustrate the performance of the PML-II formulation
(3.2) with (2.13) under different choices of the parameter $R>0$. We measure
the error (2.18) on the physical domain $I=(-4,4)$ by choosing

$R=1,\quad\mbox{or}\quad R=1/\varepsilon,\quad\mbox{or}\quad
R=1/\varepsilon^{2},$

with $\sigma_{0}=3$ fixed and some values of the layer thickness $\delta$. We
show in Figure 13 the error (2.18) of the PML-II (3.2) for approximating (3.1)
at $t=4$ under different $\varepsilon$, where the last choice is clearly the
most accurate.

$\begin{array}[]{c}\psfig{height=113.81102pt,width=170.71652pt}\psfig{height=113.81102pt,width=170.71652pt}\\\
\psfig{height=113.81102pt,width=170.71652pt}\end{array}$ Figure 13. PML-II
$e^{2}_{\textrm{pml}}$-error with respect to $\varepsilon$ for fixed
$\sigma_{0}=3$ and different parameter values $\delta$: $\sigma=\sigma_{P}$
with strategy $R=1$ (top left), $R=\varepsilon^{-1}$ (top right) and
$R=\varepsilon^{-2}$ (bottom).

Next, we test the performance of the FD-FP method (3.6) for approximating the
PML-II formulation (3.2) for different $\varepsilon$ under the choice
$R=R^{\varepsilon}=1/\varepsilon^{2}$. In these tests, we fix $\sigma_{0}=3$
and $\delta=4/8$ for the layer. The reference solution is obtained by using
the FD-FP method with a very fine mesh, e.g. $\tau=10^{-4},\,h=1/2048$ for
each $\varepsilon$. The temporal and spatial errors (2.17) of the FD-FP method
under different $\varepsilon$ at $t=4$ are shown in Figure 14. Moreover, by
taking $\tau=10^{-4},\,h=1/128$ as mesh sizes in this test, we show in Figure
15 the condition number (in spectral norm) of the iterative matrix
$G^{\varepsilon}$ as a function of $\varepsilon$ for the FD-FP method (3.6).

$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 14. Temporal (left) and spatial (right) accuracy of FD-FP:
$e^{n,\infty}_{\textrm{pml}}$-error for $\sigma=\sigma_{P}$ with different
$\varepsilon$.
$\begin{array}[]{c}\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 15. The condition number of the iterative matrix $G^{\varepsilon}$ of
FD-FP for PML-II: polynomial case $\sigma=\sigma_{P}$ and Bermúdez type
$\sigma=\sigma_{B_{2}}$.

Based on the numerical results in Figures 13-15, we can conclude on the
following observations for the PML-II formulation and the FD-FP method with
polynomial type absorption function $\sigma_{P}$:

  * 1)

With the strategy $R=1/\varepsilon^{2}$, the accuracy of the PML-II is uniform
with respect to $\varepsilon$ as far as we could test. In contrast, if one
takes $R=1$ or $R=1/\varepsilon$, the error of PML-II increases dramatically
as $\varepsilon$ decreases (see Figure 13). Overall, when $\varepsilon$ is
small, the accuracy of PML-II with $\sigma_{P}$ is much more sensitive than
for the classical scaling case to the parameters $R$ or saying equivalently
the strength $R\sigma_{0}$ and the layer size $\delta$.

  * 2)

For $0<\varepsilon\leq 1$, the condition number of the iterative matrix
$G^{\varepsilon}$ in the FD-FP scheme (3.6) stays uniformly bounded and small
(see Figure 15), so the linear system in (3.6) can be solved accurately and
efficiently. For each fixed $\varepsilon$, we still have the second-order
temporal accuracy and the near spectral accuracy in space for the FD-FP method
(see Figure 14). When $\varepsilon$ decreases with fixed mesh size $\tau$ and
$h$, the approximation error of FD-FP in both time and space increases
quickly. The $\varepsilon$-dependence of the temporal accuracy is due to the
$O(\varepsilon^{-2})$ oscillation frequency in time, and the spatial accuracy
is affected by the fast reflections at the boundary of the layer as we
explained before. A future goal is therefore to construct a numerical
algorithm that is uniformly accurate for $0<\varepsilon\leq 1$, which would
certainly make the PML more efficient.

PML-II with Bermúdez type absorption function. Let us consider
$\sigma(x)=\sigma_{B_{k}}(x)$ as defined in (2.15) for the PML-II (3.2). To
test the PML error, we take the PML-II (3.2) with $\sigma=\sigma_{B_{2}}$,
i.e. $k=2$ in (2.15), for $\sigma_{0}=3$ fixed. By using some different values
of $R>0$ and $\delta>0$, we show in Figure 16 the PML approximation error
(2.18) on the physical domain $I=(-4,4)$ at $t=4$. With
$\sigma=\sigma_{B_{2}},\,R=1,\,\delta=4/8$, we plot in Figure 17 the numerical
discretization error (2.17) of the FD-FP method at $t=4$ in time and in space.
Under $\tau=10^{-4},\,h=1/128$, the corresponding condition number of the
matrix $G^{\varepsilon}$ as a function of $\varepsilon$ is also shown in
Figure 15. The number of iterations needed for the convergence of GMRES (for
threshold $\epsilon=10^{-10}$ without restart) with the preconditioner (3.7)
for solving (3.6) at $n=1$ under $\tau=0.02,\ \delta=3/8,\ R=1$ is given in
Figure 18 (left) for different $\varepsilon>0$ and values of $h$. The number
of iterations needed for solving (3.6) under $\tau=2\times
10^{-4},\,h=1/128,\,\delta=3/8,\ R=1$ until time $t=6$ is given in Figure 18
(right). Let us remark that, without this preconditioner, the iterative solver
may not converge when $\varepsilon$ gets smaller.

$\begin{array}[]{c}\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 16. PML-II $e^{2}_{\textrm{pml}}$-error with respect to $\varepsilon$
for fixed $\sigma_{0}=3$ and different $\delta$ or $R$ for
$\sigma=\sigma_{B_{2}}$.
$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 17. Temporal (left) and spatial (right) accuracy of FD-FP:
$e^{n,\infty}_{\textrm{pml}}$-error for $\sigma=\sigma_{B_{2}}$ with different
$\varepsilon$.
$\begin{array}[]{cc}\psfig{height=113.81102pt,width=170.71652pt}&\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 18. Number of iterations needed by GMRES (fixed threshold
$\epsilon=10^{-10}$) with the preconditioner for (3.6) at $n=1$ (left) and for
(3.6) until $t=6$ (right).

It can be seen from Figures 16-18 that

  * 1)

The numerical discretization error of the FD-FP method for PML-II (3.2) with
the Bermúdez absorption function (2.15) is close to the case of the polynomial
choice above (cf. Figures 14 and 17). The temporal and spatial error grows
dramatically as $\varepsilon$ decreases, but with $\tau,\,h\to 0$, the FD-FP
method still converges and the condition number of $G^{\varepsilon}$ is again
uniformly bounded for $0<\varepsilon\leq 1$. The preconditioner (3.7) still
works very well for $0<\varepsilon<1$, which improves significantly the
efficiency of the GMRES solver, since it exhibits a convergence rate
independent of the mesh refinement $h$. Moreover, the number of iterations
stays small and bounded in time, and is relatively insensitive to
$\varepsilon$.

  * 2)

The PML-II (3.2) with the Bermúdez absorption function (2.15) is much more
accurate than that of the classical polynomial choice (2.13) (cf. Figures 16
and 13). With the Bermúdez function, the PML-II is still not sensitive to the
parameters for the layer, i.e. $R,\,\sigma_{0},\,\delta$. More importantly,
the error appears to be rather uniform with respect to $\varepsilon\in(0,1]$
as far as we could test. Thus, we can simply take fixed
$R,\,\sigma_{0},\,\delta=O(1)$ for all $\varepsilon\in(0,1]$. While, we remark
that for smaller $\varepsilon$, it is getting harder to eliminate the impact
from the numerical error of the FD-FP solver due to error behaviour observed
in Figure 17.

Tests for PML-I. At last, we report some numerical results for the PML-I
formulation (3.11) when approximating the NKGE (3.1) and illustrate its
drawbacks compared with PML-II (3.2). We use the same numerical example (3.5)
and take $\alpha=0$ in (3.11). The PML error (2.18) of the PML-I formulation
(3.11) at $t=4$ with respect to $\varepsilon$ is given in Figure 19 under
$\sigma_{0}=8$ and different parameter values $\delta$. It is clear from the
numerical results in Figure 19 that the PML-I formulation (3.11) with the
classical absorption function (2.13) is very sensitive to the parameters:
$\sigma_{0},\,\delta$ and $\varepsilon$. When $\varepsilon$ becomes small, the
error of PML-I is much larger than that of the PML-II (cf. Figure 19 and
Figures 13 and 16). In particular, for small values of $\varepsilon$, the
accuracy of PML-I can barely be improved by increasing $\delta$. One must keep
enlarging the strength parameter $\sigma_{0}$, while the stiffness brought by
large $\sigma$ in (3.11) will consequently increase the difficulty for the
numerical solver.

$\begin{array}[]{c}\psfig{height=113.81102pt,width=170.71652pt}\psfig{height=113.81102pt,width=170.71652pt}\\\
\psfig{height=113.81102pt,width=170.71652pt}\end{array}$ Figure 19. PML-I
$e^{2}_{\textrm{pml}}$-error with respect to $\varepsilon$ for different
$\delta$: $\sigma_{0}=8$ (top left), $\sigma_{0}=8/\varepsilon$ (top right)
and $\sigma_{0}=8/\varepsilon^{2}$ (bottom).

##  4\. Extension to the 2D rotating NKGE

In this section, we shall apply the PML-II formulation to simulate the
vortices dynamics in a two-dimensional rotating NKG model. The rotating NKGE
has been introduced to model the dynamics of the cosmic superfluid in a
rotating frame, and in two space dimensions it reads [35, 42]:

$\displaystyle\partial_{tt}\Psi-\Delta\Psi+\Psi+\lambda|\Psi|^{2}\Psi-2i\Omega
L_{z}\partial_{t}\Psi-\Omega^{2}L_{z}^{2}\Psi=0,\quad\mathbf{x}=(x,y)^{\intercal}\in{\mathbb{R}}^{2},\
t>0,$
$\displaystyle\Psi(\mathbf{x},0)=\psi_{0}(\mathbf{x}),\quad\partial_{t}\Psi(\mathbf{x},0)=\psi_{1}(\mathbf{x}),\quad\mathbf{x}\in{\mathbb{R}}^{2},$

where the solution $\Psi=\Psi(\mathbf{x},t)$ is a complex-valued scalar field,
$\Omega>0$ is the given angular velocity and
$L_{z}=-i(x\partial_{y}-y\partial_{x})$ is the angular momentum operator. As
proposed in [35], by introducing the rotating matrix

$A(t):=\begin{pmatrix}\cos(\Omega t)&\sin(\Omega t)\\\ -\sin(\Omega
t)&\cos(\Omega t)\end{pmatrix},\quad t\geq 0,$

and a change of variable to the _rotating Lagrangian coordinates_ :

$u(\tilde{\mathbf{x}},t)=\Psi(\mathbf{x},t),\quad\mathbf{x}=A(t)\tilde{\mathbf{x}},$

the two angular momentum terms can be eliminated and the model becomes the
NKGE (removing the $\tilde{\,}$ for conciseness) [35]:

$\left\\{\begin{split}&\partial_{tt}u(\mathbf{x},t)-\Delta
u(\mathbf{x},t)+u(\mathbf{x},t)+\lambda|u(\mathbf{x},t)|^{2}u(\mathbf{x},t)=0,\quad\mathbf{x}\in{\mathbb{R}}^{2},\
t>0,\\\
&u(\mathbf{x},0)=\psi_{0}(\mathbf{x}),\quad\partial_{t}u(\mathbf{x},0)=\Omega\nabla\psi_{0}(\mathbf{x})\cdot\binom{y}{-x}+\psi_{1}(\mathbf{x}),\quad\mathbf{x}\in{\mathbb{R}}^{2}.\end{split}\right.$
(4.1)

Then, we can directly extend the PML-II formulation to the above initial value
problem which leads to

$\left\\{\begin{split}&\partial_{tt}u(\mathbf{x},t)-\frac{1}{1+R\sigma(x)}\partial_{x}\left(\frac{1}{1+R\sigma(x)}\partial_{x}u(\mathbf{x},t)\right)-\frac{1}{1+R\sigma(y)}\partial_{y}\left(\frac{1}{1+R\sigma(y)}\partial_{y}u(\mathbf{x},t)\right)\\\
&\hskip
85.35826pt+u(\mathbf{x},t)+\lambda|u(\mathbf{x},t)|^{2}u(\mathbf{x},t)=0,\quad
t>0,\ \mathbf{x}\in I^{*}=(-L^{*},L^{*})^{2},\\\
&u(\mathbf{x},0)=\psi_{0}(\mathbf{x}),\quad\partial_{t}u(\mathbf{x},0)=\Omega\nabla\psi_{0}(\mathbf{x})\cdot\binom{y}{-x}+\psi_{1}(\mathbf{x}),\quad\mathbf{x}\in
I^{*},\\\ &u(-L^{*},y,t)=u(L^{*},y,t),\quad t\geq 0,\quad
y\in[-L^{*},L^{*}],\\\ &u(x,-L^{*},t)=u(x,L^{*},t),\quad t\geq 0,\quad
x\in[-L^{*},L^{*}],\end{split}\right.$ (4.2)

where $L^{*}=L+\delta$. With

$S_{x}=-\frac{1}{1+R\sigma(x)}\partial_{x}\left(\frac{1}{1+R\sigma(x)}\partial_{x}\right),\quad
S_{y}=-\frac{1}{1+R\sigma(y)}\partial_{y}\left(\frac{1}{1+R\sigma(y)}\partial_{y}\right),$

for short, the finite-difference time discretization for (4.2) then reads

$\displaystyle\frac{u^{n+1}-2u^{n}+u^{n-1}}{\tau^{2}}+\frac{S_{x}+S_{y}}{2}(u^{n+1}+u^{n-1})+\frac{1}{2}(u^{n+1}+u^{n-1})+\lambda|u^{n}|^{2}u^{n}=0,\quad
n\geq 1,$ $\displaystyle u^{1}=u_{0}+\tau
v_{0}-\frac{\tau^{2}}{2}\left[(S_{x}+S_{y})u_{0}+u_{0}+\lambda|u_{0}|^{2}u_{0}\right],$

where $u^{n}=u^{n}(\mathbf{x})\approx u(\mathbf{x},t)$ for $n\geq 1$ and
$u_{0}=u(\mathbf{x},0),\ v_{0}=\partial_{t}u(\mathbf{x},0)$. By applying the
Fourier pseudo-spectral discretization in the $x$\- and $y$-directions to the
equations, we obtain the FD-FP scheme. By further incorporating the previous
equations with the preconditioned GMRES and FFT in 2D, we can implement the
FD-FP scheme in the same manner as before:

$\left\\{\begin{array}[]{l}\displaystyle u^{n+1}=-u^{n-1}+w^{n},\quad n\geq
1,\\\
\displaystyle\mathcal{P}Gw^{n}=\mathcal{P}\left[\frac{2}{\tau^{2}}u^{n}-\lambda|u^{n}|^{2}u^{n}\right],\end{array}\right.\quad$

where

$G=\left[\left(\frac{1}{\tau^{2}}+\frac{1}{2}\right)\mathbb{I}+\frac{S_{x}+S_{y}}{2}\right],\quad\mathcal{P}=\left(\frac{1}{\tau^{2}}+\frac{1}{2}-\frac{\partial_{xx}}{2}-\frac{\partial_{yy}}{2}\right)^{-1}.$

To illustrate the accuracy of the PML-II approach, we consider the following
example. We take in (4.2) the parameters

$L=4,\quad\lambda=3,\quad\Omega=2,\quad R=1,$

and four initially separated vortices

$\psi_{0}(\mathbf{x})=\psi_{1}(\mathbf{x})=(x-c_{0}+iy)(x+c_{0}+iy)(x+i(y-c_{0}))(x+i(y+c_{0}))\mathrm{e}^{-(x^{2}+y^{2})/2},$

with $c_{0}=1.32$. The profile of the initial data for the NKGE is shown in
Figure 20.

$\begin{array}[]{cc}\psfig{height=128.0374pt,width=128.0374pt}&\psfig{height=128.0374pt,width=128.0374pt}\end{array}$
Figure 20. Contour plot of the initial data: $|\psi_{0}(\mathbf{x})|$ (left)
and $\mathrm{Arg}(\psi_{0}(\mathbf{x}))$ (right) for (4.1) or (4.2).

We use the Bermúdez function $\sigma=\sigma_{B_{2}}$ (2.15) as the absorption
function and choose $\sigma_{0}=3$ and $\delta=0.5$. By applying the FD-FP
similarly as in Section 2.2, we solve the PML equation (4.2) accurately and
look at the solution $u_{\textrm{pml}}$ on the physical domain $I=(-L,L)^{2}$.
As a reference solution, we solve (4.1) directly on the sufficiently large
domain $(-16,16)^{2}$ and compare the exact solution $u$ with the PML
solution. The dynamics of the exact solution $u(\mathbf{x},t)$ and the PML
solution $u_{\textrm{pml}}(\mathbf{x},t)$ are shown on the domain $I$ in
Figure 21. The corresponding trajectory of the decaying energy $H_{I}(t;w)$:

$H_{I}(t;w):=\int_{I}\left[|\partial_{t}w(x,t)|^{2}+|\nabla
w(x,t)|^{2}+|w(x,t)|^{2}+\frac{\lambda}{2}|w(x,t)|^{4}\right]dx,\quad t\geq
0,$

with $w=u$ or $u_{\textrm{pml}}$ is reported in Figure 22. For both the
solution and the energy, we see an excellent agreement between both solutions
that cannot be reached when using standard periodic boundary conditions,
unless taking a huge spatial domain.

Exact solution: $|u(\mathbf{x},t)|$ (top) and $\mathrm{Arg}(u(\mathbf{x},t))$
(bottom)

$\begin{array}[]{ccc}\psfig{height=128.0374pt,width=128.0374pt}&\psfig{height=128.0374pt,width=128.0374pt}&\psfig{height=128.0374pt,width=128.0374pt}\\\
\psfig{height=128.0374pt,width=128.0374pt}&\psfig{height=128.0374pt,width=128.0374pt}&\psfig{height=128.0374pt,width=128.0374pt}\end{array}$

PML solution: $|u_{\textrm{pml}}(\mathbf{x},t)|$ (top) and
$\mathrm{Arg}(u_{\textrm{pml}}(\mathbf{x},t))$ (bottom)

$\begin{array}[]{ccc}\psfig{height=128.0374pt,width=128.0374pt}&\psfig{height=128.0374pt,width=128.0374pt}&\psfig{height=128.0374pt,width=128.0374pt}\\\
\psfig{height=128.0374pt,width=128.0374pt}&\psfig{height=128.0374pt,width=128.0374pt}&\psfig{height=128.0374pt,width=128.0374pt}\end{array}$

Figure 21. Contour plots of exact solution $|u(\mathbf{x},t)|$ and
$\mathrm{Arg}(u(\mathbf{x},t))$ of (4.1) and the PML solution
$|u_{\textrm{pml}}(\mathbf{x},t)|$ and
$\mathrm{Arg}(u_{\textrm{pml}}(\mathbf{x},t))$ of (4.2) at times $t=2,4,6$ on
the domain $I$.
$\begin{array}[]{c}\psfig{height=113.81102pt,width=170.71652pt}\end{array}$
Figure 22. The energy decay of the 2D example: $H_{I}(t)$ with respect to
time.

##  5\. Conclusion

In this work, we investigated numerically two kinds of PML formulations for
approximating the nonlinear Klein-Gordon equation: a first-order formulation
borrowed from the study for general wave equations (PML-I), and a second-order
formulation proposed as analogy for the Helmholtz and nonlinear Schrödinger
equations (PML-II). For both PML formulations, we developed efficient and
accurate pseudo-spectral schemes for the numerical discretizations: i) a
pseudo-spectral exponential integrator scheme for PML-I and ii) a linearly
implicit preconditioned pseudo-spectral finite-difference scheme for PML-II.
To obtain a spectral accuracy, the PML absorption functions were locally
regularized. For the NKGE in classical scaling, the PML-II with smoothed
Bermúdez type absorption functions and the resulting pseudo-spectral finite-
difference scheme are proved to be the most accurate and efficient, owning to
the observation that it is less subject to the parameter tuning problems for
the layer. Such advantage also distinguishes PML-II for approximating NKGE in
the non-relativistic scaling, where the PML solution is found to be not
sensitive to the small parameter $\varepsilon$ arising in the non-relativistic
limit regime. We also extend the method to the accurate computation of vortex
dynamics for a rotating NKGE.

## Acknowledgements

Research conducted within the context of the Sino-French International
Associated Laboratory for Applied Mathematics - LIASFMA. X. Antoine thanks the
LIASFMA funding support of the Université de Lorraine. X. Zhao is partially
supported by the Natural Science Foundation of Hubei Province No. 2019CFA007
and the NSFC 11901440.

## References

  * [1] D. Appelö, T. Hagstrom, G. Kreiss, Perfectly matched layers for hyperbolic systems: general formulation, well-posedness, and stability, SIAM J. Appl. Math. 67 (2006), pp. 1-23.
  * [2] D. Appelö, G. Kreiss, Application of a perfectly matched layer to the nonlinear wave equation, Wave Motion 44 (2007), pp. 531-548.
  * [3] X. Antoine, A. Arnold, C. 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), pp.729-796.
  * [4] X. Antoine, W. Bao and C. Besse, Computational methods for the dynamics of nonlinear Schrödinger and Gross-Pitaevskii equations, (A Feature Article) Comput. Phys. Commun. 184 (2013), pp. 2621-2633.
  * [5] X. Antoine, R. Duboscq, Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates, J. Comput. Phys. 258 (2014), pp.509-523.
  * [6] X. Antoine, R. Duboscq, GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: computation of stationary solutions, Comput. Phys. Commun. 185 (2014), pp. 2969-2991.
  * [7] X. Antoine, R. Duboscq, GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: dynamics and stochastic simulations, Comput. Phys. Commun. 193 (2015), pp. 95-117.
  * [8] X. Antoine and R. Duboscq, Modeling and computation of Bose-Einstein Condensates: stationary states, nucleation, dynamics, stochasticity, in Nonlinear Optical and Atomic Systems: at the Interface of Mathematics and Physics, Lecture Notes in Mathematics, 2146, pp. 49-145, Springer, 2015.
  * [9] X. Antoine, E. Lorin, Q. Tang, A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations, Molecular Physics 115 (2017), pp. 1861-1879.
  * [10] X. Antoine, C. Geuzaine, Q. Tang, Perfectly matched layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates, Commun. Nonlinear Sci. Numer. Simulat. 90 (2020), pp. 105406.
  * [11] W. Bao, Y. Cai, Mathematical theory and numerical methods for Bose-Einstein condensation, Kinet. Relat. Mod. 6 (2013), pp. 1-135.
  * [12] W. Bao, Y. Cai, X. Zhao, A uniformly accurate multiscale time integrator pseudospectral method for the Klein–Gordon equation in the nonrelativistic limit regime, SIAM J. Numer. Anal. 52 (2014), pp. 2488-2511.
  * [13] W. Bao, X. Dong, Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime, Numer. Math. 120 (2012), pp. 189-229.
  * [14] W. Bao, X. Dong, X. Zhao, An exponential wave integrator sine pseudospectral method for the Klein–Gordon–Zakharov system, SIAM J. Sci. Comput. 35 (2013), A2903-A2927.
  * [15] W. Bao, X. Zhao, Comparison of numerical methods for the nonlinear Klein-Gordon equation in the nonrelativistic limit regime, J. Comput. Phys. 398 (2019), pp. 108886.
  * [16] W. Bao, X. Zhao, A uniformly accurate multiscale time integrator spectral method for the Klein-Gordon-Zakharov system in the high-plasma-frequency limit regime, J. Comput. Phys. 327 (2016), pp. 270-293.
  * [17] G. Bao, H. Wu, Convergence analysis of the perfectly matched layer problems for time-harmonic Maxwell’s equations, SIAM J. Numer. Anal. 43 (2005), pp. 2121–2143.
  * [18] J.P. Bérenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 114 (1994), pp. 185-200.
  * [19] A. Bermúdez, L. Hervella-Nieto, A. Prieto, R. Rodríguez, An exact bounded PML for the Helmholtz equation, C. R. Acad. Sci. Paris, Ser. I 339 (2004), pp. 803-808.
  * [20] A. Bermúdez, L. Hervella-Nieto, A. Prieto, R. Rodríguez, An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems, J. Comput. Phys. 223 (2007), pp. 469-488.
  * [21] A. Bermúdez, L. Hervella-Nieto, A. Prieto, R. Rodríguez, An exact bounded perfectly matched layer for time-harmonic scattering problems, SIAM J. Sci. Comput. 30 (2007), pp. 312-338.
  * [22] D.D. Băinov, E. Minchev, Nonexistence of global solutions of the initial-boundary value problem for the nonlinear Klein-Gordon equation, J. Math. Phys. 36 (1995), pp.756-762.
  * [23] Ph. Chartier, N. Crouseilles, M. Lemou, and F. Méhats, Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations, Numer. Math. 129 (2015), pp. 211-250.
  * [24] F. Collino, P. Monk, The perfectly matched layer in curvilinear coordinates, SIAM J. Sci. Comput. 19 (1998), pp. 2061-2090.
  * [25] Z.M. Chen, H.J. Wu, An adaptive finite element method with perfectly matched absorbing layers for the wave scattering by periodic structures, SIAM J. Numer. Anal. 41 (2003), pp. 799–826.
  * [26] A.S. Davydov, Quantum Mechanics, 2nd Edition. Pergamon Press 1976.
  * [27] X Dong, Z Xu, X Zhao, On time-splitting pseudospectral discretization for nonlinear Klein-Gordon equation in nonrelativistic limit regime, Commun. Comput. Phys. 16 (2014), pp. 440-466.
  * [28] J. Ginibre, G. Velo, The global Cauchy problem for the nonlinear Klein-Gordon equation-II, Ann. Inst. Henri Poincaré 6 (1989), pp. 15-35.
  * [29] T. Hagstrom, Perfectly matched layers for hyperbolic systems with applications to the linearized Euler equations, in: Cohen et al. (Ed.), Mathematical and Numerical Aspects of Wave Propagation, Proceedings of the Waves 2003, Springer Verlag, 2003, pp. 125-129.
  * [30] M. Hochbruck, A. Ostermann, Exponential integrators, Acta Numer. 19 (2010), pp. 209-286.
  * [31] F.Q. Hu, X.D. Li, D.K. Lin, Absorbing boundary conditions for nonlinear Euler and Navier–Stokes equations based on the perfectly matched layer technique, J. Comput. Phys. 227 (2008), pp. 4398-4424.
  * [32] S. Johnson, Notes on Perfectly Matched Layers (PMLs), Tech. report, MIT, 2010.
  * [33] S. Machihara, K. Nakanishi, T. Ozawa, Nonrelativistic limit in the energy space for nonlinear Klein-Gordon equations, Math. Ann. 322 (2002), pp. 603-621.
  * [34] N. Masmoudi, K. Nakanishi, From nonlinear Klein-Gordon equation to a system of coupled nonlinear Schrödinger equations, Math. Ann. 324 (2002), pp. 359-389.
  * [35] N. Mauser, Y. Zhang, X. Zhao, On the rotating nonlinear Klein-Gordon equation: non-relativistic limit and numerical methods, Multiscale Model. Simul. 18 (2020), pp. 999-1024.
  * [36] P.J. Pascual, S. Jiménez, L. Vázquez, Numerical simulations of a nonlinear Klein-Gordon model. Applications. Computational Physics (Granada, 1994). Lecture Notes in Physics, vol. 448. pp. 211-270. Springer, Berlin (1995).
  * [37] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed. SIAM, 2003.
  * [38] J. Shen, T. Tang, L. Wang, Spectral Methods: Algorithms, Analysis and Applications, Springer, 2011.
  * [39] K. Schratz, X. Zhao, On the comparison of asymptotic expansion techniques for the nonlinear Klein-Gordon equation in the nonrelativistic limit regime, Discrete Contin. Dyn. Syst. Ser. B 25 (2020), pp. 2841-2865.
  * [40] X.F. Sun, Z.H. Jiang, X.W. Hu, G. Zhuang, J.F. Jiang, W.X. Guo, Perfectly matched layer absorbing boundary condition for nonlinear two-fluid plasma equations, J. Comput. Phys. 286 (2015), pp. 128-142.
  * [41] S. Weinberg, The Quantum Theory of Fields, Volume 1: Foundations. Cambridge University Press, 1995.
  * [42] C. Xiong, M. Good, Y. Guo, X. Liu, K. Huang, Relativistic superfluidity and vorticity from the nonlinear Klein-Gordon equation, Phys. Rev. D 90 (2014), pp. 125019.
  * [43] C. Zheng, A perfectly matched layer approach to the nonlinear Schrödinger wave equations, J. Comput. Phys. 227 (2007), pp. 537-556.