Fast Calculation of Energy and Mass
preserving solutions of Schr¨odinger–Poisson
systems on unbounded domains
Matthias Ehrhardt ∗,1and Andrea Zisowsky 1
Institut f¨ur Mathematik, Technische Universit¨at Berlin, Strasse des 17. Juni 136,
D–10623 Berlin, Germany
Abstract
This paper deals with the numerical solution of the time–dependent Schr¨odinger–
Poisson system in the spherically symmetric case. Since the problem is posed on an
unbounded domain one has to introduce artificial boundary conditions to confine
the computational domain. The main topic of this work is the construction of a
so–called discrete transparent boundary condition (TBC) for a Crank–Nicolson-
type predictor–corrector scheme for solving the Schr¨odinger–Poisson system. This
scheme has the property of mass and energy conservation exactly on the discrete
level. We propose different strategies for the discrete TBC and present an efficient
implementation. Finally, a numerical example illustrate the findings and shows the
comparison results between the different approaches.
Key words: Schr¨odinger–Poisson system, Schr¨odinger equation, finite differences,
discrete transparent boundary conditions, difference equation, Newton potential
PACS: 02.70.Bf, 31.15.Fx, 95.30.Cq
∗Corresponding author.
URLs: http://www.math.tu-berlin.de/~ehrhardt/ (Matthias Ehrhardt),
http://www.math.tu-berlin.de/~zisowsky/ (Andrea Zisowsky).
1Supported by the DFG Research Center “Mathematics for key technologies”
(FZT 86) in Berlin.
Preprint submitted to Journal of Computational and Applied Mathematics8 October 2004
1 Introduction
In many applications in quantum mechanics one wants to calculate the evo-
lution of an ensemble of particles over long time. This computations include
the solution of the single particle Schr¨odinger equation obtained from a mean
field approximation using Coulomb potentials [16]. The transient Schr¨odinger–
Poisson problem describes the time evolution of the wave function ψunder
the force of the self–consistent potential Vcaused by the charged electrons. It
is an appropriate model for semiconductor heterostructures (cf. [16] and the
references therein). We note that Schr¨odinger–Poisson systems appear in dif-
ferent applications, e.g. electron confinement in quantum nanostructures [14]
or as a description for the helium ground state in astrophysical applications
[19].
1.1 The Schr¨odinger–Poisson system. The transient Schr¨odinger–Poisson
system (SPS) associated with a single particle system in vacuum reads for the
complex–valued wave function ψ(x, t) and the electrostatic potential V(x, t):
i~∂tψ=−~2
2m∆xψ+V ψ, x ∈R3, t > 0,(1a)
∆xV=−γ n, x ∈R3, t > 0,(1b)
where n=|ψ(x, t)|2denotes the expected particle density for a pure quantum
state and γ > 0 (repulsive case) or γ < 0 (attractive case) depending on the
considered type of Coulomb force. Here ~denotes the Planck constant and m
is the particle mass. Throughout this paper we are interested in the attractive
case where the Schr¨odinger–Poisson system describes the time evolution of an
electron in a polar crystal (a polaron) under the assumption that the phonon
cloud or lattice vibrations behave classically. However, our results also hold in
the repulsive case. Equations (1) are supplied with some initial data ψ(x, 0) =
ψI(x) and the decay conditions
lim
|x|→∞ ψ(x, t) = 0,lim
|x|→∞ V(x, t) = 0.
The self–consistent potential Vcreated by the charged electrons is obtained
as solution to (1b), and can be written explicitly as
V(x, t) = γ
4πZR3
n(x0, t)
|x−x0|dx0.(2)
We remark that the existence and uniqueness of solutions to the Schr¨odinger–
Poisson system were analyzed in [7], [12]. Recently, a transient Schr¨odinger–
Poisson system in dimension 2 or 3 with transparent boundaries was stud-
ied [5] (the one–dimensional case was already treated in [4]). Finally, we note
2
that the Schr¨odinger–Poisson system is sometimes called Schr¨odinger–Newton
equations; the stationary spherical–symmetric case was investigated numeri-
cally in [18] and analytically in [24].
1.2 The spherically symmetric Schr¨odinger–Poisson system. Since we
want to keep the numerical effort to a minimum we only consider the case of
a spherically symmetric initial condition:ψ(x, 0) = ψI(r). It can be shown
that ψ(x, t) is invariant under rotations and therefore a radial function at any
time. For convenience we introduce the reduced wave function u(r, t) by
ψ(x, t) = 1
√4π
u(r, t)
r,(3)
and define the effective charge φ(r, t) = rV (x, t) which becomes using (2)
φ(r, t) = γ
4πZr
0|u(r0, t)|2dr0+Z∞
r
r
r0|u(r0, t)|2dr0.(4)
Differentiating (4) twice with respect to the radial coordinate r, the SPS
reduces then to
i~∂tu=−~2
2m∂2
ru+φ
ru, r > 0, t > 0,(5a)
∂2
rφ=−γ
4π|u|2
r, r > 0, t > 0,(5b)
together with the homogeneous Dirichlet conditions at the origin
u(0, t) = 0, φ(0, t) = 0,
and the decay conditions
lim
r→∞ u(r, t) = 0,lim
r→∞ φ(r, t) = γ
4π.
1.3 The conserved Quantities. The practically most important conserved
quantities usually are the mass of particles and the total energy. The mass of
particles is simply the L2–norm of uand therefore we can choose as normal-
ization condition Z∞
0|u(r, t)|2dr = 1, t > 0.(6)
The conserved total energy is given by
E(t) = EKIN(t) + EINT(t) + EPOT(t), t > 0,(7)
where the kinetic, interaction and potential energies are
EKIN(t) = ~2
2mZ∞
0|∂ru(r, t)|2dr, t > 0,(8a)
3
EINT(t) = Z∞
0
φ(r, t)
r|u(r, t)|2dr, t > 0,(8b)
EPOT(t) = 2π
γZ∞
0∂rφ(r, t)2dr, t > 0.(8c)
When solving the SPS numerically it is desirable that the discrete L2–norm
(mass of particles) and the discrete total energy are preserved exactly by the
numerical scheme because discretization errors in the conservation laws accu-
mulate. While it is easy to conserve the discrete L2–norm (yielding an unitary
evolution) it turns out that it is not so trivial to preserve simultaneously the
total energy
E(t) = ~2
2mZ∞
0|∂ru(r, t)|2dr +1
2Z∞
0
φ(r, t)
r|u(r, t)|2dr, t > 0.(9)
In this paper we consider the second order predictor–corrector scheme of
Ringhofer and Soler [20] where the discretization is based on the Crank–
Nicolson scheme. The conservation of mass and energy is achieved by intro-
ducing a phase modulation in the corrector step.
Remark 1 We remark that one can get an arbitrarily high (even) order mass–
conservative scheme for the Schr¨odinger equation by using the diagonal Pad´e
approximations to the exponential [3]. The Crank–Nicolson scheme corresponds
to second order, and the fourth order is known in the ODE literature as Ham-
mer and Hollingsworth method [11].
The outline of the paper is as follows. First we review briefly in §2 the ap-
proach of Ringhofer and Soler since their presentation in [20] is rather abstract
and details concerning a concrete implementation, e.g. appropriate boundary
conditions, are omitted. Afterwards in §3 we construct a so–called discrete
transparent boundary condition (DTBC) for a Schr¨odinger equation with a
Newton–type potential term, and discuss different approaches to obtain dis-
crete asymptotic solutions. We present in §4 an efficient implementation by
the sum–of–exponentials ansatz that reduces the computational effort for im-
plementing the DTBC significantly. Finally, we illustrate the results with a
numerical example comparing the different proposed versions of the DTBC.
2 The Numerical Schemes
In this section we first review the nonlinear Crank–Nicolson scheme from
[20] and show its conservation properties. Afterwards we turn to the linear
predictor–corrector scheme proposed by Ringhofer and Soler and analyze its
conservation properties, too. For simplicity of the presentation we will only
4
deal with uniform grids: (for a nonuniform time discretization and a general
spatial discretization cf. [20]):
u(n)
j∼u(rj, tn), φ(n)
j∼φ(rj, tn), rj=j∆r, tn=n∆t,
with 0 ≤j≤J,n≥0.
2.1 The nonlinear Crank–Nicolson scheme. The discretized SPS reads
i~D+
tu(n)
j=−~2
2mD2
ru(n+1
2)
j+φ(n+1
2)
j
rj
u(n+1
2)
j, j ≥1,(10a)
D2
rφ(n+1)
j=−γ
4πu(n+1)
j2
rj
, j ≥1,(10b)
together with the discrete boundary conditions
u(n)
0= 0,lim
j→∞ u(n)
j= 0, φ(n)
0= 0, φ(n)
J=γ
4π.(10c)
A realization of the decay condition for u(n)
jin the form of a (discrete) trans-
parent boundary condition will be the topic of §3. In (10) we have used the
standard abbreviations for the forward, and second order difference quotient:
D+
tu(n)
j=u(n+1)
j−u(n)
j
∆t, D2
ru(n)
j=u(n)
j+1 −2u(n)
j+u(n)
j−1
(∆r)2,
and the time averaging u(n+1
2)
j= (u(n+1)
j+u(n)
j)/2.
Discrete conservation properties of the nonlinear scheme. First we
want to show the discrete conservation of the mass (cf. (6))
ku(n)k2
2:= ∆r
∞
X
j=0 u(n)
j2.(11)
Here and in the sequel we will need a simple identity (“discrete product rule”)
D+
t(u(n)w(n)) = u(n+1
2)D+
t(w(n)) + w(n+1
2)D+
t(u(n)),(12)
i.e. with w(n)= ¯u(n)we get
D+
tu(n)2= 2 Re{¯u(n+1
2)D+
tu(n)}.(13)
To derive the discrete mass conservation property we consider the discretized
SPS (10) on the finite domain 1 ≤j≤J−1 and multiply (10a) with −i¯u(n+1
2)
j:
~¯u(n+1
2)
jD+
tu(n)
j=i~2
2m¯u(n+1
2)
jD2
ru(n+1
2)
j−iφ(n+1
2)
j
rju(n+1
2)
j
2
.(14)
5
Loading more pages...