Numerically Absorbing Boundary Conditions for Quantum Evolution Equations

Transparent boundary conditions for the transient Schrodinger equation on a domain Ω can be 
derived explicitly under the assumption that the given potential V is constant outside of this 
domain. In 1D these boundary conditions are non-local in time (of memory type). For the 
Crank-Nicolson finite difference scheme, discrete transparent boundary conditions are 
derived, and the resulting scheme is proved to be unconditionally stable. A numerical example 
illustrates the superiority of discrete transparent boundary conditions over existing ad-hoc 
discretizations of the differential transparent boundary conditions. 
As an application of these boundary conditions to the modeling of quantum devices, a transient 
1D scattering model for mixed quantum states is presented.


INTRODUCTION
The formulation and implementation of physically reasonable and mathematically well-posed boundary conditions (BC) is one of the big open problems and challenges for transient simulations of semiconductor devices through quantum mechanical models.This paper is concerned with the construction and discreti- zation of absorbing boundary conditions (ABC) for the Schr6dinger equation (SE) ihtt 2 At + V(x,t)t, x E v, t > O, v(x,0) (1.1) where the electrostatic potential V is assumed to be given.These BC's are derived from mathematical * E-mail: amold@math.tu-beflin.de313 considerations, and they can be used as ingredients for modeling quantum contacts.
When numerically solving a whole-space evolution problem, the computation has to be confined to a finite domain by introducing artificial BC's.If the ini- tial data is supported on this finite domain D,, one can approximate lfa, the exact solution of the whole space problem restricted to D,, by solving the original problem only on D,, together with ABC's on /)if2.If this approximate solution coincides on f with the exact solution, one refers to these BC's as transparent boundary conditions (TBC).
ABC's were first derived by Engquist and Majda for hyperbolic systems by requiring that outgoing waves can leave the computational domain without being reflected back in ( [7]).Since then, these ideas have been adapted and refined for numerous applica- tions.E.g., ABC's for the Wigner equation of quan- tum mechanics were obtained and analyzed in [16], [1 ], and they were used for transient simulations of resonant tunneling devices in 10].
The recent interest in finding reasonable BC's for the SE was mostly prompted by applications from outside of quantum mechanics.We remark that a Schr6dinger-type equation naturally appears as a small-angle approximation to the Helrnholtz equation in cylindrical coordinates.There, the radial variable r plays the role of the time in (1.1), and the axial varia- ble z the role of the spatial variable x.In optical appli- cations this model is referred to as the Fresnel equation ( [17]), and in underwater acoustics it is called the parabolic equation ([ 15]).
In previous simulations, several heuristic BC's were introduced for the SE, and they yield reasonable results for either short time calculations or a limited frequency range.These existing strategies include the introduc- tion of an artificial absorption or attenuation layer close to the boundary, which can also be interpreted as add- ing a complex potential in (1.1) (see 11 ], 19], [20]).In another approach, the wave function is fitted to a plane wave close to the boundary ( [8], [20]).
Our goal is to first derive (homogeneous) TBC's for the 2D SE under the two assumptions, that the ini- tial data xlr is supported in D,, and that the medium (here modeled by the potential V) is homogeneous and time-independent outside of .O therwise infor- mation would be lost, and the whole space evolution could not be mimicked on the finite domain For the 1D SE, the TBC was independently derived by several authors from various application fields ([ 141, [41, [9]).
As this general construction may be useful for many other applications, we will sketch this deriva- tion in 2 for the 1D and 2D Schr6dinger equation, and analyze the well-posedness of the resulting initial boundary value problems (IBVP).In 3 we derive and investigate an unconditionally stable numerical discretization of the 1D problem, based on a Crank-Nicolson finite difference scheme.As an appli- cation of these TBC's we present in 4 a transient Schr6dinger-Poisson scattering model for mixed quantum states.In this Section we will first derive the TBC's for the 2D Schr6dinger equation on the slab f2 {(x,y) 210 < x < L }.We assume that the initial data xr is supported in g2, and that the given potential is constant outside: V(x,y,t) 0 for x < O, V(x,y,t) V/ for x > L.
We first cut the original whole space problem (1.1) into three subproblems, the interior problem on D, and a left and right exterior problem.They are cou- pled by the assumption that V, Vx are continuous across the boundary.The interior problem reads h2 ihwt 2 AW + V(x,y,t)W, (x,y) E g2, > O, V(x,y, 0) xlt (x,y), (2.1) tx(O,y,t) (T_v)(O,y,t), XVx(L,y,t (T+x)(L,y,t).
Here, T__.denote the Dirichlet-to-Neumann maps at the boundaries, and they are obtained by solving the two exterior problems: O(y,t), yE IR, > 0, (T_O)(y,t) vx(O,y,t), and analogously for T/.In (2.2) we also have to require that v decays as x -oo.Since the potential is constant in the exterior problems, we can solve them explicitly by the Fourier-Laplace method and thus obtain the two boundary operators T_ needed in (2.1).
The Fourier-Laplace transform of v is O(x, k, s) v(x,y, t)e -iky e-"' dy dt, ( where we set s rl + i, , and rl > 0 is fixed, with the idea to later perform the limit rl --0.Now Since its solutions have to decrease as x ----oo, we obtain eX/k2-2is/ hxt(k, s) (2.) Hence the transformed Dirichlet-to-Neumann opera- tor T_ reads T_O(k,s) e-+ k2(k,s), h k e + ig >o, (. and T+ is calculated analogously.Now we will first discuss e 1D situation.For- mally setting k 0 in (2.6), e left TBC is obtained by inverse Laplace transfo: ,(0,tl-e-N Similly, the right TBC is derived as 2e-ie-,V+t/d Zt v(L,) dS dx..8)ese BC's e non-local in and of memory-type, thus requiring e storage of all the past history at e boundy in a numerical discretization.
We remk that e term can be interpreted as a fractional , time derivative.A simple calculation also shows that (2.7) is equivalent to the impedance boundary condition derived in [14]" nifootVX(O'Z)dz. (2.10)   v(O,t) V .eZx/t 'x From the used construction, it is clear that the 1D Schr6dinger equation with the TBC's (2.7) and (2.8)   has a solution.For regular enough initial data, e.g.x H6 (0, L), the whole space solution (x,t) will sat- isfy the TBC's at least in a weak sense.The unique- ness of the solution is, however, not trivial.In order to prove uniform boundedness of II(.,t)ll0, ) in we will need the following simple lemma.
A straight forward calculation with the ID Schr6dinger equation now shows IIv(t)ll./0,/_<IIvll/0,,./,> 0, (2.13) and this implies uniqueness of the solution to the Schr6dinger IBVP.In this calculation one needs the estimate (2.11) for the left boundary term, and an analogous estimate for the right boundary. (2.13)  reflects the fact that some of the initial mass or particle density n(x, t) --I(x,t)l 2 leaves the computational domain [0, L] during the evolution.In the whole space problem IIv(/)ll )is of course conserved.
We now return to the 2D problem on the vertical strip 0 < x< L. For this case we derived the TBC in (2.6).After inverse Fourier and Laplace transforma- tion, this TBC is a convolution in and y, generalizing the explicit 1D representation (2.7).From a numerical point of view, an implementation of this non-local (in and y) BC would be rather costly.It is therefore desirable to first approximate this TBC by a BC that is at least local in y.
In order to motivate our approximation, we first sketch an alternative, heuristic derivation of the 2D TBC (2.6) at x --0.We consider left traveling plane wave solutions to the free SE (since the potential V--0 close to the left boundary) (x,y,t) -exp[i(-v/2----k2x+ky-tot)] (2.14   with to > 0 and k fixed.At x 0 this plane wave satis- fies oV/2t-O [/x -l k21t, (2.15) and this coincides with (2.6) if we set co is.One eas- ily sees that 0, the impact angle of the plane wave at the boundary (measured towards the normal) is given by sin 2 0 hk2/(2c0).If outgoing waves hit the bound- ary almost orthogonally (0 0), it is a standard strat- egy (see [7]) to use a Taylor approximation in 0 of the symbol of the pseudo-differential operator T_ (see (2.6)): s (ih) +ik 2= {/qcos0 or 9q 1+ (2. 16) In lowest order this yields a first order ABC for the 2D SE at x 0: i9 for v(O,y, r,) dx.

Vq
(2.17) This BC essentially neglects 2D effects at the bound- ary, and the well-posedness of the resulting IBVP car- ries over from the 1D situation.On finite time intervals, this solution can be expected to be a reason- able approximation to the whole space solution, under the assumption 0 0.
In the next order approximation from (2.16) we first get after multiplying by + +v-sfltx(O,k,s) e -i s + @(0,k,s).(2.18) And by inverse Laplace-transformation, the second order ABC then reads e'-i(vt--Vyy), x-0. (2. 19)  We conjecture that the resulting IBVP is well-posed; this problem will be analyzed in a forthcoming paper.

DISCRETIZATION OF THE TRANSPARENT BOUNDARY CONDITION IN 1D
In this Section we will discuss a discretization of the 1D TBC (2.7) based on the Crank-Nicolson finite dif- ference scheme for the SE.With the uniform grid points xj--jkx, n nat, and the approximations ' -(x.,t,) this scheme reads for the whole space problem r(j +1 j) j++l 2j +' + j-+l + j+, -2xj+xj_, +wV;+1/2(j +' +jj), (3.1)   with For our analysis, one of the main advantages of this second order (in zXx and At) scheme is, that it is unconditionally stable, and it preserves the discrete L2-norm IIv"I[2 2 Ax E IV'I )- jeT/ We remark that most existing discretization schemes for the 1D Schr6dinger equation with TBC's are also based on Crank-Nicolson finite differences ([4], [13], [15]).In [17] a semi-discrete TBC for the 1D SE with time-dependent potentials in the exterior domains (i.e., V/ V/(t)) was derived.For constant potentials, the obtained method (also Crank-Nicol- son) is related to our ideas.
The delicate question is how to discretize the con- volution (2.7) with its singular kernel.In [13], Mayfield used the approximation ftu Us(O, tN --"r.) N-, ft,,+, dr, ,a e E ,g, " ,go .for the C in e equivalent foulation (2.10).For e resulting scheme she obtained e following result: Our aim is to derive a different boundaw discreti- zation of (2.7) which gives an unconditionally stable scheme and a discrete analogue of (2.13), the unifo boundedness of the continuous L2-norm.When con- sidering the discretization of TBC's, it should be a standard strategy to derive the discrete TBC's of the fully discretized problem, rather then attempting to discretize the differential TBC.Discrete TBC's not only completely avoid any numerical reflections at the boundary, but their numerical stability is often also better-behaved, as it is the case here.A compari- son of these two strategies for a ID wave propagation problem is given in [6].
To derive the discrete TBC we will now mimic the derivation from 2 on a discrete level, and we will again only consider the lbft BC.The idea is to explic- itly solve the discrete exterior problems.The Z-transform of the sequence {'}, n 0, with j considered fixed, is defined as the Laurent series Z{xj} j(z) jz -n.
The Pn only have to be evaluated at the two values go and gj, and hence the numerically stable recursion formula for the Legendre polynomials can be used.
The system matrix of the implicit scheme (3.1), (3.8) is diagonal-dominant and hence invertible.To prove stability of the scheme, one can derive a dis- crete analogue of Lemma 2.1 and (2.13).We then have the main result of this Section: Theorem 3.3 The solution of the discretized SE (3.1) with the discrete TBC's (3.8) is uniformly bounded J-I IIv"ll 2< IIvll , n _> l, (3.12)   j=l and the scheme is thus unconditionally stable.
It can also be shown that (3.8) is a consistent dis- cretization of the differential BC's (2.7), (2.8).The details of the numerical analysis will be presented in [2].Fig. shows a simulation of a right traveling Gaussian beam evolving under the free Schr6dinger equation (h 1) with the rather course discretization Ax 0.00625, At 0.00002.Discretizing the analytic TBC's via (3.2) introduces strong numerical reflec- tions, while our discrete TBC's (3.8) recovers the smooth solution. 4. A TRANSIENT SCHRODINGER-POISSON SCATTERING MODEL WITH TBC'S TBC's for the 2D steady state SE have been used in [12] to model current carrying states in quantum devices.In 1D, these BC's are the steady state ana- logue of (2.7), (2.8).A similar ID Schr6dinger-Pois- son scattering model has been analyzed mathematically in [5].
The (transient) 1D TBC's of {}2 can of course only be used for simulating the time evolution of a pure state wave function .For the application to transient semiconductor device modeling, we will now present a self-consistent Schr6dinger-Poisson model for mixed quantum states with inhomogeneous TBC's.The (homogeneous) TBC (2,7) was derived for modeling the situation where an initial wave function is supported in the computational domain [0, L], and it is leaving this domain without being reflected back.
In our 1D transient scattering model the state of the electrons is represented by the sum of two density matrices p],2(x, y, t).P describes the contribution of the initial state Or, which is supported in the consid- ered domain [0, L]2, and P2 is the density matrix for the incoming wave packets from the boundary.
We remark that the scattering model of [5] is the steady state analogue of the above model (4.2)-(4.8).
The presented model was obtained in collaboration with N. Ben Abdallah, and a mathematical analysis of this model will be given in [3].