$\begin{split}\widetilde{F}^{j}=\widetilde{U}^{j}&+\frac{\alpha_{1}}{2\kappa_{1}}\boldsymbol{e}^{(1)}_{0}\otimes\left(\widetilde{\boldsymbol{\mathcal{B}}}^{j+1}_{r}+\widetilde{\boldsymbol{\mathcal{B}}}^{j+1}_{l}\right)^{\intercal}\mathcal{D}^{\intercal}_{2}+\frac{\alpha_{1}}{2(\kappa_{1}+1)}\boldsymbol{e}^{(1)}_{1}\otimes\left(\widetilde{\boldsymbol{\mathcal{B}}}^{j+1}_{r}-\widetilde{\boldsymbol{\mathcal{B}}}^{j+1}_{l}\right)^{\intercal}\mathcal{D}^{\intercal}_{2}\\\
&+\frac{\alpha_{2}}{2\kappa_{2}}\mathcal{D}_{1}\left(\widetilde{\boldsymbol{\mathcal{B}}}^{j+1}_{t}+\widetilde{\boldsymbol{\mathcal{B}}}^{j+1}_{b}\right)\otimes\left(\boldsymbol{e}_{0}^{(2)}\right)^{\intercal}+\frac{\alpha_{2}}{2(\kappa_{2}+1)}\mathcal{D}_{1}\left(\widetilde{\boldsymbol{\mathcal{B}}}^{j+1}_{t}-\widetilde{\boldsymbol{\mathcal{B}}}^{j+1}_{b}\right)\otimes\left(\boldsymbol{e}_{1}^{(2)}\right)^{\intercal}+\widetilde{G}^{j+1},\end{split}$
(168)

where $\mathcal{D}_{k}:1-\alpha_{k}^{-2}\partial^{2}_{y_{k}},\;k=1,2,$ are
symbolic representations of the second order differential operators333We would
like to mention that the second order derivative of the history functions will
be handled by first expanding them in the Legendre basis and then using the
backward recursion in the frequency space to obtain the expansion coefficients
for $\partial^{2}_{y_{k}}{\mathcal{B}}^{j+1}(y_{k}),\;k=1,2$. For more details
the reader is referred to Shen [24, Chap. 3]. and

$\widetilde{G}^{j+1}=(\widetilde{g}^{j+1}_{p_{1},p_{2}})=\begin{cases}\frac{\alpha_{1}\alpha_{2}}{4\kappa_{1}(\kappa_{2}+1)}\left(\mathcal{C}^{j+1}_{rt}+\mathcal{C}^{j+1}_{lb}+\mathcal{C}^{j+1}_{rb}+\mathcal{C}^{j+1}_{lt}\right),&p_{1}=0,p_{2}=0,\\\
\frac{\alpha_{1}\alpha_{2}}{4\kappa_{1}(\kappa_{2}+1)}\left(\mathcal{C}^{j+1}_{rt}-\mathcal{C}^{j+1}_{lb}-\mathcal{C}^{j+1}_{rb}+\mathcal{C}^{j+1}_{lt}\right),&p_{1}=0,p_{2}=1,\\\
\frac{\alpha_{1}\alpha_{2}}{4\kappa_{2}(\kappa_{1}+1)}\left(\mathcal{C}^{j+1}_{rt}-\mathcal{C}^{j+1}_{lb}+\mathcal{C}^{j+1}_{rb}-\mathcal{C}^{j+1}_{lt}\right),&p_{1}=1,p_{2}=0,\\\
\frac{\alpha_{1}\alpha_{2}}{4(\kappa_{1}+1)(\kappa_{2}+1)}\left(\mathcal{C}^{j+1}_{rt}+\mathcal{C}^{j+1}_{lb}-\mathcal{C}^{j+1}_{rb}-\mathcal{C}^{j+1}_{lt}\right),&p_{1}=1,p_{2}=1,\\\
0,&p_{1}=2,\ldots,N_{1},\;p_{2}=2,\ldots,N_{2},\end{cases}$ (169)

The variational formulation in (163) now takes the form

$\begin{split}&-\alpha_{1}^{-2}\left(\partial^{2}_{y_{1}}w^{j+1},\theta_{p_{1},p_{2}}\right)_{\Omega_{i}^{\text{ref}}}-\alpha^{-2}_{2}\left(\partial^{2}_{y_{2}}w^{j+1},\theta_{p_{1},p_{2}}\right)_{\Omega_{i}^{\text{ref}}}+\left(w^{j+1},\theta_{p_{1},p_{2}}\right)_{\Omega_{i}^{\text{ref}}}=\left(f^{j},\theta_{p_{1},p_{2}}\right)_{\Omega_{i}^{\text{ref}}}.\end{split}$
(170)

The inner product $(f^{j},\theta_{p_{1},p_{2}})_{\Omega_{i}^{\text{ref}}}$ can
be computed from the knowledge of the Legendre coefficients
$\widetilde{F}^{j}$ using the specific form of the boundary adapted basis. In
B, it is shown that the aforementioned operation can be achieved via
_quadrature matrices_ ($Q_{j}$) given by

$Q_{j}=B_{j}^{\intercal}=\begin{pmatrix}1&0&b^{(j)}_{0}&&&&\\\
&1&0&b^{(j)}_{1}&&&\\\ &&1&0&b^{(j)}_{2}&&\\\ &&&\ddots&\ddots&\ddots&\\\
&&&&1&0&b^{(j)}_{N-2}\\\ \end{pmatrix}\in\mathbb{C}^{(N-1)\times(N+1)},\quad
j=1,2.$ (171)

The linear system thus becomes

$\begin{split}&\alpha_{1}^{-2}S_{1}^{\intercal}\widehat{W}^{j+1}M_{2}+\alpha^{-2}_{2}M_{1}^{\intercal}\widehat{W}^{j+1}S_{2}+M_{1}^{\intercal}\widehat{W}^{j+1}M_{2}=Q_{1}\Gamma\widetilde{F}^{j}\Gamma
Q_{2}^{\intercal}\equiv\widehat{F},\\\
&\left(\alpha_{1}^{-2}M_{2}^{\intercal}\otimes
S_{1}^{\intercal}+\alpha_{2}^{-2}S_{2}^{\intercal}\otimes
M_{1}^{\intercal}+M_{2}^{\intercal}\otimes
M_{1}^{\intercal}\right)\widehat{\boldsymbol{w}}^{j+1}=\widehat{\boldsymbol{f}}^{j},\end{split}$
(172)

where $\otimes$ denotes the tensor product and
$\widehat{\boldsymbol{w}},\;\widehat{\boldsymbol{f}}$ denote the column
vectors obtained by stacking columns of matrices $\widehat{W},\;\widehat{F}$,
respectively, below one another. The mass and stiffness matrices are defined
in (161) and the linear system is then solved using LU-decomposition method.
In order to understand the sparsity pattern of the system matrix, we proceed
as follows: Taking into account the symmetric nature of the mass matrices
(noting stiffness matrices are already diagonal), we have

$\left(\alpha_{1}^{-2}M_{2}^{\intercal}\otimes
S_{1}^{\intercal}+\alpha_{2}^{-2}S_{2}^{\intercal}\otimes
M_{1}^{\intercal}+M_{2}^{\intercal}\otimes
M_{1}^{\intercal}\right)=\left(\alpha_{2}^{-2}S_{2}+M_{2}\right)\otimes\left(\alpha_{1}^{-2}S_{1}+M_{1}\right)-\alpha_{2}^{-2}\alpha_{1}^{-2}\left(S_{2}\otimes
S_{1}\right).$ (173)

From the property of Kronecker product, it follows that the system matrix is
symmetric block pentadiagonal because the individual matrices in the Kronecker
product are symmetric pentadiagonal in nature. The upper and lower bandwidth
works out to be $2N_{1}$. The cost of LU decomposition can be worked out to be
$\mathop{\mathscr{O}}\left(N^{2}_{1}N_{\text{dim}}\right)$ where
$N_{\text{dim}}=(N_{1}-1)(N_{2}-1)$. The lower and upper triangular matrices
inherit the bandedness of the system matrix provided pivoting is not required
[23]. Under this assumption, the resulting complexity of solving the linear
system (excluding the cost of LU decomposition) becomes
$\mathop{\mathscr{O}}\left(N_{1}N_{\text{dim}}\right)$. A thorough study of
this issue is beyond the scope of this paper. For a range of values of $\Delta
t$, we note that MATLAB arrives at banded lower and upper triangular matrices
which is consistent with lack of pivoting.

###  3.3 Algorithmic summary

In this section, we present a brief algorithmic summary for the CQ and NP
methods in context of BDF1 based temporal discretization of the boundary maps
and also the interior problem. For the TR case, the algorithmic steps are
largely similar.

####  3.3.1 CQ–BDF1

Let us recall that, on the boundary segments $\Gamma_{a_{1}}$ and
$\Gamma_{a_{2}}$, the auxiliary function is denoted by
$\varphi_{a_{1}}(y_{2},\tau_{1},\tau_{2})$ and
$\varphi_{a_{2}}(y_{1},\tau_{1},\tau_{2})$, respectively, while the auxiliary
function at the corners is denoted by
$\varphi_{a_{1}a_{2}}(\tau_{1},\tau_{2})$ where
$a_{1}\in\\{l,r\\},\;a_{2}\in\\{b,t\\}$. Let $N$ be the number of spatial
nodes along $y_{1}$ as well as $y_{2}$, i.e., $N_{1}=N_{2}=N$. Let the number
of time-steps be denoted by $N_{t}$. The temporal variables $t$, $\tau_{1}$
and $\tau_{2}$ are discretized with the same step-size $\Delta t$ so that the
maximum time $T=N_{t}\Delta t$. The storage of the auxiliary field on the
boundary segments require four 3D-arrays of size $N\times N_{t}\times N_{t}$
and the storage of auxiliary field at corner points require four 2D-arrays of
size $N_{t}\times N_{t}$ which corresponds to
$\varphi_{a_{1}a_{2}}(\tau_{1},\tau_{2})$. The algorithmic steps for the
solving the IBVP are enumerated below.

  1. Step 1:

Advance the following IVPs corresponding to auxiliary functions by one-step,
say, $j\rightarrow(j+1)$ using the previously computed values of the auxiliary
field on the segments (as described in Fig. 2):

$\begin{split}&-\alpha^{-2}_{1}\partial_{y_{1}}^{2}\varphi^{j+1,q}_{a_{2}}+\varphi^{j+1,q}_{a_{2}}=\varphi^{j,q}_{a_{2}},\quad
q=0,1,\ldots,j,\\\
&-\alpha^{-2}_{2}\partial_{y_{2}}^{2}\varphi^{m,j+1}_{a_{1}}+\varphi^{m,j+1}_{a_{1}}=\varphi^{m,j}_{a_{1}},\quad
m=0,1,\ldots,j,\end{split}$

using the TBCs at the corner points given by

$(\partial_{n_{1}}+\alpha_{1})\varphi^{j+1,q}_{a_{2}}(y_{a_{1}})=-\alpha_{1}\mathcal{B}^{j+1,q}_{a_{2},a_{1}},\quad(\partial_{n_{2}}+\alpha_{2})\varphi^{m,j+1}_{a_{1}}(y_{a_{2}})=-\alpha_{2}\mathcal{B}^{m,j+1}_{a_{1},a_{2}},$

where the history functions are computed as

$\mathcal{B}^{j+1,q}_{a_{2},a_{1}}=\sum_{k=1}^{j+1}\omega_{k}\varphi^{j+1-k,q}_{a_{1}a_{2}},\quad\mathcal{B}^{m,j+1}_{a_{1},a_{2}}=\sum_{k=1}^{j+1}\omega_{k}\varphi^{m,j+1-k}_{a_{1}a_{2}}.$

Let the cost of solving the linear 1D system of the kind (142) be denoted by
$\Theta_{1}(N_{\text{dim}})$ where the dimension is $N_{\text{dim}}=N-1$. The
cost of the steps outlined here works out to be
$4(j+1)\Theta_{1}(N-1)+4(j+1)^{2}$ where the second part is the cost of
computing the history functions. The memory requirement is
$4(j+1)^{2}N+4(j+1)^{2}$ size array. Therefore, the complexity scales with
$j$. Moreover, for large number of time-steps, i.e., $N_{t}\gg N$, the term
$(j+1)^{2}$ starts dominating after a point.

  2. Step 2:

Compute the history functions on the boundary segments of the domain
$\Omega_{i}$ as

$\mathcal{B}^{j+1}_{a_{1}}(y_{2})=\sum_{k=1}^{j+1}\omega_{k}\varphi^{j+1-k,j+1}_{a_{1}}(y_{2}),\quad\mathcal{B}^{j+1}_{a_{2}}(y_{1})=\sum_{k=1}^{j+1}\omega_{k}\varphi^{j+1,j+1-k}_{a_{2}}(y_{1}).$

The cost of this step is $\mathop{\mathscr{O}}\left((j+1)N\right)$.

  3. Step 3:

Solve the linear system for the interior problem to compute the solution at
$(j+1)$-th time step

$-\left(\alpha^{-2}_{1}\partial^{2}_{y_{1}}+\alpha^{-2}_{2}\partial^{2}_{y_{2}}\right)u^{j+1}+u^{j+1}=u^{j},\quad\rho=1/\Delta
t.$

using the following Robin-type TBCs at discrete level

$\partial_{n_{1}}u^{j+1}+\alpha_{1}u^{j+1}=-\alpha_{1}\mathcal{B}^{j+1}_{a_{1}}(y_{2}),\quad\partial_{n_{2}}u^{j+1}+\alpha_{2}u^{j+1}=-\alpha_{2}\mathcal{B}^{j+1}_{a_{2}}(y_{2}).$

Let the cost of solving the linear 2D system of the kind (172) be denoted by
$\Theta_{2}(N_{\text{dim}})$ where the dimension is
$N_{\text{dim}}=(N-1)^{2}$. The cost of solving the linear system in this step
works out to be
$\Theta_{2}(N_{\text{dim}})=\mathop{\mathscr{O}}\left(N^{3}\right)$ provided
pivoting is not needed.

  4. Step 4:

Update the 2D and 3D arrays storing the values of the auxiliary field with the
$(j+1)$-th time step values using the fact that
$\varphi^{j+1,j+1}(y_{1},y_{2})=u^{j+1}(y_{1},y_{2})$.

For moderate value of $N_{t}$, i.e.,
$N_{t}=\mathop{\mathscr{O}}\left(N\right)$, the cost of computation at each
time-step is dominated by
$\Theta_{2}(N_{\text{dim}})=\mathop{\mathscr{O}}\left(N^{3}\right)$ which
corresponds to the linear system for the interior problem. The memory
requirement for the auxiliary functions in this case scales as
$\mathop{\mathscr{O}}\left(N^{3}\right)$ which is quite prohibitive.

In the regime $N_{t}\gg N$, say,
$N_{t}=\mathop{\mathscr{O}}\left(N^{2}\right)$, the cost of computation scales
as $\mathop{\mathscr{O}}\left(N^{4}\right)$ which corresponds to cost of
computing the history functions while the memory requirement scales as
$\mathop{\mathscr{O}}\left(N^{5}\right)$ .

####  3.3.2 NP–BDF1

Let us recall that, on the boundary segments $\Gamma_{a_{1}}$ and
$\Gamma_{a_{2}}$, the auxiliary fields are denoted by
$\varphi_{k,a_{1}}(y_{2},\tau_{1},\tau_{2})$ and
$\varphi_{k,a_{2}}(y_{1},\tau_{1},\tau_{2})$, respectively, while the
auxiliary field at the corners is denoted by
$\psi_{k,k^{\prime},a_{1},a_{2}}(\tau_{1},\tau_{2})$ where
$a_{1}\in\\{l,r\\},\;a_{2}\in\\{b,t\\}$ and $k,k^{\prime}=1,2,\ldots,M$. Let
$N$ be the number of spatial nodes along $y_{1}$ as well as $y_{2}$, i.e.,
$N_{1}=N_{2}=N$. The storage of the auxiliary fields on the boundary segments
require four 2D-arrays of size $M\times N$ and the storage of auxiliary field
at corner points require four 2D-arrays of size $M\times M$ where $M(<N)$ is
the order of the diagonal Padé approximants used in the NP method. The
evolution of these fields are effectively local in time which circumvents the
need to store the history of these fields in time. The algorithmic steps for
the solving the IBVP are enumerated below.

  1. Step 1:

Advance the following IVPs corresponding to auxiliary fields by one step ,
say, $j\rightarrow(j+1)$ using the previously computed values of the auxiliary
field on the segments (as described in Fig. 3):

$-\alpha^{-2}_{1}\partial_{y_{1}}^{2}\varphi^{j+1,j}_{k,a_{2}}+\varphi^{j+1,j}_{k,a_{2}}=\varphi^{j,j}_{k,a_{2}},\quad-\alpha^{-2}_{1}\partial_{y_{2}}^{2}\varphi^{j,j+1}_{k,a_{1}}+\varphi^{j,j+1}_{k,a_{1}}=\varphi^{j,j}_{k,a_{1}},\quad
k=1,2,\ldots,M,$

together with the boundary conditions

$\partial_{n_{2}}\varphi_{k,a_{1}}^{j,j+1}+\alpha_{2}\varpi\varphi_{k,a_{1}}^{j,j+1}=-\alpha_{2}\mathcal{B}^{j+1}_{k,a_{1},a_{2}},\quad\partial_{n_{1}}\varphi_{k,a_{2}}^{j+1,j}+\alpha_{1}\varpi\varphi_{k,a_{2}}^{j+1,j}=-\alpha_{1}\mathcal{B}^{j+1}_{k,a_{2},a_{1}},$

where the history functions are given by

$\mathcal{B}^{j+1}_{k,a_{1},a_{2}}=\sum_{k^{\prime}=1}^{M}\Gamma_{k^{\prime}}\psi_{k,k^{\prime},a_{1},a_{2}}^{j,j},\quad\mathcal{B}^{j+1}_{k,a_{2},a_{1}}=\sum_{k^{\prime}=1}^{M}\Gamma_{k^{\prime}}\psi_{k,k^{\prime},a_{2},a_{1}}^{j,j}=\sum_{k^{\prime}=1}^{M}\Gamma_{k^{\prime}}\psi_{k^{\prime},k,a_{1},a_{2}}^{j,j}.$

These boundary condition are taken from Sec. 3.1 where precise definition of
the quantities are provided. The cost of the steps outlined here works out to
be $4M\Theta_{1}(N-1)+4M^{2}$ where the second part is the cost of computing
the history functions. The memory requirement is $4MN+4M^{2}$ size array.

  2. Step 2:

Compute the history functions on the boundary segments of the domain
$\Omega_{i}$ as

$\mathcal{B}^{j+1}_{a_{1}}(y_{2})=\sum_{k=1}^{M}\Gamma_{k}\varphi_{k,a_{1}}^{j,j+1},\quad\mathcal{B}^{j+1}_{a_{2}}(y_{1})=\sum_{k=1}^{M}\Gamma_{k}\varphi_{k,a_{2}}^{j+1,j}.$

The cost of this step is $\mathop{\mathscr{O}}\left(MN\right)$.

  3. Step 3:

Solve the linear system for the interior problem to compute the solution at
$(j+1)$-th time step

$-\left(\alpha^{-2}_{1}\partial^{2}_{y_{1}}+\alpha^{-2}_{2}\partial^{2}_{y_{2}}\right)u^{j+1}+u^{j+1}=u^{j},\quad\rho=1/\Delta
t.$

using the following Robin-type TBCs at discrete level

$\partial_{n_{1}}u^{j+1}+\alpha_{1}\varpi
u^{j+1}=-\alpha_{1}\mathcal{B}^{j+1}_{a_{1}}(y_{2}),\quad\partial_{n_{2}}u^{j+1}+\alpha_{2}\varpi
u^{j+1}=-\alpha_{2}\mathcal{B}^{j+1}_{a_{2}}(y_{1}).$

The cost of solving the linear system in this step works out to be
$\Theta_{2}(N_{\text{dim}})=\mathop{\mathscr{O}}\left(N^{3}\right)$ provided
pivoting is not needed.

  4. Step 4:

Update the 2D arrays storing the values of the auxiliary fields with the
$(j+1)$-th time step values using the following update relations

$\left\\{\begin{aligned}
&\varphi^{j+1,j+1}_{k,a_{1}}=\frac{1}{(1+\overline{\eta}^{2}_{k})}\varphi^{j,j+1}_{k,a_{1}}+\frac{1/\rho}{(1+\overline{\eta}^{2}_{k})}u^{j+1},\\\
&\varphi^{j+1,j+1}_{k,a_{2}}=\frac{1}{(1+\overline{\eta}^{2}_{k})}\varphi^{j+1,j}_{k,a_{2}}+\frac{1/\rho}{(1+\overline{\eta}^{2}_{k})}u^{j+1},\end{aligned}\right.\quad\text{and}\quad\left\\{\begin{aligned}
&\psi^{j,j+1}_{k,k^{\prime},a_{1}a_{2}}=\frac{1}{(1+\overline{\eta}^{2}_{k^{\prime}})}\psi^{j,j}_{k,k^{\prime},a_{1},a_{2}}+\frac{1/\rho}{(1+\overline{\eta}^{2}_{k^{\prime}})}\varphi_{k,a_{1}}^{j,j+1},\\\
&\psi^{j+1,j+1}_{k,k^{\prime},a_{1}a_{2}}=\frac{1}{(1+\overline{\eta}^{2}_{k})}\psi^{j,j+1}_{k,k^{\prime},a_{1},a_{2}}+\frac{1/\rho}{(1+\overline{\eta}^{2}_{k})}\varphi_{k^{\prime},a_{2}}^{j+1,j+1}.\end{aligned}\right.$

From the discussion above, it is evident that the overall cost of computation
is completely independent of the time-step and the DtN-maps get implemented at
the cost of $\mathop{\mathscr{O}}\left(MN\right)$. The memory requirement for
the auxiliary functions also remain independent of time-step with an estimate
of $\mathop{\mathscr{O}}\left(MN\right)$. Therefore, the overall complexity is
dominated by
$\Theta_{2}(N_{\text{dim}})=\mathop{\mathscr{O}}\left(N^{3}\right)$ which
corresponds to the linear system for the interior problem.

Figure 4: The figure shows the evolution of the relative energy content as
defined in (177) of the chirped-Gaussian and the Hermite-Gaussian profiles
considered in Table 1. Here the computational domain is
$\Omega_{i}=(-10,10)^{2}$.

##  4 Numerical Experiments

In this section, we will carry out the extensive numerical tests to showcase
the accuracy of the numerical schemes developed in this work to solve the IBVP
in (14). We start by analyzing the behaviour of the exact solutions for the
IBVP under consideration followed by studying the error evolution behaviour
and convergence analysis of the numerical schemes.

###  4.1 Exact solutions

The exact solutions admissible for our numerical experiments are such that the
initial profile must be effectively supported within the computational domain.
We primarily consider the wavepackets which are a modulation of a Gaussian
envelop so that the requirement of effective initial support can be easily
met. In this class of solutions, we consider the Chirped-Gaussian and the
Hermite-Gaussian profiles.

Table 1:  Chirped-Gaussian profile with $A_{0}=2$ and
$c_{0}\in\\{4,8,12,16\\}$.
$G\left(\boldsymbol{x},t\right)=A_{0}\sum_{j=1}^{n}G(\boldsymbol{x},t;\boldsymbol{a}_{j},\boldsymbol{b}_{j},\boldsymbol{c}_{j}),\quad\boldsymbol{c}_{j}=c_{0}(\cos\theta_{j},\sin\theta_{j})$
---
Type |  $n$ |  $(\boldsymbol{a}_{j}\in\mathbb{R}^{2})_{j=1}^{n}$ |  $(\boldsymbol{b}_{j}\in\mathbb{R}^{2})_{j=1}^{n}$ |  $\boldsymbol{\theta}\in\mathbb{R}^{n}$
IA,
IB |  $2$ |  $\boldsymbol{a}_{1}=(1/2.5,1/2.4)$,
$\boldsymbol{a}_{2}=(1/2.3,1/2.2)$ |  $\boldsymbol{b}_{j}=(1/2)\boldsymbol{1}$ |  IA: $\boldsymbol{\theta}_{A}=(0,\pi)$,
IB: $\boldsymbol{\theta}_{B}=\boldsymbol{\theta}_{A}+(\pi/4)\boldsymbol{1}$
IIA,
IIB |  $4$ |  $\boldsymbol{a}_{1}=(1/2.5,1/2.4)$,
$\boldsymbol{a}_{2}=(1/2.3,1/2.2)$,
$\boldsymbol{a}_{3}=(1/2.7,1/2.6)$,
$\boldsymbol{a}_{4}=(1/2.2,1/2.5)$ |  $\boldsymbol{b}_{j}=(1/2)\boldsymbol{1}$ |  IIA: $\boldsymbol{\theta}_{A}=(0,\pi/2,\pi,3\pi/2)$,
IIB: $\boldsymbol{\theta}_{B}=\boldsymbol{\theta}_{A}+(\pi/4)\boldsymbol{1}$

####  4.1.1 Chirped-Gaussian profile

Using the function defined by

$\mathcal{G}(x,t;a,b)=\frac{1}{\sqrt{1+4i(a+ib)t}}\exp\left[-\frac{(a+ib)}{1+4i(a+ib)t}x^{2}\right],\quad
a>0,\;b\in\mathbb{R},$ (174)

one can define the family of solutions referred to as chirped-Gaussian profile
by

$G(\boldsymbol{x},t;\boldsymbol{a},\boldsymbol{b},\boldsymbol{c})=\mathcal{G}(x_{1}-c_{1}t,t;a_{1},b_{1})\mathcal{G}(x_{2}-c_{2}t,t;a_{2},b_{2})\exp\left(+i\frac{1}{2}\boldsymbol{c}\cdot\boldsymbol{x}-i\frac{1}{4}\boldsymbol{c}\cdot\boldsymbol{c}\,t\right),$
(175)

where $\boldsymbol{a}\in\mathbb{R}^{2}_{+}$ determines the effective support
of the profile at $t=0$, $\boldsymbol{b}\in\mathbb{R}^{2}$ is the _chirp_
parameter and $\boldsymbol{c}\in\mathbb{R}^{2}$ is the velocity of the
profile. Using linear combination, one can further define a more general
family of solutions with parameters $A_{0},c_{0}\in\mathbb{R}$,
$\boldsymbol{\theta}\in\mathbb{R}^{n}$,
$(\boldsymbol{a}_{j}\in\mathbb{R}^{2}_{+})_{j=1}^{n}$ and
$(\boldsymbol{b}_{j}\in\mathbb{R}^{2})_{j=1}^{n}$ given by

$G\left(\boldsymbol{x},t;\;c_{0},A_{0};\;(\boldsymbol{a}_{j})_{j=1}^{n},(\boldsymbol{b}_{j})_{j=1}^{n},\boldsymbol{\theta}\right)=A_{0}\sum_{j=1}^{n}G(\boldsymbol{x},t;\boldsymbol{a}_{j},\boldsymbol{b}_{j},\boldsymbol{c}_{j}),\quad\boldsymbol{c}_{j}=c_{0}(\cos\theta_{j},\sin\theta_{j}).$
(176)

Let the constant vector $(1,1,\ldots)\in\mathbb{R}^{n}$ be denoted by
$\boldsymbol{1}$, then the specific values of the parameters of the solutions
used in the numerical experiments can be summarized as in Table 1. The energy
content of the profile within the computational domain $\Omega_{i}$ over time
is

$E_{\Omega_{i}}(t)=\left.{\int_{\Omega_{i}}|G(\boldsymbol{x},t)|^{2}d^{2}\boldsymbol{x}}\right/{\int_{\Omega_{i}}|G(\boldsymbol{x},0)|^{2}d^{2}\boldsymbol{x}},\quad
t\geq 0.$ (177)

The profiles are chosen with non-zero speed $c_{0}$ so that the field hits the
boundary of $\Omega_{i}$, In case of a square computational domain, the type
‘A’ class of solutions are directed to the segments of the boundary normally
while type ‘B’ class of solutions are directed to the corners. The behaviour
of the $E_{\Omega_{i}}(t)$ for $t\in[0,5]$ is shown in Fig. 4.

Table 2:  Hermite-Gaussian profile with $A_{0}=2$ and
$c_{0}\in\\{4,8,12,16\\}$.
$G\left(\boldsymbol{x},t\right)=A_{0}\sum_{j=1}^{n}G(\boldsymbol{x},t;\boldsymbol{m}_{j},\boldsymbol{a}_{j},\boldsymbol{c}_{j}),\;\boldsymbol{c}_{j}=c_{0}(\cos\theta_{j},\sin\theta_{j})$
---
Type |  $n$ |  $(\boldsymbol{a}_{j}\in\mathbb{R}^{2})_{j=1}^{n}$ |  $(\boldsymbol{m}_{j}\in\mathbb{N}_{0}^{2})_{j=1}^{n}$ |  $\boldsymbol{\theta}\in\mathbb{R}^{n}$
IA,
IB |  $2$ |  $\boldsymbol{a}_{1}=(1/2.5,1/2.4)$,
$\boldsymbol{a}_{2}=(1/2.3,1/2.2)$ |  $\boldsymbol{m}_{1}=(1,2)$,
$\boldsymbol{m}_{2}=(2,1)$ |  $\boldsymbol{\theta}_{A}=(0,\pi)$,
$\boldsymbol{\theta}_{B}=\boldsymbol{\theta}_{A}+(\pi/4)\boldsymbol{1}$
IIA,
IIB |  $4$ |  $\boldsymbol{a}_{1}=(1/2.5,1/2.4)$,
$\boldsymbol{a}_{2}=(1/2.3,1/2.2)$,
$\boldsymbol{a}_{3}=(1/2.7,1/2.6)$,
$\boldsymbol{a}_{4}=(1/2.2,1/2.5)$ |  $\boldsymbol{m}_{1}=(1,2)$,
$\boldsymbol{m}_{2}=(2,1)$,
$\boldsymbol{m}_{3}=(2,1)$,
$\boldsymbol{m}_{4}=(1,2)$ |  $\boldsymbol{\theta}_{A}=(0,\pi/2,\pi,3\pi/2)$,
$\boldsymbol{\theta}_{B}=\boldsymbol{\theta}_{A}+(\pi/4)\boldsymbol{1}$
Figure 5: The figure shows the behaviour of error involved in the
discretization of the boundary map alone for the chirped-Gaussian profile with
different values of the speed ‘c’ (see Table 1). The numerical parameters and
the labels are described in Sec. 4.2 where the error is quantified by (184).
Figure 6: The figure shows the behaviour of error involved in the
discretization of the boundary map alone for the Hermite-Gaussian profile with
different values of the speed ‘c’ (see Table 2). The numerical parameters and
the labels are described in Sec. 4.2 where the error is quantified by (184).
Figure 7: The figure shows a comparison of evolution of error in the numerical
solution of the IBVP (14) with the various TBCs for the chirped-Gaussian
profile with different values of the speed ‘c’ (see Table 1). The numerical
parameters and the labels are described in Sec. 4.3 where the error is
quantified by (185).

####  4.1.2 Hermite-Gaussian profile

Consider the class of normalized Hermite-Gaussian functions defined by

$\mathcal{G}_{m}(x,t;a)=\gamma_{m}^{-1}H_{m}\left(\frac{\sqrt{2a}\,x}{w(t)}\right)\sqrt{\frac{\mu(t)}{a}}\exp\left[-\mu(t)x^{2}-i\,m\theta(t)\right],\quad
m\in\mathbb{N}_{0}=\\{0,1,\ldots\\},$ (178)

where $a>0$,

$w(t)=\sqrt{1+(4at)^{2}},\quad\frac{1}{\mu(t)}=\frac{1}{a}+i{4t}=\frac{1}{a}w(t)\exp[i\theta(t)],$
(179)

and the normalization factor is given by
$\gamma^{2}_{m}=2^{m}(m!)\sqrt{\pi}(2a)^{-1/2}$.

The Hermite polynomials are evaluated using the following relations

$H_{n+1}(x)=2xH_{n}(x)-2nH_{n-1}(x),\quad H_{n+1}^{\prime}(x)=2(n+1)H_{n}(x),$
(180)

with $H_{0}(x)=1$ and $H_{1}(x)=2x$, These relations can be used to compute
the first spatial derivative of the Hermite-Gaussian functions:

$\partial_{x}\mathcal{G}_{m}(x,t;a)=-\sqrt{(m+1)a}\mathcal{G}_{m+1}(x,t;a)+\sqrt{ma}\mathcal{G}_{m-1}(x,t;a),\quad
m>0,$ (181)

with $\partial_{x}\mathcal{G}_{0}(x,t;a)=-\sqrt{a}\mathcal{G}_{1}(x,t;a)$.
This is required for testing our Dirichlet-to-Neumann maps. Using these
functions, we can define a family of solutions referred to as Hermite-Gaussian
profile by

$G(\boldsymbol{x},t;\boldsymbol{m},\boldsymbol{a},\boldsymbol{c})=\mathcal{G}_{m_{1}}(x_{1}-c_{1}t,t;a_{1})\mathcal{G}_{m_{2}}(x_{2}-c_{2}t,t;a_{2})\exp\left(+i\frac{1}{2}\boldsymbol{c}\cdot\boldsymbol{x}-i\frac{1}{4}\boldsymbol{c}\cdot\boldsymbol{c}\,t\right),$
(182)

where $\boldsymbol{m}\in\mathbb{N}_{0}^{2}$ is the order parameter,
$\boldsymbol{a}\in\mathbb{R}^{2}_{+}$ determines the effective support of the
profile at $t=0$ and $\boldsymbol{c}\in\mathbb{R}^{2}$ is the velocity of the
profile. Once again, using linear combination, one can further define a more
general family of solutions with parameters $A_{0},c_{0}\in\mathbb{R}$,
$\boldsymbol{\theta}\in\mathbb{R}^{n}$,
$(\boldsymbol{m}_{j}\in\mathbb{N}^{2}_{0})_{j=1}^{n}$ and
$(\boldsymbol{a}_{j}\in\mathbb{R}^{2}_{+})_{j=1}^{n}$ given by

$G\left(\boldsymbol{x},t;\;c_{0},A_{0};\;(\boldsymbol{m}_{j})_{j=1}^{n},(\boldsymbol{a}_{j})_{j=1}^{n},\boldsymbol{\theta}\right)=A_{0}\sum_{j=1}^{n}G(\boldsymbol{x},t;\boldsymbol{m}_{j},\boldsymbol{a}_{j},\boldsymbol{c}_{j}),\quad\boldsymbol{c}_{j}=c_{0}(\cos\theta_{j},\sin\theta_{j}).$
(183)

As in the last case, we set $\boldsymbol{1}=(1,1,\ldots)\in\mathbb{R}^{n}$,
then the specific values of the parameters of the solutions used in the
numerical experiments can be summarized as in Table 2. The energy content of
the profile within the computational domain $\Omega_{i}$ over time is
$E_{\Omega_{i}}(t)$ as defined by (177) with the appropriate profile in the
integrand. The profiles are chosen with non-zero speed $c_{0}$ so that the
field hits the boundary of $\Omega_{i}$. In case of a square computational
domain, the type ‘A’ class of solutions are directed to the segments of the
boundary normally while type ‘B’ class of solutions are directed to the
corners. The behaviour of the $E_{\Omega_{i}}(t)$ for $t\in[0,5]$ is shown in
Fig. 4.

Figure 8: The figure shows a comparison of evolution of error in the numerical
solution of the IBVP (14) with the various TBCs for the Hermite-Gaussian
profile with different values of the speed ‘c’ (see Table 2). The numerical
parameters and the labels are described in Sec. 4.3 where the error is
quantified by (185).

###  4.2 Tests for DtN maps

Prior to the full implementation of the numerical scheme, one can test the
accuracy of the DtN boundary maps where we quantify the error on
$\partial\Omega_{i}$ by

$e(t_{j})=\left(\int_{\partial\Omega_{i}}\left|\partial_{n}u(\boldsymbol{x},t_{j})-[\partial_{n}u(\boldsymbol{x},t_{j})]_{\text{num.}}\right|^{2}d\sigma(\boldsymbol{x})\right)^{1/2},\quad
t_{j}\in[0,T_{max}],$ (184)

where $d\sigma(\boldsymbol{x})$ denotes the measure along the boundary of the
computational domain $\Omega_{i}$. Along the spatial dimension, the integral
above will be approximated by Gauss quadrature over LGL-points 444The
quadrature nodes and the corresponding weights are computed using quad
precision in order to avoid any numerical artifacts due to round-off errors..
The quantity $[\partial_{n}u(\boldsymbol{x},t_{j})]_{\text{num.}}$ refers to
the numerical approximation of the Neumann datum from the Dirichlet datum. All
the intermediate quantities involved in this process is computed according to
a numerical scheme developed as part of the numerical realization of the DtN
map.

In this work, we have developed two discrete versions of the DtN maps, namely,
CQ and NP. Each of these methods have a variant determined by the choice of
one-step method used in the temporal discretization so that the complete list
of methods to be tested can be labelled as CQ–BDF1, NP–BDF1 (corresponding to
the one-step method BDF1), and, CQ–TR, NP–TR (corresponding to the one-step
method TR). For the Padé approximant based method ‘NP’, we distinguish the
diagonal approximant of order $20$ from that of $50$ via the labels ‘NP20’ and
‘NP50’. The numerical parameters used in this section are summarized in Table
3. The exact solutions used are described in Sec. 4.1. Let us note that the
accuracy of the DtN maps do not reflect the accuracy of the numerical solution
computed with the TBCs on account of the fact that the interior problem
together with the discretized DtN maps may turn out to be unstable.

Table 3:  Numerical parameters for testing DtN maps Computational domain ($\Omega_{i}$) |  $(-10,10)\times(-10,10)$
---|---
Maximum time ($T_{max}$) |  $2$
No. of time-steps ($N_{t}$) |  $1000+1$
Time-step ($\Delta t$) |  $2\times 10^{-3}=T_{max}/(N_{t}-1)$
Number of LGL-points ($N+1$) |  $200$

The numerical results for error on $\partial\Omega_{i}$ corresponding to the
chirped-Gaussian and the Hermite-Gaussian profiles are shown in Fig. 5 and
Fig. 6, respectively. It can be seen that the diagonal Padé approximant based
methods, namely, NP20 and NP50, perform equally well as compared to that of
the convolution quadrature based method. Moreover, the TR methods perform
better than the BDF1 methods which is clear from the error peaks in figures 5
and 6. The accuracy deteriorates with increase in the value of $c$. For faster
moving profiles corresponding to IA and IIA type solutions, NP-20 becomes
slightly less accurate compared to CQ and NP-50. There is no significant
difference in the results for error on $\partial\Omega_{i}$ corresponding to
the chirped-Gaussian profile and the Hermite-Gaussian profile.

Figure 9: The figure shows a comparison of evolution of error in the numerical
solution of the IBVP (14) with the TBCs developed using novel Padé approach
for the chirped-Gaussian profile with different values of the speed ‘c’ (see
Table 1). The numerical parameters and the labels are described in Sec. 4.3
where the error is quantified by (185). Figure 10: The figure shows a
comparison of evolution of error in the numerical solution of the IBVP (14)
with the TBCs developed using novel Padé approach for the Hermite-Gaussian
profile with different values of the speed ‘c’ (see Table 2). The numerical
parameters and the labels are described in Sec. 4.3 where the error is
quantified by (185).

###  4.3 Tests for evolution error

In this section, we consider the IBVP in (14) where the initial condition
corresponds to the exact solutions described in Sec. 4.1. The numerical
solution is labelled according to that of the boundary maps described in Sec.
4.2. The error in the evolution of the profile computed numerically is
quantified by the relative $\mathsf{L}^{2}(\Omega_{i})$-norm:

$e(t_{j})=\left.\left(\int_{\Omega_{i}}\left|u(\boldsymbol{x},t_{j})-[u(\boldsymbol{x},t_{j})]_{\text{num.}}\right|^{2}d^{2}\boldsymbol{x}\right)^{1/2}\,\right/\left(\int_{\Omega_{i}}\left|u_{0}(\boldsymbol{x})\right|^{2}d^{2}\boldsymbol{x}\right)^{1/2},$
(185)

for $t_{j}\in[0,T_{max}]$. The integral above will be approximated by Gauss
quadrature over LGL-points.

For Padé based discretization (NP) of the DtN maps, it is possible to carry
out evolution over larger number of time-steps. The Table 4 lists the
numerical parameters for the NP20 and NP50 methods. The CQ based methods
cannot be tested with the same parameters on account of the increasing memory
requirements at each time-step. Therefore, we introduce a separate test-case
for the sake of comparison with CQ methods by restricting $N_{t}=1000+1$ and
$T_{max}=2$ while other parameters are same as that of Table 4.

The numerical results for the evolution error on $\Omega_{i}$ corresponding to
the chirped-Gaussian and the Hermite-Gaussian profiles are shown in Fig. 7 and
Fig. 8, respectively. It can be seen that the diagonal Padé approximants based
method, namely, NP20 and NP50, perform equally well as compared to that of the
convolution quadrature based method. Moreover, the TR methods perform better
than the BDF1 methods which is clear from the error peaks in figures 7 and 8.
The accuracy deteriorates with increase in the value of $c$. For faster moving
profiles corresponding to IA and IIA type solutions, NP-20 becomes slightly
less accurate compared to CQ and NP-50. There is no significant difference in
the results for evolution error on $\Omega_{i}$ corresponding to the chirped-
Gaussian profile and the Hermite-Gaussian profile.

The numerical results for the evolution error on $\Omega_{i}$ particularly for
NP methods corresponding to the chirped-Gaussian and Hermite-Gaussian profiles
are shown in Fig. 9 and Fig. 10, respectively. Note that the results for NP
are superior than that of CQ on account of the reduced size ($\Delta t$) of
the time-step which was possible because increasing the number of time-steps
did not significantly increase the computation time as well as memory
requirements unlike that of CQ555Attempts to lower the step-size for the same
maximum duration of time in case of CQ failed due to out-of-memory errors. We
note that these problems may severely limit the practical utility of the CQ
method..

Table 4: Numerical parameters for studying the evolution error for the NP methods Computational domain ($\Omega_{i}$) |  $(-10,10)\times(-10,10)$
---|---
Maximum time ($T_{max}$) |  $5$
No. of time-steps ($N_{t}$) |  $5000+1$
Time-step ($\Delta t$) |  $10^{-3}=T_{max}/(N_{t}-1)$
Number of LGL-points ($(N+1)\times(N+1)$) |  $200\times 200$

In order to visualize the reflections at the boundary, we resort to contour
plots depicting the quantity $f_{\text{mag.}}\log_{10}|u(\boldsymbol{x},t)|$
where the scalar quantity $f_{\text{mag.}}$ is adjusted to resolve the contour
lines properly. The specific snapshots correspond to the equispaced energy
levels in the logarithmic scale that is remaining in the computational domain.
We restrict ourselves to presenting the contour plots only for the case of
type IIB class of solutions which are directed towards the corners of the
domain on account of being most complex among four different class of
solutions. The contour plots depicting corresponding to the chirped-Gaussian
($f_{\text{mag.}}=4$) and the Hermite-Gaussian ($f_{\text{mag.}}=8$) profiles
are shown in Fig. 13 and Fig. 14, respectively. The top row corresponds to the
exact solution for reference. From these plots, we conclude that the results
are similar for each of the types of solutions. There are minimal visible
reflections in the contour plots proving the accuracy of the developed
numerical schemes to implement the TBCs. Moreover, the TR based methods better
resolve the contour levels in comparison to the BDF1 based methods.

Table 5: Numerical parameters for studying the convergence error for the NP methods Computational domain ($\Omega_{i}$) |  $(-10,10)\times(-10,10)$
---|---
Maximum time ($T_{max}$) |  $5$
Set of no. of time-steps ($\mathbb{N}_{t}$) |  $\\{2^{8},2^{9},\ldots,2^{18}\\}$
Time-step |  $\\{\Delta t=T_{max}/(N_{t}-1),\;N_{t}\in\mathbb{N}_{t}\\}$
Number of LGL-points ($(N+1)\times(N+1)$) |  $200\times 200$

###  4.4 Tests for convergence

In this section, we analyze the convergence behaviour of the numerical schemes
for the IBVP in (14) where the initial condition corresponds to the exact
solutions described in Sec. 4.1. The numerical solution is labelled according
to that of the boundary maps described in Sec. 4.2.

The error used to study the convergence behaviour is quantified by the maximum
relative $\mathsf{L}^{2}(\Omega_{i})$-norm:

$e=\max\left\\{e(t_{j})|\;j=0,1,\ldots,N_{t}-1\right\\},$ (186)

for $t_{j}\in[0,T_{max}]$. For Padé based discretization (NP) of the DtN maps,
it is possible to analyze the convergence behaviour over large number of time-
steps. The Table 5 lists the numerical parameters for the NP20 and NP50
methods. The CQ based methods cannot be tested with the same parameters on
account of the increasing memory requirements at each time-step.

Figure 11: The figure depicts the convergence behaviour of the numerical
solution of the IBVP (14) with the TBCs developed using novel Padé approach
for the chirped-Gaussian profile with different values of the speed ‘c’ (see
Table 1). The numerical parameters and the labels are described in Sec. 4.4
where the error is quantified by (186). Figure 12: The figure depicts the
convergence behaviour of the numerical solution of the IBVP (14) with the TBCs
developed using novel Padé approach for the Hermite-Gaussian profile with
different values of the speed ‘c’ (see Table 2). The numerical parameters and
the labels are described in Sec. 4.4 where the error is quantified by (186).

The numerical results for the convergence behaviour on $\Omega_{i}$
corresponding to the chirped-Gaussian and the Hermite-Gaussian profiles are
shown in Fig. 11 and Fig. 12, respectively. It can be seen that the diagonal
Padé approximant based method, namely, NP20 and NP50, shows stable behaviour
for the parameters specified in Table 5. In the log-log scale, we can clearly
identify the error curve to be a straight line before it plateaus. The orders
can be recovered from the slope which we found consistent with the order of
the underlying one-step method. The TR methods perform better than the BDF1
methods which is obvious from the slope of the error curves in figures 11 and
12. We note that the error background increases with increase in the value of
$c$. This is on account of the increased oscillatory behaviour of the solution
in $t$ as $c$ increases and a method with higher order than $2$ would be
required to achieve the same accuracy as that of the slow moving profiles
keeping the step-size same. Note that NP20 methods starts plateauing after a
certain number of time-steps unlike NP50. This can be attributed to the
approximate nature of the TBCs with a lower order of diagonal approximants
chosen. Finally, let us observe that there is no significant difference in the
results for convergence on $\Omega_{i}$ corresponding to the chirped-Gaussian
profile and the Hermite-Gaussian profile.

##  5 Conclusion

In this paper, we presented a new effectively local approximation for the
numerical realization of the TBCs on a rectangular computational domain for
the free Schrödinger equation in $\mathbb{R}^{2}$. Our starting point was a
numerically tractable formulation of the transparent boundary operator,
$\sqrt{\partial_{t}-i\triangle_{\Gamma}}$, where $\triangle_{\Gamma}$ is the
Laplace-Beltrami operator [4, 5], which is nonlocal in time as well as space.
This formulation of the TBCs introduces an auxiliary function which happens to
satisfy the one-dimensional Schrödinger equation along the boundary segments
of the computational domain. It addressed the nonlocality in space, however,
the nonlocality in time remained a serious bottleneck on account of the
growing computational complexity with number of time-steps which becomes
prohibitive after a point. To overcome this challenge, we developed a novel
Padé approximant based rational approximation for the $1/2$-order time-
fractional derivative. This approach circumvents the need to store the entire
history of the field on the boundary segments and allows us to achieve a
computational complexity that remains independent of the number of time-steps.
Let us emphasize that our approach is significantly different from that of
Menza [7] which fails at the corners of the rectangular domain. Further, we
also present a convolution quadrature based scheme which happens to be
computationally expensive but it is well-known as a golden standard for 1D
problems. As far as time marching is concerned, we note that one-step methods
are best suited for our purpose and any higher-order method can be readily
used to implement our algorithm. These aspects will be further explored in a
future publication.

For the spatial discretization, we use a Legendre-Galerkin spectral method
with a new boundary adapted basis which leads to a banded linear system for
the 1D problems on the boundary segments as well as for the interior problem.
The time-discrete version of the boundary condition are formulated as a Robin-
type condition in a manner reminiscent of the 1D problem [3]. A compatible
boundary-lifting procedure to homogenize the Robin-type boundary conditions is
also presented. In order to demonstrate the effectiveness of the methods
developed, we carried out numerical tests to study the behaviour of evolution
of error as well as to verify the order of convergence empirically. If the
number of spatial grid points in every dimension is
$\mathop{\mathscr{O}}\left(N\right)$, then the computational complexity of our
novel Padé method for the DtN map scales as
$\mathop{\mathscr{O}}\left(MN\right)$ where $M$ is the order of the diagonal
Padé approximant. Let us emphasize that this complexity is independent of the
number of time-steps. The overall complexity thus becomes dominated by the
cost of the linear solver for interior problem which is given by
$\Theta_{2}(N_{\text{dim}})=\mathop{\mathscr{O}}\left(N^{3}\right)$. A
theoretical proof of the stability and convergence of the numerical scheme is
beyond the scope of the paper. These issues are deferred to a future
publication.

Figure 13: The figure shows contour plots of $4\log_{10}|u(\boldsymbol{x},t)|$
with various TBCs for the chirped-Gaussian profile (see Table 1) for $c=12$ at
different instants of time. The numerical parameters and the labels are
described in Sec. 4.3. Figure 14: The figure shows contour plots of
$8\log_{10}|u(\boldsymbol{x},t)|$ with various TBCs for the Hermite-Gaussian
profile (see Table 2) for $c=8$ at different instants of time. The numerical
parameters and the labels are described in Sec. 4.3.

The ideas presented in this paper can be readily extended to three-dimensions.
Feschenko and Popov have already shown that the ideas in 2D [4] can be
extended to 3D [25]. The formulation of the operator
$\sqrt{\partial_{t}-i\triangle_{\Gamma}}$ on the boundary of a cuboidal
computational domain can be obtained along the same lines as that presented in
[5]. For all the cases discussed above, it remains to determine how
competitive are PMLs. Such a comparison should however be made for a more
general type of Schrödinger equation which is identical to the free
Schrödinger on the exterior domain
$\Omega_{e}=\mathbb{R}^{d}\setminus\overline{\Omega}_{i},\;d=2,3$.

## Acknowledgements

The first author thanks CSIR India (grant no. 09/086(1431)/2019-EMR-I) for
providing the financial assistance.

## References

  * [1] Y. S. Kivshar, G. P. Agrwal, Optical Solitons: From Fibers to Photonic Crystals, 1st Edition, Academic Press, San Diego, California, 2003.
  * [2] D. Lee, S. T. McDaniel, Ocean Acoustic Propagation by Finite Difference Methods, Modern Applied Mathematics and Computer Science, Pergmon Press, New York, 1988.
  * [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, Comm. Comput. Phys. 4 (4) (2008) 729–796.
URL https://hal.archives-ouvertes.fr/hal-00347884

  * [4] R. M. Feshchenko, A. V. Popov, Exact transparent boundary condition for the parabolic equation in a rectangular computational domain, J. Opt. Soc. Am. A 28 (3) (2011) 373–380.  doi:10.1364/JOSAA.28.000373.
  * [5] V. Vaibhav, On the nonreflecting boundary operators for the general two dimensional Schrödinger equation, J. Math. Phys. 60 (1) (2019) 011509.  doi:10.1063/1.5030875.
  * [6] L. D. Menza, Absorbing boundary conditions on a hypersurface for the Schrödinger equation in a half-space, Appl. Math. Lett. 9 (4) (1996) 55–59.  doi:10.1016/0893-9659(96)00051-1.
  * [7] L. D. Menza, Transparent and absorbing boundary conditions for the Schrödinger equation in a bounded domain, Numer. Funct. Anal. Optim. 18 (7) (1997) 759–775.  doi:10.1080/01630569708816790.
  * [8] A. Schädle, Non-reflecting boundary conditions for the two-dimensional Schrödinger equation, Wave Motion 35 (2) (2002) 181–188.  doi:10.1016/S0165-2125(01)00098-1.
  * [9] H. Han, Z. Huang, Exact artificial boundary conditions for Schrödinger equation in $\mathbb{R}^{2}$, Comm. Math. Sci. 2 (1) (2004) 79–94.
URL https://projecteuclid.org/euclid.cms/1250880210

  * [10] S. Ji, Y. Yang, G. Pang, X. Antoine, Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains, Comput. Phys. Commun. 222 (2018) 84–93.  doi:10.1016/j.cpc.2017.09.019.
  * [11] X. Antoine, C. Besse, P. Klein, Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential. Part I: Construction and a priori estimates, Math. Models Methods Appl. Sci. 22 (10) (2012) 1250026.  doi:10.1142/S0218202512500261.
  * [12] X. Antoine, C. Besse, P. Klein, Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential. Part II: Discretization and numerical results, Numer. Math. 125 (2) (2013) 191–223.  doi:10.1007/S00211-013-0542-8.
  * [13] V. Vaibhav, Microlocal approach towards construction of nonreflecting boundary conditions, J. Comput. Phys. 272 (2014) 588–607.  doi:10.1016/j.jcp.2014.04.050.
  * [14] J. Shen, Efficient spectral-Galerkin method i. Direct solvers of second- and fourth-order equations using Legendre polynomials, SIAM J. Sci. Comput. 15 (6) (1994) 1489–1505.  doi:10.1137/0915089.
  * [15] C. Lubich, Discretized fractional calculus, SIAM J. Math. Anal. 17 (3) (1986) 704–719.  doi:10.1137/0517050.
  * [16] V. Baskakov, A. Popov, Implementation of transparent boundaries for numerical solution of the Schrödinger equation, Wave Motion 14 (2) (1991) 123–128.  doi:10.1016/0165-2125(91)90053-Q.
  * [17] B. Mayfield, Nonlocal boundary conditions for the Schrödinger equation, Ph.D. thesis, University of Rhodes Island, Providence, RI (1989).
  * [18] C. Zheng, A perfectly matched layer approach to the nonlinear Schrödinger wave equations, J. Comput. Phys. 227 (1) (2007) 537–556.  doi:10.1016/j.jcp.2007.08.004.
  * [19] 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, Communications in Nonlinear Science and Numerical Simulation 90 (2020) 105406.  doi:10.1016/j.cnsns.2020.105406.
  * [20] K. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, 1st Edition, John Wiley and Sons Inc., 1993.
  * [21] V. Vaibhav, Artificial boundary conditions for certain evolution pdes with cubic nonlinearity for noncompactly supported initial data, J. Comput. Phys. 230 (8) (2011) 3205–3229.  doi:10.1016/j.jcp.2011.01.024.
  * [22] S. Yadav, V. Vaibhav, Nonreflecting boundary condition for the free Schrödinger equation in 2d, in: 2023 Photonics and Electromagnetics Research Symposium (PIERS), 2023, pp. 328–337.  doi:10.1109/PIERS59004.2023.10221299.
  * [23] N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd Edition, Society for Industrial and Applied Mathematics, Philadelphia, 2002.  doi:10.1137/1.9780898718027.
  * [24] J. Shen, T. Tang, L. L. Wang, Spectral Methods: Algorithms, Analysis and Applications, 1st Edition, Springer Berlin Heidelberg, 2011.
  * [25] R. M. Feshchenko, A. V. Popov, Exact transparent boundary condition for the three-dimensional Schrödinger equation in a rectangular cuboid computational domain, Phys. Rev. E 88 (2013) 053308.  doi:10.1103/PhysRevE.88.053308.

##  Appendix A Fractional Operators

######  Definition 1 (Riemann-Liouville fractional integrals).

Let $\alpha>0$ and $f(t)$ be a piecewise continuous function on $(0,\infty)$,
locally integrable on any finite subinterval of $[0,\infty)$. For $t>0$, the
Riemann-Liouville fractional integral of order $\alpha$, denoted by
$\partial_{t}^{-\alpha}$, is defined by

$\partial_{t}^{-\alpha}f(t)=\frac{1}{\Gamma(\alpha)}\int^{t}_{0}(t-s)^{\alpha-1}f(s)ds.$
(187)

where $\Gamma$ denotes the Euler’s Gamma function.

######  Theorem 1 (Law of exponents).

Let $f(t)\in\mathsf{C}^{\infty}([0,\infty))$ and let $\mu,\nu>0$, then the
composition of two fractional integrals is given by

$\partial_{t}^{-\nu}[\partial_{t}^{-\mu}f(t)]=\partial_{t}^{-(\mu+\nu)}f(t)=\partial_{t}^{-\mu}[\partial_{t}^{-\nu}f(t)],\quad\forall
t\geq 0.$

######  Definition 2 (Riemann-Liouville fractional derivatives).

Let $f(t)\in\mathsf{C}^{\infty}([0,\infty))$ and let $\mu>0$. If $p$ is the
smallest integer greater than $\mu$, then the fractional derivative of $f$ of
order $\mu$ is defined by

$\partial_{t}^{\mu}f(t)=\partial^{p}_{t}[\partial_{t}^{-\nu}f(t)],\quad t>0,$
(188)

where $\nu=p-\mu>0$. For this fractional derivative to be defined at $t=0$,
the function $f(t)$ must have first $p-1$ derivatives equal to zero at $t=0$.

##  Appendix B Boundary Adapted Basis Representation: Useful Results and
Identities

In order to leverage the discrete Legendre transforms over Legendre-Gauss-
Lobatto nodes, we choose to favour storage of space dependent variables in
terms of their Legendre coefficients. As a result of this choice, all
operations involving the boundary adapted basis must be implemented via the
Legendre transform.

Let us recall the definitions: $L_{n}(y)$ denote the Legendre polynomial of
degree $n$ and the index set $\\{0,1,\ldots,N-2\\}$ is denoted by
$\mathbb{J}$. The boundary adapted basis polynomials are given by

$\phi_{p}(y)=L_{p}(y)+b_{p}L_{p+2}(y),\quad
b_{p}=-\frac{\kappa+\frac{1}{2}p(p+1)}{\kappa+\frac{1}{2}(p+2)(p+3)},\quad
p\in\mathbb{J},$ (189)

where $\kappa\in\mathbb{C}$. Let us briefly discuss the parameters that
determine this quantity: We recall from Sec. 3.2.1 that $\beta=J^{-2}$ where
$J=\partial_{y}x$ is the Jacobian corresponding to the transformation from the
reference interval $\mathbb{I}=[-1,1]$ to the actual domain $[x_{-},x_{+}]$.
We suppress the subscript which characterizes the axis along which the domain
lies for the sake of brevity. The time-stepping schemes discussed in Sec. 3
dictate that either $\rho=1/\Delta t$ (for BDF1) or $\rho=2/\Delta t$ (for
TR). The form of $\kappa$ follows from the Robin-type formulation of the
discretized version of the DtN-maps as discussed in Sec. 3.1. For CQ methods
$\kappa=\alpha=\sqrt{\rho/\beta}\exp(-i\pi/4)$ while for the NP methods, we
have $\kappa=\alpha\varpi$ where $\varpi$ is defined in (82).

Any arbitrary function in $\mathsf{L}^{2}(\mathbb{I})$ can be expanded in this
basis as follows:

$f(y)=\sum_{p=0}^{N-2}\hat{f}_{p}\phi_{p}(y)+\eta_{0}L_{0}(y)+\eta_{1}L_{1}(y)=\sum_{p^{\prime}=0}^{N}\tilde{f}_{p^{\prime}}L_{p^{\prime}}(y).$
(190)

By virtue of linear independence of Legendre polynomials, we can equate their
coefficients on both sides of the equality so that

$\left\\{\begin{aligned}
&\tilde{f}_{0}=\hat{f}_{0}+\eta_{0},\tilde{f}_{1}=\hat{f}_{1}+\eta_{1},\\\
&\tilde{f}_{p}=\hat{f}_{p}+b_{p-2}\hat{f}_{p-2},\quad p=2,3,\ldots N-2,\\\
&\tilde{f}_{N-1}=b_{N-3}\hat{f}_{N-3},\tilde{f}_{N}=b_{N-2}\hat{f}_{N-2}.\end{aligned}\right.$
(191)

The linear system can be written as

$\widetilde{\boldsymbol{f}}=B\widehat{\boldsymbol{f}}+\eta_{0}\boldsymbol{e}_{0}+\eta_{1}\boldsymbol{e}_{1}$
(192)

where

$B=\begin{pmatrix}1&&&&&\\\ 0&1&&&&\\\ b_{0}&0&1&&&\\\ &b_{1}&0&1&&\\\
&&\ddots&\ddots&\ddots&\\\ &&&b_{N-4}&0&1\\\ &&&&b_{N-3}&0\\\
&&&&&b_{N-2}\end{pmatrix}\in\mathbb{C}^{(N+1)\times(N-1)},$ (193)

Let the normalization factors be denoted by
$\gamma_{k}=\|L_{k}\|^{2}_{2}=2/(2k+1),\,k\in\mathbb{N}_{0},$ and introduce
the diagonal matrix
$\Gamma=\operatorname{diag}(\gamma_{0},\gamma_{1},\ldots,\gamma_{N})$, then
inner products of the form
$g_{p}=\left(f,\phi_{p}\right)_{\mathbb{I}},\,p\in\mathbb{J},$ can be
understood as

$g_{p}=\left(f,\phi_{p}\right)_{\mathbb{I}}=\left[\gamma_{p}\tilde{f}_{p}+b_{p}\gamma_{p+2}\tilde{f}_{p+2}\right],\quad
p\in\mathbb{J}.$ (194)

The relationship above can be compactly stated as
$\boldsymbol{g}=(g_{0},g_{1},\ldots,g_{N-2})^{\intercal}=Q\Gamma\widetilde{\boldsymbol{f}}$
where $Q=B^{\intercal}$ is referred to as the _quadrature matrix_.

Let the index set $\\{0,1,\ldots,N_{1}-2\\}$ be denoted by $\mathbb{J}_{1}$
and the index set $\\{0,1,\ldots,N_{2}-2\\}$ be denoted by $\mathbb{J}_{2}$.
The boundary adapted basis polynomials are given by

$\phi^{(j)}_{p_{j}}(y_{j})=L_{p_{j}}(y_{j})+b^{(j)}_{p_{j}}L_{p_{j}+2}(y_{j}),\quad
p_{j}\in\mathbb{J}_{j},\;j=1,2.$ (195)

where $\\{b^{(j)}_{p_{j}}\\}$ is a sequence. Any arbitrary function in
$\mathsf{L}^{2}(\Omega_{i}^{\text{ref}})$ can be expanded in this basis as
follows:

$\begin{split}f(\boldsymbol{y})&=\sum_{p_{1}\in\mathbb{J}_{1}}\sum_{p_{2}\in\mathbb{J}_{2}}\hat{f}_{p_{1},p_{2}}\phi^{(1)}_{p_{1}}(y_{1})\phi^{(2)}_{p_{2}}(y_{2})\\\
&\quad+\sum_{p_{1}\in\mathbb{J}_{1}}\left[\zeta_{p_{1},0}\phi^{(1)}_{p_{1}}(y_{1})L_{0}(y_{2})+\zeta_{p_{1},1}\phi^{(1)}_{p_{1}}(y_{1})L_{1}(y_{2})\right]\\\
&\quad+\sum_{p_{2}\in\mathbb{J}_{2}}\left[\eta_{0,p_{2}}L_{0}(y_{1})\phi^{(2)}_{p_{2}}(y_{2})+\eta_{1,p_{2}}L_{1}(y_{1})\phi^{(2)}_{p_{2}}(y_{2})\right]\\\
&\quad+\left[\kappa_{0,0}L_{0}(y_{1})L_{0}(y_{2})+\kappa_{0,1}L_{0}(y_{1})L_{1}(y_{2})+\kappa_{1,0}L_{1}(y_{1})L_{0}(y_{2})+\kappa_{1,1}L_{1}(y_{1})L_{1}(y_{2})\right]\\\
&=\sum_{p_{1}=0}^{N_{1}}\sum_{p_{2}=0}^{N_{2}}\tilde{f}_{p_{1},p_{2}}L_{p_{1}}(y_{1})L_{p_{2}}(y_{2}).\end{split}$
(196)

For $p_{1}=2,3,\ldots N_{1}-2$ and $p_{2}=2,3,\ldots N_{2}-2$, equating the
coefficients on both sides of equality so that

$\tilde{f}_{p_{1},p_{2}}=\hat{f}_{p_{1},p_{2}}+b^{(1)}_{p_{1}-2}\hat{f}_{p_{1}-2,p_{2}}+b^{(2)}_{p_{2}-2}\hat{f}_{p_{1},p_{2}-2}+b^{(1)}_{p_{1}-2}b^{(2)}_{p_{2}-2}\hat{f}_{p_{1}-2,p_{2}-2},$
(197)

The linear system can be written as

$\begin{split}\widetilde{F}&=B_{1}\widehat{F}B_{2}^{\intercal}+B_{1}\boldsymbol{\zeta}_{0}\otimes\left(\boldsymbol{e}^{(2)}_{0}\right)^{\intercal}+B_{1}\boldsymbol{\zeta}_{1}\otimes\left(\boldsymbol{e}^{(2)}_{1}\right)^{\intercal}+\boldsymbol{e}^{(1)}_{0}\otimes\boldsymbol{\eta}_{0}B_{2}^{\intercal}+\boldsymbol{e}^{(1)}_{1}\otimes\boldsymbol{\eta}_{1}B_{2}^{\intercal}\\\
&\quad+\kappa_{0,0}\boldsymbol{e}^{(1)}_{0}\otimes\left(\boldsymbol{e}^{(2)}_{0}\right)^{\intercal}+\kappa_{0,1}\boldsymbol{e}^{(1)}_{0}\otimes\left(\boldsymbol{e}^{(2)}_{1}\right)^{\intercal}+\kappa_{1,0}\boldsymbol{e}^{(1)}_{1}\otimes\left(\boldsymbol{e}^{(2)}_{0}\right)^{\intercal}+\kappa_{1,1}\boldsymbol{e}^{(1)}_{1}\otimes\left(\boldsymbol{e}^{(2)}_{1}\right)^{\intercal}\end{split}$
(198)

where

$B_{j}=\begin{pmatrix}1&&&&&\\\ 0&1&&&&\\\ b^{(j)}_{0}&0&1&&&\\\
&b^{(j)}_{1}&0&1&&\\\ &&\ddots&\ddots&\ddots&\\\ &&&b^{(j)}_{N-4}&0&1\\\
&&&&b^{(j)}_{N-3}&0\\\
&&&&&b^{(j)}_{N-2}\end{pmatrix}\in\mathbb{C}^{(N+1)\times(N-1)},$ (199)

The inner products of the form
$g_{p_{1},p_{2}}=\left(f,\phi^{(1)}_{p_{1}}\phi^{(2)}_{p_{2}}\right)_{\Omega_{i}^{\text{ref}}}$
can be understood as

$\begin{split}g_{p_{1},p_{2}}=\left(f,\phi^{(1)}_{p_{1}}\phi^{(2)}_{p_{2}}\right)_{\Omega_{i}^{\text{ref}}}&=\gamma_{p_{1}}\gamma_{p_{2}}\tilde{f}_{p_{1},p_{2}}+b^{(1)}_{p_{1}}\gamma_{p_{1}+2}\gamma_{p_{2}}\tilde{f}_{p_{1}+2,p_{2}}\\\
&\qquad+b^{(2)}_{p_{2}}\gamma_{p_{1}}\gamma_{p_{2}+2}\tilde{f}_{p_{1},p_{2}+2}+b^{(1)}_{p_{1}}b^{(2)}_{p_{2}}\gamma_{p_{1}+2}\gamma_{p_{2}+2}\tilde{f}_{p_{1}+2,p_{2}+2}.\end{split}$
(200)

The relationship above can be compactly stated as
$G=(g_{p_{1},p_{2}})=Q_{1}\Gamma\widetilde{F}\Gamma Q_{2}^{\intercal}$ where
$Q_{j}=B_{j}^{\intercal}$ are referred to as the _quadrature matrices_.