# Absorbing boundary conditions for the time-dependent Schrödinger-type equations in R^3

Absorbing boundary conditions are presented for three-dimensional time-dependent Schrödinger-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ödinger 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.

## Authors

• 2 publications
• 2 publications
05/28/2020

### Time-dependent acoustic scattering from generalized impedance boundary conditions via boundary elements and convolution quadrature

Generalized impedance boundary conditions are effective, approximate bou...
11/04/2021

### Local Compatibility Boundary Conditions for High-Order Accurate Finite-Difference Approximations of PDEs

We describe a new approach to derive numerical approximations of boundar...
12/28/2021

### Reduced order modeling with time-dependent bases for PDEs with stochastic boundary conditions

Low-rank approximation using time-dependent bases (TDBs) has proven effe...
06/18/2020

### Bérenger's PML on Rectangular Domains for Pauli's Equations

This article proves the well posedness of the boundary value problem tha...
09/02/2020

### A numerical study of third-order equation with time-dependent coefficients: KdVB equation

In this article we present a numerical analysis for a third-order differ...
04/06/2020

### Stable Boundary Conditions and Discretization for PN Equations

A solution to the linear Boltzmann equation satisfies an energy bound, w...
04/11/2022

### Linear Moment Models to Approximate Knudsen Layers

We propose a well-posed Maxwell-type boundary condition for the linear m...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 three approaches: The exterior complex scaling (ECS)

Scrinzi2014; NICOLAIDES197811; SIMON1979211; mccurdy2004solving, where the coordinate outside a fixed radius is scaled into the complex plane, the mask function method marques2003; tafipolsky2006; Wang201107, where the wave functions are gradually scaled to zero, and complex absorbing potentials (CAP) riss1996investigation; Riss_1993; MUGA2004357, which is introduced outside the boundary by adding a complex potential to the Hamiltonian. 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). 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 that 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.

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

This paper is conerned with the time-dependent Schrödinger equations in a domain ,

 i∂∂tϕj(x,t)=^Hϕj(x,t),x∈Ωj=1,…N. (1)

Here we have scaled the Planck constant and mass to unity. is the imaginary unit. is the one-particle wave function for an electron or nucleon in the -dimensional quantum system. For convenience, let us denote

to be a vector-valued function that collects all the wave functions. In general, the Hamiltonian operator can be written as

, with being the total potential of the system. Eq. (1), when posed in three dimensions, 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.

We suppose that is a subdomain of interest in . Rather than solving the Schrödinger equations in the entire domain we take into account the influence from the surrounding region, here denoted by , by deriving ABCs at the boundary of . . In the exterior region , we assume that the operator does not explicitly depend on . More specifically, , or more generally a constant, for any .

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 . 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 be the set of grid points with a constant spacing of in the interior and grid points in the exterior . . can be infinite. In real-space methods octopus, the Hamiltonian operator in (1) is usually approximated by finite-difference formulas. 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 , can be approximated by,

 ∂2∂x2ψ(x,t)≈−ψ−2+16ψ−1−30ψ0+16ψ1−ψ212h2.

Here with 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 , and . 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:

 (2)

where and . and . Following the same partition, the discretization of the Hamiltonian operator H is structured as follows,

 H=[HI,IHI,I\!IHI\!I,IHI\!I,I\!I]. (3)

Since the Schrödinger equation (1) has been descretize in space, we use to denote the time derivatives hereafter.

Our goal is to simplify the second equation in (2), while retaining the first equation. and are respectively the discretized Hamiltonian operator in and . can be a nonlinear operator. and are the off-diagonal blocks representing the coupling between and . 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 . To this end, we denote as the boundary of , such that

 Γ={xj∈ΩI| if there exists xk∈ΩI\!I such that Hjk≠0}. (4)

The width of depends on the width of the finite-difference stencil. We define as a vector formed by all for , where is the number of grid points in . The vector can be reordered so that the first components are associated with the grid points in , i.e.,

 ϕI=[ϕΓϕI∖Γ]. (5)

The boundary region is defined in such a way that there is no direct coupling between the points in the interior of and those in . This is reflected in the off-diagonal block of the Hamiltonian,

 HI\!I,I=[HI\!I,ΓHI\!I,I∖Γ]=[HI\!I,Γ0]. (6)

In general, this procedure can be carried out by defining a restriction operator to extract the components of a function that correspond to grid points at the boundary from a function defined in . Namely,

 ϕΓ=EϕI. (7)

With the reordering in (5), the matrix can be explicitly expressed as . Further, (6) can now be simply written as . Since , we also have,

 HI,I\!I=ETHΓ,I\!I. (8)

We are now set to formulate the ABC. We first consider the influence of the wave functions in on the wave functions in , and we define

 fΓ=HΓ,I\!IϕI\!I. (9)

 ˙ϕI(t)=−iHI,IϕI(t)−iETfΓ(t). (10)

At this point, we observe that what is needed to compute in time is For this purpose, let us take the Laplace transform of the second equation in (2). Due to the assumption that one has,

 isΦI\!I(s)=HI\!I,I\!IΦI\!I(s)+HI\!I,IΦI(s), (11)

where and . By defining , equation (11) simply reads

 ˜HI\!I,I\!I(s)ΦI\!I(s)=−HI\!I,I% ΦI(s). (12)

Formally, its solution can be expressed as

 ΦI\!I(s)=−˜H−1I\!I,I\!IHI\!I% ,IΦI(s). (13)

Due to the locality of , the right hand side of (13) can be further reduced. Notice that,

 HI\!I,IΦI=HI\!I,ΓEΦ% I=HI\!I,ΓΦΓ.

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

 ΦI\!I(s)=−˜H−1I\!I,I\!IHI\!I% ,ΓΦΓ(s). (14)

Since is a matrix, it is in general impractical to solve the linear system (14) numerically. However, the main observation from (10) is that we only need to keep the computation in . We let be the Laplace transform of . By left-multiplying the equation above by , we obtain,

 FΓ(s)=K(s)ΦΓ(s). (15)

The matrix-valued function , which will play a key role in the numerical approximation, is given by,

 K(s)=−HΓ,I\!I˜H−1I\!I,I\!I(s)HI% \!I,Γ=−HΓ,I\!I[HI\!I,I\!I−isI]−1HI\!I,Γ. (16)

The mapping is precisely the Dirichlet-to-Neumann (DtN) map of the problem (2) in the Laplace domain. Although 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,

 fΓ(t)=∫t0κ(t−τ)ϕΓ(τ)dτ. (17)

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

 ˙ϕI(t)=−iHI,IϕI(t)−iET∫t0κ(t−τ)EϕI(τ)dτ, (18)

where is the real-time kernel function which corresponds to in the Laplace domain.

The main difficulties in implementing the transparent boundary condition (17) are: i) there is no analytical expression for or in general, other than the formula (16) in terms of a large-dimensional matrix, in which case it is expensive to evaluate or repeatedly. ii) the direction evaluation of the time-convolution integral adds up quickly to the computational cost: Since the kernel function 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 can be viewed as the low-dimensional output, and corresponds to the control quantity. In the reduced-order literature, it is sometimes convenient to write (11) into the alternative form

 i(s−s0)(HI\!I,I\!I−is0I)−1ΦI\!I(s)=ΦI\!I(s)+(HI\!I,I\!I−is0I)−1HI\!I,IΦI(s), (19)

where is a pre-selected scalar. In this case, the kernel function becomes

 K(s)=−HΓ,I\!I(I−i(s−s0)(HI\!I,I\!I−is0I)−1)−1(HI\!I,I\!I−is0I)−1HI\!I,Γ. (20)

If we denote and , the Taylor expansion of the kernel function around reads

 K(s)=−HΓ,I\!I(C+AC(s−s0)+A2C(s−s0)2+A3C(s−s0)3+⋯). (21)

The coefficients (moments) in the expansion are

, , , etc. One way to approximate the kernel function 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 . However, for the problem considered here, the higher powers of are much more difficult to compute.

## Iii Approximation of the discrete DtN map

Consider the following rational functions:

 Rm,m(s)=(sm−sm−1B0−⋯−Bm−1)−1(sm−1A0+⋯+Am−1). (22)

and are matrices to be determined. Here the integer will be referred to as the order of the approximation.

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

 FΓ(s)=Rm,m(s)ΦΓ(s). (23)

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

 f(m)Γ=B0f(m−1)Γ+⋯+Bm−1fΓ+A0ϕ(m−1)Γ+⋯+Am−1ϕΓ, (24)

assuming appropriate initial conditions. 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,

 ⎧⎪⎨⎪⎩∂∂tϕI(t)=−iHI,IϕI(t)−iHI,I\!IϕI\!I(t),f(m)Γ=B0f(m−1)Γ+⋯+Bm−1fΓ+A0ϕ(m−1)Γ+⋯+Am−1ϕΓ. (25)

The remaining question is how to determine the coefficients and . This is done by interpolation. Namely we match to at certain points. In principle we need points of 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 and first-order approximation . 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 can be done very efficiently, but the higher order inverse is very difficult.

### iii.1 Zeroth order approximation.

For the zeroth order approximation, becomes a constant matrix, here denoted by . We restrict to be a Hermitian matrix. The corresponding dynamics is

 ˙ϕI(t)=−iHI,IϕI(t)−iETMEϕI(t). (26)
###### Theorem III.2 (Stability condition of the zero order approximation).

The dynamics (Eq. (26)) is stable if is chosen as for an arbitrary positive value of ().

###### Proof.

We define a Lyapunov functional as

 W(t)=ϕ∗I(t)ϕI(t). (27)

We examine the derivative of the Lyapunov functional,

 dWdt =ddtϕ∗IϕI+ϕ∗IddtϕI (28) =i(ϕ∗IHI,I+ϕ∗ΓM∗E)ϕI−iϕ∗I(HI,IϕI+ETMϕΓ) =iϕ∗Γ(M∗−M)ϕΓ.

This indicates that if has a negative definite imaginary part, then . To verify this, we note that since , always has a negative imaginary part as long as 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 to approximate the DtN map. The first-order approximated dynamics is

 ⎧⎪ ⎪⎨⎪ ⎪⎩ddtϕI(t)=−iHI,IϕI(t)−iETfΓ(t),ddtfΓ(t)=BfΓ(t)+AEϕI(t). (29)

To verify the stability, we introduce a Lyapunov functional,

 W(t)=ϕ∗I(t)ϕI(t)+f∗Γ(t)QfΓ(t), (30)

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

 dWdt= ddtϕ∗IϕI+ϕ∗IddtϕI+ddtf∗ΓQfΓ+f∗ΓQddtfΓ (31) = (iϕ∗IHI,I+if∗ΓE)ϕI+ϕ∗I(−iHI,Iϕ%I−iETfΓ) +(f∗ΓB∗+ϕ∗IETA∗)QfΓ+f∗ΓQ(BfΓ+AEϕI) = f∗Γ(iI+QA)EϕI+ϕ∗% IET(−iI+A∗Q)fΓ+f∗Γ(B∗Q+QB)fΓ.

At this point, we observe that if , and 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 and are determined by the interpolations , where , and , where is a any positive real number.

###### Proof.

From

 limλ→0ddλR1,1=limλ→0ddλK,

we have

 A=−iHΓ,I\!IHI\!I,Γ.

Now we pick . is symmetric positive definite. The Lyapunov functional (Eq. (30)) is thus positive definite. Let us denote . When , we have,

 B=s1I−AK−11.

Hence, . Next, we show that has a negative definite imaginary part. We start with

 QB+B∗Q =2s1(HΓ,I\!IHI\!I,Γ)−1+i(K−11−(K∗1)−1), (32) =2s1(HΓ,I\!IHI\!I,Γ)−1+i(K∗1)−1(K∗1−K1)K−11.

To continue, we denote . As a result, , and

 K∗1−K1 =−HΓ,I\!I(HI\!I,I\!I+is1I)−1HI\!I,Γ+HΓ,I\!I(HI\!I,I\!I−is1I)−1HI\!I,Γ, (33) =2is1P∗1P1.

Therefore we have,

 QB+B∗Q =2s1(K∗1)−1(K∗1(HΓ,I\!IH% I\!I,Γ)−1K1−P∗1P1)K−11 (34) =2s1(K∗1)−1(P∗1HΓ,I\!I(HΓ,I\!IHI\!I,Γ)−1HI\!I,ΓP1−P∗1P1)K−11, =2s1(K∗1)−1P∗1(HI\!I,Γ(HΓ,I\!IHI\!I,Γ)−1HΓ,I\!I−I)P1K−11.

Since is an orthogonal projection matrix, is symmetric negative semi-definite. As a result, is negative semi-definite. Therefore, the dynamics (Eq. (29)) is stable. ∎

The interpolation scheme suggested in the theorem involves an interpolation point and another point toward . Our numerical tests indicate that the approximation is still stable when 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 , and the dynamics augmented with this ABC is given by,

 {˙ϕI(t)=−iHI,IϕI(t)−iETfΓ(t),¨fΓ(t)=B1˙fΓ(t)+B0fΓ(t)+A1E˙ϕI(t)+A0EϕI(t). (35)

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 . The system is a 1D free electron,

 i∂∂tψ(x,t)=ˆHψ(x,t),ˆH=−∂2x2 in R, (36)

with the initial condition .

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

 i∂∂tψ(x,t)=ˆHψ(x,t),ˆH=−∇22 in R3, (37)

with the initial condition

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

 i∂∂tφj(x,t)=ˆHφj(x,t) in R3, for j=1,…,A (38)

with the initial condition determined from the ground state. The form of is given by (47) later in this section. The 3D TDHF model is a system of nonlinear 3D time-dependent Schrödinger equations. The Hamiltonian depends on the one-particle wave functions.

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

 ϕ(n+1)=Uϕ(n), (39)

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

 UE(t,t′)=exp(−iΔtH(t′)). (40)

One widely used method is the Crank-Nicholson scheme,

 UCN=(1+iΔt/2H)−1(1−iΔt/2H). (41)

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,

 U5=I−iHΔt−12H2(Δt)2+16H3(Δt)3−i124H4(Δt)4. (42)

Clearly, the operator is not unitary. However, we will choose 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

 ψex(x,t)=√ii−2texp(−k0(x−xc)+k20t−i(x−xc)2i−2t), (43)

assuming the initial condition

 ψ0(x)=exp(ik0(x−xc)−(x−xc)2). (44)

The initial condition is localized around , which is the center of the wave packet. is the wavenumber. In this test, we set and . The exact solution, , propagates to the right when . 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 and the exterior domain is . The Laplacian operator is discretized by the five-point scheme with grid spacing of .

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 four interpolation points, , , , for the second-order approximation and two points , and for the first-order approximation.

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 first-order ABC qualitatively captures the profile of the exact solution, with some errors when the wave reaches the boundary. The second-order ABC provides a much more accurate solution.

In Fig. 3, the total electron number (-norm of the wave function in the computational domain) is presented. Before , 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 -norm over time is less than .

### 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 . The 3D Laplace operator is approximated by the 7-point finite difference scheme in each direction with uniform grid spacing of . Consequently, there are interior points and exterior points in each axis direction. For each interpolation point, the DtN map is a dense matrix. The coefficients of the first-order ABC are computed by interpolation of . Four interpolation points are used for the second-order ABC.

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

 ψex=(ii−2t)32exp⎛⎝−i(x21+x22+x23)−k0x1+12k20ti−2t⎞⎠ (45)

with . 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 . In the first test, we observe how the total electron density in 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 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 , 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)

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

 iℏ∂∂tφj(t)=H(t)φj(t),j=1,…,A, (46)

where is the time-dependent one-body HF Hamiltonian. This one-body HF Hamiltonian can be explicitly written as

 H=−ℏ22mΔ+34t0ρ+316t3ρ2+Wy+WC. (47)

Here and are the coefficients of the Skyrme interactions flocard1978three; Skyrme1977; maruhn1977time. The one-body Hamiltonian depends on the nucleon density, given by

 ρ(r,t)=A∑j=1|φj(r,t)|2. (48)

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, is the Yukawa potential

 Wy(r)=V0∫drexp−|r−r′|/a|r−r′|/aρ(r′), (49)

and is the Coulomb potential given by

 WC(r)=e2∫dr′1|r−r′|ρp(r′), (50)

where and are the coefficients of Yukawa interactions. is the proton density. In practice, the Yukawa and Coulomb potentials are calculated by solving the corresponding Poisson and Helmholtz problems, respectively,

 (51)

We consider the above TDHF model (38) in the infinite region . In practice, the computational domain (usually a box) is only a small part of the entire . We make the further assumption that

 ψj=0 for i=1,…,A and ρ=0 in ΩI% \!I, at t=0. (52)

Our approach starts from the discrete model. We denote is the discretization of the wave function in . is the discretization of the operator . 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 in at .

As an application, we study the nuclear reaction of the system in the infinite space . 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 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 . 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 to create a head-on collision. Here, 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 . In each time step, we also need to perform the self-consistent iteration due to the nonlinearity pueyo2018; castro2004propagators.

In the numerical experiment, is discretized with grid spacing of fm. The Laplacian operator is approximated by the 7-point scheme in each spatial direction. The region of interest is . We employ the solution of a larger system () 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 MeV from the total energy, and emitted 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): ; (B): ; (C): ; (D): ; (E): ; (F): fixed boundary condition.

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

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, 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 .

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

 s0=argmin∥K(s)−R1,1(s)∥H2, (53)

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 corresponding to the discretized operator satisfies

 ∑j˜Hi,j˜Gjk=δi,k,for any % integer i,k.

To shorten the derivation, the notations of and follow the convention of and respectively. Furthermore, we introduce the boundary set in ,

 Σ={xj∈ΩI\!I| if there exists xk∈ΩI such that Hjk≠0}. (54)

and matrices

 ˜GΓ,Σ=[˜Gij]i∈Γ,j∈Σ,˜GΣ,Σ=[˜Gij]i∈Σ,j∈Σ,HΓ,Σ=[Hij]i∈Γ,j∈Σ,HΣ,Γ=[Hij]i∈Σ,j∈Γ.

We start with the discretized operator in . By (11), the discretized operator satisfies

 ˜HI\!I,I\!IΦI\!I+˜HI% \!I,IΦI=0. (55)

The definition of the discrete Green’s function implies that

 ˜GI\!I,I\!I˜HI\!I,I\!IΦI\!I+˜GI\!I,I˜HI,I\!IΦI\!I=ΦI\!I. (56)

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

 ΦI\!I=˜GI\!I,IHI,I\!IΦI\!I−˜GI\!I,I\!IHI\!I,IΦI, (57)

where the fact that has been used. By multiplying both sides of (57) by , we have

 HΓ,I\!IΦI\!I=HΓ,I\!I˜GI\!I,I\!IHI,I\!IΦI\!I−HΓ,%II˜