Absorbing boundary conditions for the time-dependent Schr\"odinger-type equations in $\mathbb R^3$

TL;DR

This paper develops and analyzes absorbing boundary conditions for 3D time-dependent Schrödinger equations, improving computational efficiency and accuracy in quantum simulations with arbitrary geometries.

Contribution

It introduces a semi-discrete derived boundary condition compatible with finite-difference schemes and provides stability proofs and practical tests.

Findings

Boundary conditions are stable for zeroth and first order.

Effective in reducing boundary reflections in benchmark problems.

Accurate in energy and nucleon density for Hartree-Fock models.

Abstract

Absorbing boundary conditions are presented for three-dimensional time-dependent Schr\"odinger-type of equations as a means to reduce the cost of the quantum-mechanical calculations. The boundary condition is first derived from a semi-discrete approximation of the Schr\"odinger equation with the advantage that the resulting formulas are automatically compatible with the finite-difference scheme and no further discretization is needed in space. The absorbing boundary condition is expressed as a discrete Dirichlet-to-Neumann (DtN) map, which can be further approximated in time by using rational approximations of the Laplace transform to enable a more efficient implementation. This approach can be applied to domains with arbitrary geometry. The stability of the zeroth order and first order absorbing boundary conditions is proved. We tested the boundary conditions on benchmark problems. The effectiveness is further verified by a time-dependent Hartree-Fock model with Skyrme interactions. The accuracy in terms of energy and nucleon density is examined as well.

Equations

$$i\frac{\partial }{\partial t}\varphi (\mathbf{x},t)=\hat{H}\varphi (\mathbf{x},t),\mathbf{x}\in \Omega .$$

$$\frac{{\partial }^2}{\partial x^2}\psi (x,t)\approx \frac{-{\psi }_{-2}+16{\psi }_{-1}-30{\psi }_0+16{\psi }_1-{\psi }_2}{12h^2}.$$

$$\{\begin{array}{rl}i{\dot{\varphi }}_{\text{I}}(t)& =H_{\text{I,I}}{\varphi }_{\text{I}}(t)+H_{\text{I,II}}{\varphi }_{\text{II}}(t),\\ i{\dot{\varphi }}_{\text{II}}(t)& =H_{\text{II,I}}{\varphi }_{\text{I}}(t)+H_{\text{II,II}}{\varphi }_{\text{II}}(t),\end{array}$$

$$H=\left[\begin{array}{cc}{H}_{\text{I,I}}& H_{\text{I,II}}\\ H_{\text{II,I}}& H_{\text{II,II}}\end{array}\right].$$

$$\Gamma =\{{\mathbf{x}}_j\in {\Omega }_{\text{I}}|\text{ if there exists }{\mathbf{x}}_k\in {\Omega }_{\text{II}}\text{ such that }{H}_{jk}\ne 0\}.$$

$${\mathbf{\varphi }}_{\text{I}}=\left[\begin{array}{c}{\mathbf{\varphi }}_{\Gamma }\\ {\mathbf{\varphi }}_{\text{I}\backslash \Gamma }\end{array}\right].$$

$$H_{\text{I\!I,I}}=\big{[}H_{\text{I\!I,}\Gamma}\;\;H_{\text{I\!I},\text{I}\backslash\Gamma}\big{]}=\big{[}H_{\text{I\!I,}\Gamma}\;\;0\big{]}.$$

$${\mathbf{\varphi }}_{\Gamma }=E{\mathbf{\varphi }}_{\text{I}}.$$

$$H_{\text{I,II}}=E^T{H}_{\Gamma ,\text{II}}.$$

$${\mathbf{f}}_{\Gamma }=H_{\Gamma ,\text{II}}{\mathbf{\varphi }}_{\text{II}}.$$

$${\dot{\mathbf{\varphi }}}_{\text{I}}(t)=-i{H}_{\text{I,I}}{\mathbf{\varphi }}_{\text{I}}(t)-i{E}^T{\mathbf{f}}_{\Gamma }(t).$$

$$is{\mathbf{\Phi }}_{\text{II}}(s)=H_{\text{II,II}}{\mathbf{\Phi }}_{\text{II}}(s)+H_{\text{II,I}}{\mathbf{\Phi }}_{\text{I}}(s),$$

$${\tilde{H}}_{\text{II,II}}(s){\mathbf{\Phi }}_{\text{II}}(s)=-H_{\text{II,I}}{\mathbf{\Phi }}_{\text{I}}(s).$$

$${\mathbf{\Phi }}_{\text{II}}(s)=-{\tilde{H}}_{\text{II,II}}^{-1}{H}_{\text{II,I}}{\mathbf{\Phi }}_{\text{I}}(s).$$

$$H_{\text{II,I}}{\mathbf{\Phi }}_{\text{I}}=H_{\text{II,}\Gamma }E{\mathbf{\Phi }}_{\text{I}}=H_{\text{II,}\Gamma }{\mathbf{\Phi }}_{\Gamma }.$$

$${\mathbf{\Phi }}_{\text{II}}(s)=-{\tilde{H}}_{\text{II,II}}^{-1}{H}_{\text{II,}\Gamma }{\mathbf{\Phi }}_{\Gamma }(s).$$

$${\mathbf{F}}_{\Gamma }(s)=K(s){\mathbf{\Phi }}_{\Gamma }(s).$$

$$K(s)=-H_{\Gamma,\text{I\!I}}\widetilde{H}_{\text{I\!I,I\!I}}^{-1}(s)H_{\text{I\!I,}\Gamma}=-H_{\Gamma,\text{I\!I}}\big{[}H_{\text{I\!I,I\!I}}-isI\big{]}^{-1}H_{\text{I\!I,}\Gamma}.$$

$${\mathbf{f}}_{\Gamma }(t)={\int }_0^t\kappa (t-\tau ){\mathbf{\varphi }}_{\Gamma }(\tau )d\tau .$$

$${\dot{\mathbf{\varphi }}}_{\text{I}}(t)=-i{H}_{\text{I,I}}{\mathbf{\varphi }}_{\text{I}}(t)-i{E}^T{\int }_0^t\kappa (t-\tau )E{\mathbf{\varphi }}_{\text{I}}(\tau )d\tau ,$$

$$i(s-s_0)(H_{\text{II,II}}-i{s}_{0}I{)}^{-1}{\mathbf{\Phi }}_{\text{II}}(s)={\mathbf{\Phi }}_{\text{II}}(s)+(H_{\text{II,II}}-i{s}_{0}I{)}^{-1}{H}_{\text{II,I}}{\mathbf{\Phi }}_{\text{I}}(s),$$

$$K(s)=-H_{\Gamma ,\text{II}}(I-i(s-s_0)(H_{\text{II,II}}-i{s}_{0}I{)}^{-1}{)}^{-1}(H_{\text{II,II}}-i{s}_{0}I{)}^{-1}{H}_{\text{II,}\Gamma }.$$

$$K(s)=-H_{\Gamma ,\text{II}}(C+\mathcal{A}C(s-s_0)+{\mathcal{A}}^{2}C(s-s_0{)}^2+{\mathcal{A}}^{3}C(s-s_0{)}^3+\cdots ).$$

$$R_{m,m}(s)=(s^m-s^{m-1}{B}_0-\cdots -B_{m-1}{)}^{-1}(s^{m-1}{A}_0+\cdots +A_{m-1}).$$

$${\mathbf{F}}_{\Gamma }(s)=R_{m,m}(s){\mathbf{\Phi }}_{\Gamma }(s).$$

$${\mathbf{f}}_{\Gamma }^{(m)}=B_0{\mathbf{f}}_{\Gamma }^{(m-1)}+\cdots +B_{m-1}{\mathbf{f}}_{\Gamma }+A_0{\mathbf{\varphi }}_{\Gamma }^{(m-1)}+\cdots +A_{m-1}{\mathbf{\varphi }}_{\Gamma },$$

$$\{\begin{array}{rl}& \frac{\partial }{\partial t}{\mathbf{\varphi }}_{\text{I}}(t)=-i{H}_{\text{I,I}}{\mathbf{\varphi }}_{\text{I}}(t)-i{H}_{\text{I,II}}{\mathbf{\varphi }}_{\text{II}}(t),\\ & {\mathbf{f}}_{\Gamma }^{(m)}=B_0{\mathbf{f}}_{\Gamma }^{(m-1)}+\cdots +B_{m-1}{\mathbf{f}}_{\Gamma }+A_0{\mathbf{\varphi }}_{\Gamma }^{(m-1)}+\cdots +A_{m-1}{\mathbf{\varphi }}_{\Gamma }.\end{array}$$

$${\dot{\mathbf{\varphi }}}_{\text{I}}(t)=-i{H}_{\text{I,I}}{\mathbf{\varphi }}_{\text{I}}(t)-i{E}^{T}ME{\mathbf{\varphi }}_{\text{I}}(t).$$

$$W(t)={\mathbf{\varphi }}_{\text{I}}^{\ast }(t){\mathbf{\varphi }}_{\text{I}}(t).$$

$$\frac{dW}{dt}$$

Full text

I Introduction

Quantum-mechanical simulations, expressed in terms of the Schrödinger equation, provide a fundamental description of physical properties in chemistry and condensed-matter physics martin2004 ; helgaker2014 . On the other hand, due to the small length scale, an outstanding challenge is the size of the system that one can simulate, even with the rapid growth of computing power. Computer simulations often face the scenarios where the electrons are being emitted out of the computational domain, e.g., the photoionization process. Recently there has been a great deal of renewed interest in such issues in the quantum transport problem xie2014complex , optical response of molecules yabana2006real , open quantum systems Hellumns1994 , etc.

One approach to address the aforementioned issue is by absorbing boundary conditions (ABCs) which reduce the problem to a much smaller computational domain and yield results as if the simulation is being performed in a much larger (or unbounded) domain. Rather than simply removing the exterior region, the ABCs provide an efficient approach to mimic the influence from the surrounding environment. There are several different approaches to construct and implement ABCs, most of which involve the derivation and approximation of the Dirichlet-to-Neumann (DtN) map. They are also commonly referred to as non-reflecting, transparent, or radiating boundary conditions. By replacing the surrounding region with the ABCs, the computational effort can be focused on simulating the region of interest.

There seemed to be separate development of ABCs for quantum-mechanical simulations. For example, there is a great deal of progress in developing ABCs in physics and chemistry. These ABCs can be roughly classified into four approaches: 1) the exterior complex scaling (ECS) Scrinzi2014 ; NICOLAIDES197811 ; SIMON1979211 ; mccurdy2004solving , where the coordinate outside a fixed radius is scaled into the complex plane; 2) the mask function method marques2003 ; tafipolsky2006 ; Wang201107 , where the wave functions are gradually scaled to zero; 3) complex absorbing potentials (CAP) riss1996investigation ; Riss_1993 ; MUGA2004357 , which is introduced outside the boundary by adding a complex potential to the Hamiltonian; and 4) coordinate scaling, which scales the spatial coordinate by a time-dependent factor Soloviev1985 ; sidky2000 to make wave functions compatible with the physical boundary conditions. These methods are not constructed to directly approximate the exact boundary condition. Rather, they are very much goal-oriented. Their direct aim is to absorb the electrons that move outside the computational domain. These methods in general are easy to implement, and are quite robust in practice. Meanwhile, the effectiveness depends heavily on the choice of the parameters, e.g., in the CAP methods, the size of the boundary region and the magnitude of the potential. In some cases, the efficiency and accuracy (reflections) are even debatable for decades Riss_1998 ; He2007 ; Tao2009 ; Scrinzi2010 ; Ge1998 . The ABCs can also be introduced for the density-matrix in the Liouville von-Neumann equation li2019absorbing .

In the applied math community, much effort has been focused on one-dimensional problems, with a few works for multi-dimensional problems with flat boundary, for a comprehensive review, see antoine2008review ; han2013 . There, a different path is followed: The exact ABCs are first derived by solving an exterior problem, many of which are written as a Dirichlet-to-Neumann (DtN) or Neumann-to-Dirichlet (NtD) map for the continuous time-dependent Schrödinger equation (TDSE) han2005finite (or temporally discrete model antonie2003 , spatially discrete model alonso2002 ; alonso2003 , and fully discrete model arnold1998 ). The same technique has been used by Ermolaev et al. ermolaev1999 ; ermolaev2000 to treat electron dynamics under a laser field, where the exact boundary condition, referred to as an integral boundary condition, is derived with the help of the Green’s function for a free electron. Then, the exact ABCs, which involve a time-convolution, are approximated to avoid the repeated calculations of the integrals. The most predominant method is by finite sums of exponentials in the real-time space arnold2003 ; Jiang2004 or rational functions in the Laplace domain or the Fourier space fevens1999 ; shibata1991 ; szeftel2004 ; xu2006absorbing to facilitate a fast evaluation. The convergence, stability, and efficiency of these methods have been thoroughly investigated antoine2008review ; han2013 . The most extensively studied case is the one-dimensional case. The extensions to high dimensional cases can be found in jiang2008 ; antoine2004 ; Alpert2002 ; han2004 ; Xu2007 ; but these methods are often limited only to special geometries of the interior domain. Another classical strategy is the perfectly matched layer (PML) berenger1994 , where one first constructs a buffer layer so that the outgoing waves in the computational domain are exactly preserved (perfect matching). The most common approach is to introduce a complex stretching of the spatial coordinate to derive a modified equation in the buffer zone, and then the resulting models are discretized simultaneously in the implementation. The PML has been applied to the nonlinear Schrödinger equation by Zheng zheng2007 . Similar to the CAP methods, the parameters in the PML method have to be calibrated a priori.

In the three-dimensional case, one challenge in implementing the ABCs stems from the fact that the kinetic energy in the Schrödinger equation often needs to be approximated by high-order finite-difference schemes, either to gain sufficient accuracy, or to reduce the computational cost by using larger grid size. Therefore, the ABCs that are developed from continuous PDEs, especially those that are expressed in terms of the normal derivatives of the solution at the boundary bayliss1980radiation ; higdon1986absorbing , are difficult to be integrated with the discretization in the computational domain. For instance, due to the large finite-difference stencil, there are several layers of grids at the boundary, where the ABCs need to be applied, which is non-trivial. To overcome this difficulty, we formulate the ABCs based on a semi-discrete approximation of the TDSE obtained from high-order finite-difference schemes in the entire domain. The resulting boundary condition does not need to be further discretized in space. It can be readily implemented. Meanwhile, the time-convolution in the corresponding DtN map is formulated and approximated by rational functions in Laplace domain. This approximation reduces the exact DtN map to ordinary differential equations in the time domain, which is much more efficient in practice. We also tackle the important issue that when the domains is of arbitrary geometry, there is no simple representations of the DtN map, e.g., in terms of a pseudo-differential operator antonie2004 ; szeftel2004 . In this case, it is expressed as a matrix-function which involves the Hamiltonian in the exterior. We will employ the discrete boundary-element method Li2012 to evaluate the DtN map in the rational interpolation.

We will demonstrate that the first-order ABC obtained from this approach correspond to a CAP method. This offers an interesting connection between the DtN map and an existing method. But the Hamiltonian that represents the absorbing potential is computed from the exact ABC, rather than empirically built (e.g., a diagonal matrix). We will also show that the second-order ABC is more general than CAP. In most cases, it is more accurate than the first-order approximation.

This paper is organized as follows: In section II, the solution of the exterior problem, which leads to the exact ABC, is formulated. We approximate the exact boundary condition based on its Laplace transform in section III. We discuss the choices of the interpolation conditions such that the dynamics is stable. In section IV, both 1D and 3D numerical simulations are performed to verify the stability and effectiveness. Further, we test the ABCs on a time-dependent Hartree-Fock (TDHF) model with localized interactions.

II Formulation of a discrete DtN map for Schrödinger-type of equations

This paper is concerned with the time-dependent Schrödinger equations in a domain $\Omega\subset\mathbb{R}^{d}$,

[TABLE]

Here we have scaled the Planck constant and mass to unity. $i$ is the imaginary unit. $\phi(\bm{x},t)$ is the wave function for an electron or nucleon in the $d$-dimensional quantum system. In general, the Hamiltonian operator can be written as $\hat{H}=-\frac{\nabla^{2}}{2}+V(\bm{x},t)$, with $V$ being the total potential of the system. A Schrödinger-like equation can be derived from the exact many-body Schrödinger equation with various mean-field approximations of the many-body wave function, such as the Hartree-Fock approach baerends1973self and the time-dependent density-functional theory marques2006time . For simplicity, we illustrate the formulation with the simplest model, and demonstrate the applications of other Schrödinger-type equations in section IV.3.

We suppose that $\Omega_{\text{I}}$ is a subdomain of interest in $\Omega$. Rather than solving the Schrödinger equations in the entire domain $\Omega,$ we take into account the influence from the surrounding region, here denoted by $\Omega_{\text{I\!I}}$, by deriving ABCs at the boundary of $\Omega_{\text{I}}$. $\Omega=\Omega_{\text{I}}\cup\Omega_{\text{I\!I}}$. In the exterior region $\Omega_{\text{I\!I}}$, we assume that the operator $\hat{H}$ does not explicitly depend on $\bm{x}$. More specifically, $V(\bm{x})=0$, or more generally a constant, for any $\bm{x}\in\Omega_{\text{I\!I}}$.

Most ABCs have been formulated in the continuous case szeftel2004 ; antonie2003 ; baskakov1991implementation ; Hellumns1994 , which has to be further discretized in order to be combined with the finite-difference or finite-element approximation in the interior of $\Omega_{\text{I}}$. For easier implementation, our approach starts directly with a spatially discrete model. Namely, we first consider the semi-discrete approximation of (1) in the entire domain, and regard it as a full/exact model. This has been done for the acoustic wave equations in givoli1998discrete . But efficient evaluations of the ABCs were not fully discussed there.

Let $\{\bm{x}_{j}\}$ be the set of $n_{\text{I}}$ grid points with a constant spacing of $h$ in the interior $\Omega_{\text{I}}$ and $n_{\text{I\!I}}$ grid points in the exterior $\Omega_{\text{I\!I}}$. $n_{\text{I}}\ll n_{\text{I\!I}}$. $n_{\text{I\!I}}$ can be infinite. In real-space methods, the Hamiltonian operator in (1) can be approximated by a representation over numerical grids (e.g., octopus ), finite elements (e.g., motamarri2013higher ), or atom centered basis set. These spatial approximations reduce the problem to a discrete model in space. Here we take the simplest approach, the finite difference method, to demonstrate the formulation of the ABC. To ensure high order accuracy, these finite difference methods often employ large stencils. For example, a one dimensional Laplacian operator, acting on a continuous function $\psi(x,t)$, can be approximated by,

[TABLE]

Here $\psi_{j}(t)=\psi(x+jh,t)$ with $h$ being the grid spacing. In the high-dimensional case, this finite-difference formula can be applied in each direction beck2000real .

Let us introduce some notations to make the derivation more transparent. These notations are largely adopted from domain decomposition methods quarteroni1999domain . To begin with, the indices of the nodal points are sorted as $\bm{x}_{1},\dots,\bm{x}_{n_{\text{I}}}\in\Omega_{\text{I}}$, and $\bm{x}_{n_{\text{I}}+1},\dots,\bm{x}_{n_{\text{I}}+n_{\text{I\!I}}}\in\Omega_{\text{I\!I}}$. Due to the discretization in space, the wave functions are written in vectors, and the semi-discrete model can be written in a compact form as follows:

[TABLE]

where $\bm{\phi}_{\text{I}}=[\phi(\bm{x}_{k})]_{\bm{x}_{k}\in\Omega_{\text{I}}}$ and $\bm{\phi}_{\text{I\!I}}=[\phi(\bm{x}_{k})]_{\bm{x}_{k}\in\Omega_{\text{I\!I}}}$. $\bm{\phi}_{\text{I}}\in\mathbb{C}^{n_{\text{I}}}$ and $\bm{\phi}_{\text{I\!I}}\in\mathbb{C}^{n_{\text{I\!I}}}$. $\dot{}$ means the derivative with respect to the time $t$. Following the same partition, the discretization of the Hamiltonian operator $H$ is structured as follows,

[TABLE]

Since the Schrödinger equation (1) has been discretized in space, we use $\dot{\quad}$ to denote the time derivatives hereafter.

Our goal is to simplify the second equation in (2), while retaining the first equation. $H_{\text{I,I}}$ and $H_{\text{I\!I,I\!I}}$ are respectively the discretized Hamiltonian operator in $\Omega_{\text{I}}$ and $\Omega_{\text{I\!I}}$. $H_{\text{I,I}}$ can be a nonlinear operator. $H_{\text{I,I\!I}}$ and $H_{\text{I\!I,I}}$ are the off-diagonal blocks representing the coupling between $\Omega_{\text{I}}$ and $\Omega_{\text{I\!I}}$. Notice that due to the finite-difference approximation of the kinetic energy term, all these matrices are sparse.

The boundary condition will be expressed in terms of the values of the wave functions near the boundary of $\Omega_{\text{I}}$. To this end, we denote $\Gamma$ as the boundary of $\Omega_{\text{I}}$, such that

[TABLE]

The width of $\Gamma$ depends on the width of the finite-difference stencil. We define $\bm{\phi}_{\Gamma}\in\mathbb{C}^{n_{\Gamma}}$ as a vector formed by all $\bm{\phi}_{j}$ for $\bm{x}_{j}\in\Gamma$, where $n_{\Gamma}$ is the number of grid points in $\Gamma$. The vector $\bm{\phi}_{\text{I}}$ can be reordered so that the first $n_{\Gamma}$ components are associated with the grid points in $\Gamma$, i.e.,

[TABLE]

The boundary region $\Gamma$ is defined in such a way that there is no direct coupling between the points in the interior of $\Omega_{\text{I}}$ and those in $\Omega_{\text{I\!I}}$. This is reflected in the off-diagonal block of the Hamiltonian,

[TABLE]

In general, this procedure can be carried out by defining a restriction operator $E$ to extract the components of a function that correspond to grid points at the boundary $\Gamma$ from a function defined in $\Omega_{\text{I}}$. Namely,

[TABLE]

With the reordering in (5), the matrix $E\in\mathbb{R}^{n_{\Gamma}\times n_{\text{I\!I}}}$ can be explicitly expressed as $E=[I_{n_{\Gamma}}\;\;0]$. Further, (6) can now be simply written as $H_{\text{I\!I,}\Gamma}=H_{\text{I\!I,I}}E^{T}$. Since $EE^{T}=I_{n_{\Gamma}}$, we also have,

[TABLE]

We are now set to formulate the ABC. We first consider the influence of the wave functions in $\Omega_{\text{I\!I}}$ on the wave functions in $\Omega_{\text{I}}$, and we define

[TABLE]

Then the first equation reads,

[TABLE]

Here $\dot{}$ denotes the derivative with respect to the time variable $t$.

At this point, we observe that what is needed to compute $\bm{\phi}_{\text{I}}$ in time is $\bm{f}_{\Gamma}(t).$ For this purpose, let us take the Laplace transform of the second equation in (2). Due to the assumption that $\Phi_{\text{I\!I}}(0)=0,$ one has,

[TABLE]

where $\bm{\Phi}_{\text{I\!I}}(s)=\mathcal{L}\{\bm{\phi}_{\text{I\!I}}\}(s)$ and $\bm{\Phi}_{\text{I}}(s)=\mathcal{L}\{\bm{\phi}_{\text{I}}\}(s)$. By defining $\widetilde{H}_{\text{I\!I,I\!I}}(s)=H_{\text{I\!I,I\!I}}-isI$, equation (11) simply reads

[TABLE]

Formally, its solution can be expressed as

[TABLE]

Due to the locality of $H_{\text{I\!I,I\!I}}$, the right hand side of (13) can be further reduced. Notice that,

[TABLE]

The first step is from the identity (6), and the second step used (7). This simplifies the solution to,

[TABLE]

Since $H_{\text{I\!I,I\!I}}$ is a $n_{\text{I\!I}}\times n_{\text{I\!I}}$ matrix, it is in general impractical to solve the linear system (14) numerically. However, the main observation from (10) is that we only need $\bm{f}_{\Gamma}$ to keep the computation in $\Omega_{\text{I}}$. We let $\bm{F}_{\Gamma}$ be the Laplace transform of $\bm{f}_{\Gamma}$. By left-multiplying the equation above by $H_{\Gamma,\text{I\!I}}$, we obtain,

[TABLE]

The matrix-valued function $K(s):\mathbb{R}\to\Omega_{\Gamma}\times\Omega_{\Gamma}$, which will play a key role in the numerical approximation, is given by,

[TABLE]

The mapping $K(s)$ is precisely the Dirichlet-to-Neumann (DtN) map of the problem (2) in the Laplace domain. Although $K(s)$ still involves the inverse of a large matrix, it is much easier to compute than (14) due to the fact that only a small number of entries are needed. This is known as selected inversion linalg . The detailed solution to this problem, which involves a discrete boundary element method, will be discussed in Section III and Appendix A.

In the time domain, the DtN map becomes a time convolution,

[TABLE]

With the DtN map, we incorporate (17) into (10), and we arrive at a reduced problem (2),

[TABLE]

where $\kappa(t)$ is the real-time kernel function which corresponds to $K(s)$ in the Laplace domain.

The main difficulties in implementing the transparent boundary condition (17) are: i) there is no analytical expression for $K(s)$ or $\kappa(t)$ in general, other than the formula (16) in terms of a large-dimensional matrix, in which case it is expensive to evaluate $K(s)$ or $\kappa(t)$ repeatedly. ii) the direct evaluation of the time-convolution integral adds up quickly to the computational cost: Since the kernel function $\kappa(t)$ does not have compact support in time, long time integration is required to evaluate the exact boundary condition.

In light of these concerns, we will introduce further approximations in the next section.

Remark II.1**.**

Our procedure for reducing the problem is reminiscent of reduced-order modeling Bai2002 ; benner2015survey . More specifically, the quantity $\bm{f}_{\Gamma}$ can be viewed as the low-dimensional output, and $\bm{\phi}_{\Gamma}$ corresponds to the control quantity. In the reduced-order literature, it is sometimes convenient to write (11) into the alternative form

[TABLE]

where $s_{0}\in\mathbb{C}$ is a pre-selected scalar. In this case, the kernel function becomes

[TABLE]

If we denote ${\mathcal{A}}=i(H_{\text{I\!I,I\!I}}-is_{0})^{-1}$ and $C=(H_{\text{I\!I,I\!I}}-is_{0}I)^{-1}H_{\text{I\!I,}\Gamma}$, the Taylor expansion of the kernel function $K(s)$ around $s=s_{0}$ reads

[TABLE]

The coefficients (moments) in the expansion are $M_{0}=-H_{\Gamma,\text{I\!I}}C$, $M_{1}=-H_{\Gamma,\text{I\!I}}{\mathcal{A}}C$, $M_{2}=-H_{\Gamma,\text{I\!I}}{\mathcal{A}}^{2}C$, etc. One way to approximate the kernel function $K(s)$ is by rational functions with the same moments, which is known as moment matching Bai2002 ; baker1996pade . One may also notice that (21) is connected to a Krylov subspace $\mathcal{K}({\mathcal{A}},C)$. However, for the problem considered here, the higher powers of ${\mathcal{A}}$ are much more difficult to compute.

III Approximation of the discrete DtN map

Consider the following rational functions:

[TABLE]

$A_{0},\cdots,A_{m-1}$ and $B_{0},\cdots,B_{m-1}$ are $n_{\Gamma}\times n_{\Gamma}$ matrices to be determined. Here the integer $m\geq 0$ will be referred to as the order of the approximation.

The rational approximation reduces the DtN map in the Laplace domain to

[TABLE]

One advantage of the rational approximation is that in the time domain, the dynamics can be represented by an ODE,

[TABLE]

assuming appropriate initial conditions. The superscript (m) denotes the $m$-th derivative. Now, the non-local time-convolution in the DtN map (Eq. (18)) is replaced by a linear ODE system, which is much more efficient in practical implementations. The reduced model (Eq. (2)) is replaced by,

[TABLE]

The remaining question is how to determine the coefficients $A_{i}$ and $B_{i}$. This is done by interpolation. Namely we match $R_{m,m}(s)$ to $K(s)$ at certain points. In principle we need $2m$ points of $(s_{i},K(s_{i}))$ to determine these coefficients. In general, one may hope that the accuracy would improve as the number of interpolation points increases. However, a more subtle issue is the stability of the system (25). The semi-discrete models are quite similar to molecular dynamics models, for which the stability of ABCs has been analyzed in li2009stability . In particular, the Lyapunov functional approach is quite useful. In this paper, we will present the guaranteed stability of the zeroth-order $(m=0)$ and first-order approximation $(m=1)$. We did not find a simple proof of stability for higher order approximations, and we will resort to our numerical simulations to examine the stability property.

Remark III.1**.**

We did not pursue a one-point Padé approximation, or the Krylov subspace projections, which has been overwhelmingly successful in reduced-order modeling Bai2002 . The main reason is that in our case, a selected inversion of $(H_{\text{I\!I,I\!I}}-is_{0})$ can be done very efficiently, but the higher order inverse is very difficult.

III.1 Zeroth order approximation.

For the zeroth order approximation, $R_{0,0}$ becomes a constant matrix, here denoted by $M$. We restrict $M$ to be a Hermitian matrix. The corresponding dynamics is

[TABLE]

Theorem III.2** (Stability condition of the zero order approximation).**

The dynamics (Eq. (26)) is stable if $M$ is chosen as $K(s_{0})$ for an arbitrary positive value of $s_{0}$ ($s_{0}\geq 0$).

Proof.

We define a Lyapunov functional as

[TABLE]

We examine the derivative of the Lyapunov functional,

[TABLE]

This indicates that if $M$ has a negative definite imaginary part, then $\frac{dW}{dt}\leq 0$. To verify this, we note that since $-H_{\text{I,I\!I}}[H_{\text{I\!I,I\!I}}-is_{0}I]^{-1}H_{\text{I\!I,I}}=-H_{\text{I,I\!I}}(H_{\text{I\!I,I\!I}}^{2}+s_{0}^{2}I)^{-1}(H_{\text{I\!I,I\!I}}+is_{0}I)H_{\text{I\!I,I}}$, $K(s_{0})$ always has a negative imaginary part as long as $s_{0}$ is real and positive. Therefore, the dynamics (26) is stable. ∎

This shows that the stability does not depend on the choice of the interpolation point.

Remark III.3**.**

The stability analysis also reveals that the additional term in (26) has an imaginary part that acts as a complex absorbing potential (CAP). However, unlike the commonly used CAP methods riss1996investigation ; Riss_1993 ; MUGA2004357 , this term is derived from the DtN map, and it is dependent of the mesh size, finite difference formulas, as well as the geometry of the domain. Further, the higher order approximations of the DtN map, which will be presented next, provide a more general form that goes beyond the CAP methods.

III.2 First-order approximation.

For the first-order approximation, we use the rational function $R_{1,1}=(sI-B)^{-1}A$ to approximate the DtN map. The first-order approximated dynamics is

[TABLE]

To verify the stability, we introduce a Lyapunov functional,

[TABLE]

where $Q$ is a positive definite matrix to be determined. The derivative of the above Lyapunov functional is

[TABLE]

At this point, we observe that if $iI+QA=0$, and $B^{*}Q+QB$ has negative definite real part, then the dynamics (29) is stable. The following theorem provides a guideline for the choice of the interpolation points to ensure stability.

Theorem III.4** (Stability condition of the first order approximation).**

The dynamics (Eq. (29)) is stable if the coefficients $A$ and $B$ are determined by the interpolations $\lim_{\lambda\to 0}\frac{d}{d\lambda}R_{1,1}^{\prime}=\lim_{\lambda\to 0}\frac{d}{d\lambda}K$, where $\lambda=s^{-1}$, and $R_{1,1}(s_{1})=K(s_{1})$, where $s_{1}$ is a any positive real number.

Proof.

From

[TABLE]

we have

[TABLE]

Now we pick $Q=(H_{\Gamma,\text{I\!I}}H_{\text{I\!I,}\Gamma})^{-1}$. $Q$ is symmetric positive definite. The Lyapunov functional (Eq. (30)) is thus positive definite. Let us denote $K_{1}=K(s_{1})$. When $R_{1,1}(s_{1})=K(s_{1})$, we have,

[TABLE]

Hence, $QB=-is_{1}A^{-1}+iK_{1}^{-1}$. Next, we show that $QB$ has a negative definite imaginary part. We start with

[TABLE]

To continue, we denote $P_{1}=(H_{\text{I\!I,I\!I}}-is_{1}I)^{-1}H_{\text{I\!I,}\Gamma}$. As a result, $K_{1}=-H_{\text{I\!I,}\Gamma}P_{1}$, and

[TABLE]

Therefore we have,

[TABLE]

Since $H_{\text{I\!I,}\Gamma}(H_{\Gamma,\text{I\!I}}H_{\text{I\!I,}\Gamma})^{-1}H_{\Gamma,\text{I\!I}}$ is an orthogonal projection matrix, $H_{\text{I\!I,}\Gamma}(H_{\Gamma,\text{I\!I}}H_{\text{I\!I,}\Gamma})^{-1}H_{\Gamma,\text{I\!I}}-I$ is symmetric negative semi-definite. As a result, $QB+B^{*}Q$ is negative semi-definite. Therefore, the dynamics (Eq. (29)) is stable. ∎

The interpolation scheme suggested in the theorem involves an interpolation point $s_{1}>0$ and another point toward $s_{2}\to+\infty$. Our numerical tests indicate that the approximation is still stable when $s_{2}$ is finite.

III.3 Second-order approximation

The rational function approximation can be extended to higher order. We briefly describe the second-order approximation here, in which case the DtN map is approximated by $R_{2,2}=(s^{2}I-sB_{1}-B_{0})^{-1}(sA_{1}+A_{0})$, and the dynamics augmented with this ABC is given by,

[TABLE]

In the above equation, $\ddot{}$ denotes the second derivative with respect to the time. We have not found a simple proof that provides the stability condition of the dynamics (35) in terms of the interpolation points. The stability will thus be demonstrated by numerical tests in Section IV.

IV Applications & numerical experiments

In this section, we test the ABCs with three examples. For each model, the details regarding the numerical tests are respectively discussed in section IV.1, section IV.2 and section IV.3. The three problems are briefly summarized as follows:

• •

A 1D time-dependent Schrödinger equation extensively used as a test example in the literature antoine2008review , although our main emphasis is on Schrödinger equation in $\mathbb{R}^{3}$. The system is a 1D free electron,

[TABLE]

with the initial condition $\psi^{0}(x)=\exp(-(x-x_{c})^{2}+ik_{0}(x-x_{c}))$.

• •

A 3D time-dependent Schrödinger equation considered in antoine2004 :

[TABLE]

with the initial condition $\psi^{0}(x)=\exp(-x_{1}^{2}-x_{2}^{2}-x_{3}^{2}+ik_{0}x_{1})$

• •

A 3D time-dependent Hartree-Fock model flocard1978three :

[TABLE]

with the initial condition $\varphi_{j}^{0}$ determined from the ground state. The form of $\widehat{H}$ is given by (49) later in this section. The 3D TDHF model is a system of nonlinear 3D time-dependent Schrödinger equations. The Hamiltonian $\hat{H}$ depends on the one-particle wave functions.

Integrators. In general, numerical integrators can be formulated as castro2004propagators

[TABLE]

where $U$ is the operator that mimics the time evolution operator. For linear problems with time-independent potential, the exact operator is a matrix exponential,

[TABLE]

One widely used method is the Crank-Nicholson scheme,

[TABLE]

For the 3D case, it is often impractical to perform the matrix inversion in the Crank-Nicholson scheme. In TDDFT castro2004propagators ; pueyo2018 ; martin2004 , one classical method is the Taylor expansion of the exact integrator,

[TABLE]

Clearly, the operator $U_{5}$ is not unitary. However, we will choose $\Delta t$ to be sufficiently small, in which case this integrator is stable and accurate marques2003 . This allows us to focus more on the performance of various ABCs.

IV.1 The 1D time-dependent Schrödinger equation

In the first test, we look at a 1D quantum system. The setting of the problem is illustrated in FIG. 1.

The analytical solution of Eq. (36) can be explicitly written as

[TABLE]

assuming the initial condition

[TABLE]

The initial condition $\psi^{0}$ is localized around $x_{c}$, which is the center of the wave packet. $k_{0}$ is the wavenumber. In this test, we set $k_{0}=5$ and $x_{c}=-6$. The exact solution, $\psi^{\text{ex}}$, propagates to the right when $k_{0}>0$. Therefore, we only need to implement an ABC on the right boundary. A Dirichlet boundary condition will be imposed on the left.

In our simulations, we pick the interior region to be $\Omega_{\text{I}}=[-12,3]$ and the exterior domain is $\Omega_{\text{I\!I}}=(\infty,-12)\cup(3,\infty)$. The Laplacian operator is discretized by the five-point scheme with grid spacing of $h=0.01$.

The evaluation of the DtN map is discussed in our previous work Wu2018 using the discrete Green’s function. The details of 1D lattice Green’s function will be discussed in the Appendix B. We select $s=20$ for the zeroth-order approximation, two points $s=10$, and $20$ for the first-order approximation, and four interpolation points, $s=10$, $11$, $20$, $21$ for the second-order approximation. The zeroth order approximation corresponds to a complex absorbing potential over two grid points at the boundary. But we point out that in practice, the latter method can be applied to a much larger buffer region. In accordance with the width of the finite-difference stencil, $4$ and $8$ extra variables are introduced in the first-order and second-order approximations, respectively.

The solution computed with the Dirichlet boundary condition is completely reflected back into the interior region when the wave function propagates to the boundary (FIG. 2). The zeroth-order approximation causes some reflection when the wave packet first arrives at the boundary, but most of the reflection is eliminated eventually. The first-order ABC qualitatively captures the transient profile of the exact solution, with some errors when the wave reaches the boundary. The second-order ABC provides a much more accurate solution. The solution with ABCs is systematically improved as the order of the approximation increases.

In Fig. 3, the total electron number ($L_{2}$-norm of the wave function in the computational domain) is presented. Before $t=0.5$, the electron number in four cases is the same. The norm of the wave function does not decay with the fixed boundary condition since all the wave function are reflected. The electron number by the first order ABC decays slower than the exact one. In the second order approximation, The maximum error in the $L_{2}$-norm over time is less than $3\times 10^{-3}$.

In the proposed ABCs, one can either improve the order of ABCs or optimize the interpolation points. In most absorbing techniques, the adjustable parameters are the width of the buffer region or the absorbing strength. It is generally difficult to quantitatively compare the absorbing properties of these methods. As a qualitative comparison, we implemented the following CAP Riss_1993

[TABLE]

in the 1D model. The size of the buffer region is fixed in our simulations. This leads to the effective Hamiltonian $H^{\text{eff}}(\eta)=H-i\eta W$ on the interval $[-16,7]$ where $\eta$ is the CAP strength. The CAP varies from $0.01$ to $10$ (Fig. 3). When $\eta$ is smaller than $0.1$, we still observed the reflection from the buffer region. When $\eta=1$, we observed the decent accuracy ($98\%$). As a comparison, this is also achieved by the proposed first-order approximation without optimizing the interpolation conditions, and without introducing a buffer region.

IV.2 A 3D time-dependent Schrödinger equation

In this section, we will test the absorbing boundary conditions for the 3D Schrödinger equation for a free electron. We restrict the computational domain in a box $[-1.5,1.5]\times[-1.5,1.5]\times[-1.5,1.5]$. The 3D Laplace operator is approximated by the 7-point finite difference scheme in each direction with uniform grid spacing of $h=0.1$. Consequently, there are $31$ interior points and $6$ exterior points in each axis direction. For each interpolation point, the DtN map $K(t)$ is a $14166\times 14166$ dense matrix. The coefficients of the first-order ABC are computed by interpolation of $s=1,2$. Four interpolation points $s=1,2,3,10$ are used for the second-order ABC.

Similar to the 1D case, we still can construct the analytical solution of Eq. (37),

[TABLE]

with $k_{0}=5$. One should notice that the difference between the analytical solution and the exact solution of the discrete model might not be small due to the large grid spacing. Therefore, we will only use the analytical solution as a reference for qualitative comparisons.

For the time integration, the step size is chosen as $\Delta t=0.001$. In the first test, we observe how the total electron density in $\Omega_{\text{I}}$ changes in time (Fig. 4). The number of electrons is almost a constant when we fix the wave function at the boundary (Dirichlet boundary condition). If we impose the ABC on the system, the number of electrons will decrease after the wave function propagates to the boundary. With the first-order ABC, about $20\%$ of the wave function in magnitude is reflected. Much more reflections are reduced by the second-order ABC. Over a longer time period, such reflections occur multiple times for both the first-order ABC and the second-order ABC. After that, almost all the electrons are emitted out of the box.

The top view of surface plots of solutions of 3D time-dependent Schrödinger equation with different boundary conditions are shown in Fig. 5. The fixed boundary condition leads to significant reflections. From $t=0.011$, it does not provide any significant results. These errors are reduced by the absorbing boundary conditions. Some reflections are still observed in the first-order approximation. The reflections in the second-order ABC are almost negligible. This experiment also shows that the proposed ABCs is not sensitive to the presence of corners and edges along the boundary. Even though the wave function propagates to the corners, no significant reflection is observed around the corners.

IV.3 3D time-dependent Hartree-Fock model (TDHF) with localized interactions

The Hartree-Fock equation for a nucleon system can be formulated from the many-body system by approximating the many-body wave function with the Slater determinant, and applying the variational principle of the Skyrme functional flocard1978three ; Skyrme1977 ; maruhn1977time with respect to the wave function. The direct procedure yields the coupled TDHF equations with localized interactions.

[TABLE]

where $H$ is the time-dependent one-body HF Hamiltonian. The one-body Hamiltonian depends on the nucleon density, given by

[TABLE]

The one-body HF Hamiltonian can be explicitly written as

[TABLE]

Here $t_{0}$ and $t_{3}$ are the coefficients of the Skyrme interactions flocard1978three ; Skyrme1977 ; maruhn1977time . Among the five terms in the one-body Hamiltonian, the first term is from the kinetic energy. The following two terms are the expectation value of the zero-range density-dependent two-body effective interaction. Furthermore, $W_{y}$ is the Yukawa potential

[TABLE]

and $W_{C}$ is the Coulomb potential given by

[TABLE]

where $V_{0}$ and $a$ are the coefficients of Yukawa interactions. $\rho_{p}$ is the proton density. In practice, the Yukawa and Coulomb potentials are calculated by solving the corresponding Poisson and Helmholtz problems, respectively,

[TABLE]

We consider the above TDHF model (38) in the infinite region $\mathbb{R}^{3}$. In practice, the computational domain $\Omega_{\text{I}}$ (usually a box) is only a small part of the entire $\Omega$. We make the further assumption that

[TABLE]

Our approach starts from the discrete model. We denote $\bm{\varphi}_{\text{I}}$ is the discretization of the wave function $\bm{\varphi}(\bm{x})$ in $\Omega_{\text{I}}$. $H_{\text{I\!I,I\!I}}$ is the discretization of the operator $-\frac{\hbar^{2}}{2m}\nabla^{2}$. $H_{\text{I,I}}$ is still a nonlinear operator containing the Coulomb potential. The notations of boundaries and DtN map follow those in Section II. In this model, we only assume $\rho=0$ in $\Omega_{\text{I\!I}}$ at $t=0$.

As an application, we study the nuclear reaction of the ${}^{16}\text{O}+^{16}\text{O}$ system in the infinite space $\mathbb{R}^{3}$. The setup and data of the numerical experiments are mainly from flocard1978three . In this model, we assume the perfect spin-isospin degeneracy for each particle, so that each spatial orbital is occupied by four nucleons. There are $32$ nucleons in total in this system. The two particles at ground state are positioned in a box away from each other and away from the boundary. The ground state is achieved by solving the static Hartree-Fock equation self-consistently in $\Omega_{\text{I}}$. We assume that there is no interaction between the two particles at the initial state. The Poisson and Helmholtz problems are solved by preconditioned conjugate gradient method using the same discretization method as the wave functions.

The initial condition is given by multiplying each orbital by the phase $e^{i\bm{k}\cdot\bm{r}}$ to create a head-on collision. Here, $\bm{k}$ should be carefully selected to ensure the particles enter into the fusion window. Namely, two particles pass the Coulomb barrier and are trapped by the nuclear force. Otherwise, the two particles move to the boundaries and the nuclear fusion will not occur. We use the same integrator as the case for three dimensional Schrödinger equation. The time step is set to be $0.001fm/c$. In each time step, we also need to perform the self-consistent iteration due to the nonlinearity pueyo2018 ; castro2004propagators .

In the numerical experiment, $\Omega_{\text{I}}$ is discretized with grid spacing of $1$ fm. The Laplacian operator is approximated by the 7-point scheme in each spatial direction. The region of interest is $[-15fm,15fm]^{3}$. We employ the solution of a larger system ($[-30fm,30fm]^{3}$) as the exact solution to examine the ABCs. The exact number of nucleons and the exact total energy are computed from the larger system restricted to the small region. The initial conditions of the two systems are the nucleon density at the ground state of the smaller system. One should notice that the ground state of the smaller system is not necessarily the ground state of the larger system. The initial condition of the larger system is not the ground state of itself. This will lead a small truncation error due to the long-range interaction in the system.

The energy conservation and mass conservation (FIG. 6) for the standard TDHF with Dirichlet boundary conditions are easily observed in our simulations. The system with ABCs released over $1$ MeV from the total energy, and emitted $0.1$ nucleons in the simulating period. More energy and nucleons are expected to be emitted over a longer period.

Our main focus here is to test how the absorbing property is influenced by the interpolation points. We present results from the following six cases: (A): $s_{1}=10^{-1},s_{2}=2\times 10^{-1}$; (B): $s_{1}=10^{0},s_{2}=2\times 10^{0}$; (C): $s_{1}=10^{1},s_{2}=2\times 10^{1}$; (D): $s_{1}=10^{2},s_{2}=2\times 10^{2}$; (E): $s_{1}=10^{-2},s_{2}=2\times 10^{-2}$; (F): fixed boundary condition.

The case when $s_{1}=1$ and $s_{2}=2$ provides the best agreement with the exact solution in terms of both the number of nucleons and total energy. The case $s_{1}=0.1$ and $s_{2}=0.2$ follows second. When the order of magnitudes of $s_{1}$ and $s_{2}$ are much larger or smaller than 1, the absorbing properties are much worse. The optimal $s_{1}$ and $s_{2}$ should be on the order of 1. From FIG 6, the absorbing property is not sensitive to the selection of $s_{1}$ and $s_{2}$, as long as we choose them in the interval $[0.1,1]$.

Further, we compare the profile of the densities with the exact solution. Results are displayed in FIG 7. In the first 5000 time steps, the error of different cases are almost the same. However, when the wave function reaches the boundaries, the Dirichlet boundary condition generates massive errors as expected. The case (A) and (B) still shows great agreement with the exact density. The case (C), (D), and (E) also improved the accuracy of density in varying degrees.

We would like to remark that the long-range potential (Coulomb and Yukawa) make the problem nonlinear in the entire space, in which case the ABCs are difficult to derive, unless further simplifications are made. However, since the two particles never propagate to the boundary, $\rho\approx 0$ outside. We expect that one can neglect the potential terms in the exterior.

V Summary and discussions

We constructed absorbing boundary conditions for time-dependent Schrödinger equations by first deriving the Dirichlet-to-Neumann map. We chose the starting point to be a semi-discrete approximation so that the resulting boundary condition can be readily implemented with the discretization in the interior seamlessly. The nonlocality in time in the exact boundary condition is treated by rational interpolations of the Laplace transform, which in the time domain turns into linear ODEs. For the zeroth and first-order approximations, the stability has been proved. The effectiveness and accuracy are also illustrated by various numerical tests. For higher-order absorbing boundary conditions, a direct proof seems rather challenging. Our numerical observations are that second-order and third-order methods are still stable for a wide range of interpolation points $s$.

In principle, the boundary conditions presented in this paper can be applied to three-dimensional problems with general geometry. The stability results do not depend on the specific configuration of the computational domain. However, a remaining challenge is to choose the optimal interpolation conditions to maximize the accuracy, i.e., minimize the reflection. For one-dimensional problems, or higher dimensional problems in a half-plane, one can compute the reflection coefficients and choose the optimal boundary condition by minimizing the total reflection li2006variational . This approach breaks down when corners and edges are present. One possible remedy is the optimal interpolation strategy from order-reduction problems anic2008interpolation ; gugercin2008 . For instance, in the current setting, one can formulate the following optimization problem

[TABLE]

for the first-order boundary condition to identify the best interpolation point. This will be pursued in separate works.

Another long-standing issue is associated with the long-range interactions among the electrons, i.e., the Coulomb interactions. This issue has been partially addressed in ehrhardt2006fast , where the Coulomb potential is replaced by an asymptotic form. But in general, all the existing methods have not be constructed to take full account of the Coulomb potential.

Appendix A Evaluation of the Dirichlet-to-Neumann map

In this section, we briefly discuss the evaluation of the kernel (Eq. (16)) in the Laplace domain. The key is to use the Green’s function to evaluate the selected inversion. We refer readers to our previous work Wu2018 ; Li2012 for the detailed derivations.

The discrete Green’s function $\widetilde{G}_{ij}$ corresponding to the discretized operator $\widetilde{H}_{ij}$ satisfies

[TABLE]

To shorten the derivation, the notations of $\widetilde{G}_{\text{I,I\!I}}$ and $\widetilde{G}_{\text{I\!I,I\!I}}$ follow the convention of $H_{\text{I,I\!I}}$ and $H_{\text{I\!I,I\!I}}$ respectively. Furthermore, we introduce the boundary set in $\Omega_{\text{I\!I}}$,

[TABLE]

and matrices

[TABLE]

We start with the discretized operator $\tilde{H}$ in $\Omega$. By (11), the discretized operator satisfies

[TABLE]

The definition of the discrete Green’s function implies that

[TABLE]

Eq. (56) and Eq. (57) can now be combined into the following equation

[TABLE]

where the fact that $\widetilde{H}_{\text{I,I\!I}}=H_{\text{I,I\!I}}$ has been used. By multiplying both sides of (58) by $H_{\Gamma,\text{I\!I}}$, we have

[TABLE]

Notice that $\bm{f}_{\Gamma}=H_{\Gamma,\text{I\!I}}\bm{\Phi}_{\text{I\!I}}$. Making use of the sparsity of $H_{\Gamma,\text{I\!I}}$ and $H_{\text{I,I\!I}}$, we obtain

[TABLE]

Therefore, the discrete Dirichlet-to-Neumann map is expressed as

[TABLE]

This shows that once the Green’s function is available, the evaluation of the DtN map only involves the operations of small matrices.

Appendix B Discrete Green’s function for 1D Schrödinger equation

Consider the one-dimensional time-dependent Schrödinger equation,

[TABLE]

By taking the Laplace transform in time ($t\rightarrow s$) and Fourier transform in space ($x\rightarrow q$), the corresponding fundamental solution solves

[TABLE]

which gives,

[TABLE]

The inverse Fourier transform yields the fundamental solution in the Laplace domain

[TABLE]

$G(s,x)$ is the continuous Green’s function of the following equation in the Laplace domain. Now, we turn to the discrete case. To obtain the discrete Green’s function, we discretize the Schrödinger equation in the Laplace domain by the five-point scheme,

[TABLE]

where $c_{-2}=-1/12$, $c_{-1}=4/3$, $c_{0}=-5/2$, $c_{1}=4/3$, and $c_{2}=-1/12$. $h$ is the grid spacing.

The characteristic equation of equation (66) is given by

[TABLE]

The four roots of the characteristic equation are given by

[TABLE]

where $a=\sqrt{9+6sh^{2}i}$, $Re(u_{1})<1$ and $Re(u_{3})<1$ for a small $s$. We take the branch cut ($u_{1}$ and $u_{3}$) to ensure the boundedness of the Green’s function. Hence, the discrete Green’s function can be explicitly expressed as

[TABLE]

where

[TABLE]

The coefficients $b_{1}$ and $b_{2}$ are determined by the conditions

[TABLE]

We make a comparison between the continuum Green’s function and analytical discrete Green’s function in FIG. 8. Both the imaginary part and real part of the two Green’s functions get close to each other when the distance gets large. This reveals that we can use the continuous Green’s function to approximate the discrete Green’s function when the distance is large, which becomes quite important especially when there is no explicit formula for the discrete Green’s function.

Appendix C First order ABC for 1D Schrödinger equation

Let us consider the 1D Schrödinger equation (Eq. (62)) on the interval $\Omega_{\text{I}}=[x_{l},x_{r}]$. We assume the discretization (66) of eq. (62). $x_{1},\cdots,x_{N}$ are $N$ grid points inside the $\Omega_{\text{I}}$. $x_{-1}$, $x_{0}$, $x_{N+1}$, and $x_{N+2}$ are the exterior points. The input and output of the DtN map, respectively, are

[TABLE]

The matrices in the DtN map are explictly written as

[TABLE]

and

[TABLE]

Therefore, the kernel function $K(s)$, which is defined by (15), can be evaluated (61). Even for the 1D model, the DtN map is not trivial. To demonstrate the idea of the proposed ABCs, we determine the coefficients of the first order ABC is by the moments of $K(s)$. More specifically, the coefficients of the first order ABC are obtained by solving

[TABLE]

for certain $s_{1}>0$. Generally, the coefficients $A$ and $B$ can not be explicitly expressed.

Appendix D Discrete Green’s function for the 3D Schrödinger equation

Consider the 3D time-dependent Schrödinger equation,

[TABLE]

By following the same procedure, the corresponding fundamental solution in the Laplace space is given by

[TABLE]

where $\bm{x}=(x_{1},x_{2},x_{3})$ and $r=\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$.

For the finite-difference approximation, we discretize the operator, e.g., by the 9-point scheme in each spatial direction.

[TABLE]

where $c_{-4}=-1/560$, $c_{-3}=8/315$, $c_{-2}=-1/5$, $c_{-1}=8/5$, $c_{0}=205/72$, $c_{1}=8/5$, $c_{2}=-1/5$, $c_{3}=8/315$, and $c_{4}=-1/560$. We denote $c_{\bm{k}}$ as the coefficients of the 3D discretization.

The discrete Green’s function can be expressed as a Fourier integral,

[TABLE]

Here $B$ refers to the Fourier domain associated with the finite-difference grid points. It is given by: $\frac{\pi}{h}[-1,1]\times\frac{\pi}{h}[-1,1]\times\frac{\pi}{h}[-1,1]$ for uniform grids.

In FIG 9, we show an example of the discrete Green’s function, obtained by a quadrature with $100\times 100\times 100$ points in the Fourier domain, compared with the continuous Green’s function.