Convergence of a Strang splitting finite element discretization for the SchrödingerPoisson equation
Abstract.
Operator splitting methods combined with finite element spatial discretizations are studied for timedependent nonlinear Schrödinger equations. In particular, the Schrödinger–Poisson equation under homogeneous Dirichlet boundary conditions on a finite domain is considered. A rigorous stability and error analysis is carried out for the secondorder Strang splitting method and conforming polynomial finite element discretizations. For sufficiently regular solutions the classical orders of convergence are retained, that is, secondorder convergence in time and polynomial convergence in space is proven. The established convergence result is confirmed and complemented by numerical illustrations.
Key words and phrases:
Nonlinear Schrödinger equations, Operator splitting methods, Finite element discretization, Stability, Local error, Convergence1991 Mathematics Subject Classification:
65J15, 65L05, 65M60 65M12, 65M152
1. Introduction and overview
We consider full discretization methods for the timedependent Schrödinger–Poisson equation, which typically arises in models of quantum transport [10, 20]. Our approach relies on a secondorder Strang splitting time discretization combined with a conforming finite element space discretization. The motivation for the proposed solution method is that separate treatment of the nonlinear part suggests the application of special solvers for the Poisson equation, which are particularly efficient in the context of an underlying finite element space discretization. For this purpose it is common to truncate the unbounded spatial domain to a sufficiently large finite domain and impose homogeneous Dirichlet boundary conditions. Indeed, the evaluation of the nonlocal convolution integral in the standard formulation generally implies a huge computational effort caused by the suitable treatment of the singular integral kernel for the evaluation on a large domain. By the splitting approach, we can separately treat the Poisson equation by appropriate methods where optimized linear solvers are available as for instance multigrid or domain decomposition methods [31, 33]. The finite element discretization additionally enables a solution on a solutionadapted nonuniform spatial grid, which can be updated in the course of the time integration [35].
Our main objective is to provide an error analysis for this full discretization, showing the expected secondorder convergence of the Strang splitting method and polynomial spatial error decay corresponding with the finite elements employed. By using Gauss–Lobatto nodes, the setup of the stiffness matrix is exact; the errors arising in the construction of the mass matrix and the righthand side are of higher order than the discretization error and therefore will not be taken into account.
Splitting methods
The computational advantages of operator splitting methods for the time integration of problems in quantum dynamics have been emphasized in recent literature. A comprehensive overview of investigations for timedependent Gross–Pitaevskii equations is given in [2], which summarizes most of the studies conducted in this field. The CrankNicholson finite difference method preserves most of the important invariants like symmetry in time, mass and energy and is unconditionally stable; however, the computational cost for this fully implicit method is considerable, and the conservation properties only hold up to the accuracy of the nonlinear solver. Semiimplicit relaxation methods which only treat the kinetic part implicitly share the conservation properties if only a cubic nonlinearity is present, but they are still computationally expensive and suffer from stability limitations. Semiimplicit finite difference schemes lose most of the desired properties. For regular solutions, timesplitting methods in conjunction with Fourier or Sinespectral methods are overall concluded to be the most successful discretization schemes; they are unconditionally stable, conserve norm, energy, and also dispersion, which is not the case for many other timestepping schemes. For nonsmooth or random spatial profiles, the spectral accuracy may be lost, however, and thus splitting methods in conjunction with finite difference spatial discretizations may be more efficient (see [7]).
Recently, full discretization of the Schrödinger–Poisson equation by splitting methods in conjunction with spectral space discretization has been investigated in [8], where the longrange interaction is approximated efficiently by nonuniform fast Fourier transform (NUFFT). The authors conclude superior accuracy and performance of their approach in particular over the Sinespectral method.
Error analysis
The stability and error behavior of operator splitting methods for the Schrödinger–Poisson equation have first been analyzed in [27]. For the structurally similar equations associated with the multiconfiguration timedependent Hartree–Fock method, a complete convergence analysis of highorder splitting methods has been given in [25]. An error analysis of splitting methods applied to the Schrödinger–Poisson equation in the semiclassical regime is provided in [11].
Finite element method
The literature on finite element spatial discretizations is vast. Finite element methods (FEM) have been widely used for electronic structure calculations, see for instance [14, 34] and [6, 13, 30, 38]. For the solution of timedependent Schrödinger equations see for example [37, 21] and the more recent contribution [22], and for atomic and molecular systems see [19] for a general review.
Truncation to a finite domain
In conjunction with the application of the finite element method, the restriction to a finite domain introduces a truncation error which we do not consider in this work. Strategies to cope with related issues have been proposed for instance in [3]. The investigation in the context of the Schrödinger–Poisson equation remains an open question.
Outline
In Sec. 2 we state the Schrödinger–Poisson equation. We specify the full discretization method and formulate our main convergence results. In Sec. 3 we provide the underlying comprehensive stability and error analysis. Our numerical illustrations given in Sec. 4 confirm the theoretical convergence result and demonstrate that also higherorder splitting methods show their expected behavior. The appendices contain proof details, important results from the literature which we rely on and auxiliary estimates used in our analysis.
2. Problem setting, discretization method, and main results
2.1. Problem setting
SchrödingerPoisson equation
We consider the timedependent SchrödingerPoisson equation for ,
(2.1a)  
where , , is a bounded domain with smooth boundary. 
We impose homogeneous Dirichlet boundary conditions and an initial condition
(2.1b) 
For the subsequent analysis we will assume that the initial state satisfies^{1}^{1}1For simplicity of notation we write instead of , etc. . The nonlocal nonlinear term describing the electrostatic selfinteraction is the solution of the Poisson equation under homogeneous Dirichlet boundary conditions,
(2.1c) 
The evolution operator associated with problem (2.1) will be denoted by , i.e.,
Abstract formulation
Introducing the operator notation
(2.2a)  
we employ a compact formulation of problem (2.1) as an abstract evolution equation  
(2.2b) 
2.2. Semidiscretization in time by the Strang splitting method
Subproblems
For the discretization of (2.2b) in time we apply exponential operator splitting methods based on the solution of two subproblems, see for instance [17, 28].

The evolution operator associated with the linear initial value problem
(2.3a) is denoted by , such that (2.3b) 
The evolution operator associated with the nonlinear initial value problem
(2.4a) is denoted by , such that (2.4b) Due to the fact that defines a realvalued function and thus the nonlinear equation (2.4a) reduces to the linear equation (2.4c) We will also employ a notation analogous to (2.4c) but with a linear evolution operator depending on and as the solution to
(2.5) such that . Clearly,
(2.6)
Strang splitting method
Our main focus is on the symmetric secondorder Strang splitting method applied to the splitting according to (2.3), (2.4). That is, for a time increment , the timediscrete solution values
are determined by the recurrence
(2.7a)  
For notational simplicity we shall employ a formal notation for the fold composition,  
(2.7b) 
Weak formulation of the subproblems
In view of full discretization (see Sec. 2.3) we consider the following weak formulations of the subproblems. For (2.3a),
(2.8) 
where we require .
For (2.4c),  
(2.9a)  
where is the solution of the Poisson equation in weak formulation,  
(2.9b) 
requiring .
In the following we use the standard denotation for the Sobolev seminorms, i.e., for , and for .
2.3. Conforming finite element discretization of the subproblems
A full discretization arises by solving both initial value subproblems (2.3) and (2.4) in their weak reformulation (2.8) and (2.9), respectively, by means of a finite element method (FEM).
Finite element space
For the space discretization of the subproblems, we choose a tessellation over subdomains , with
which are affineequivalent to a reference domain . With we associate a triplet , where the set comprises the interpolation nodes , and is the linear space spanned by the polynomial nodal basis functions of degree . We require the finite elements to be conforming and quasiuniform. As common we choose a linear indexing of the basis functions, . The subspace spanned by these functions is denoted by
(2.10) 
Finite element interpolation and projection
By we denote the nodal interpolation operator,
(2.11) 
The RayleighRitz projection is defined implicitly by the Galerkin orthogonality relation
(2.12a)  
satisfying  
(2.12b)  
By the Poincaré inequality , this also implies  
(2.12c)  
with a constant depending on . 
Remark 2.1.
Fully discrete solution and computational representation
The full discretization of the SchrödingerPoisson equation is based on solving the subproblems (2.8), (2.9) arising in the Strang splitting time discretization by means of a FEM/Galerkin space discretization. Here, the coefficients associated with the prescribed initial state are determined by interpolation,
In each substep of the time propagation by Strang splitting, subproblems of the following types (2.14), (2.15) for the solutions and are solved.
The above computations imply that the unknown coefficients satisfy systems of linear ordinary differential equations. For a more compact formulation we introduce the vectors
(2.16a)  
and the invertible symmetric matrices  
(2.16b)  
In this notation, the system (2.14b) reads
(2.17) 
with solution
(2.18) 
System (2.15c) takes the form
(2.19) 
with solution
(2.20) 
To realize the fully discrete propagation in time according to the Strang recurrence (2.7), systems of this type are alternately solved.
Finite element operators
We define the discrete Laplace operator () and its inverse () via
(2.21a)  
(2.21b) 
In particular, (2.21b) means that for the solution of we have
(2.22) 
in the sense of Remark 2.1.
Moreover, in analogy to (2.2a) we set
(2.23) 
In this notation,
For representing the solution of a system of the type
(2.24) 
we will also employ an analogous notation as for problem (2.5),
(2.25a)  
Then, analogously as in (2.6),  
(2.25b) 
For the resulting fully discrete Strang splitting solution we again write
determined by the recurrence
(2.26a)  
and we again employ a formal notation for the fold composition,  
(2.26b) 
2.4. Main results
The central interest of this paper is to establish a convergence result for the splitting finite element discretization of the Schrödinger–Poisson equation (2.1). Here we give a brief overview of the structure of our convergence proof and state the resulting theorem. The detailed convergence analysis is worked out in Sec. 3.
In order to study the global error we separate the terms associated with space and time discretization, respectively. With and , we write
(2.27) 
The first term represents the error attributable to the space discretization and the second term is the splitting error at the semidiscrete level.

The first term in (2.27) is expanded into a telescoping sum in the following way:
(2.28a) We combine a stability argument for the fully discrete splitting operator (see Sec. 3.2) with the approximation properties of the finiteelement interpolants (see Theorem C.4) and the RayleighRitzprojection (see Theorem C.5). What remains to be estimated are terms of the form , which is worked out in Sec. 3.3 (see Theorem 3.1)
This leads to the following global error bound for the full discretization.
Theorem 2.2.
The proof of Theorem 2.2 is based on the combination of Theorem 3.1 for the contribution and on Theorem 3.2 for the contribution of .
Conclusions
From Theorem 2.2, we can deduce the following convergence properties:

For an initial value , we obtain the classical convergence order in and ,

For an initial value , we obtain convergence of order respectively in space, but with a possibly reduced convergence order in time (depending in the ratio between and ),
where .
3. Convergence analysis
3.1. Global error bound
We start by separating the effects of space and time discretization, see (2.28) and (2.28b), and consider bounds in the  and norm.
By a Lady Windermere’s fan argument and the stability estimates from Sec. 3.2, the expression in (2.28) can be expressed by an bound in and an bound in , as shown in the following Theorem 3.1. The norms have already been studied in [27] and are summarized in Theorem 3.2 below. This implies the main convergence result stated in Theorem 2.2, where error bounds depending on the regularity of the initial values are given.
In our convergence theory we make use of several stability estimates and consistency results which are collected in Sec. 3.2–3.4 below. Several auxiliary results and estimates are collected in the appendix.
Theorem 3.1.
Let for , , and . Then for , the and bounds of the semidiscrete error can be bounded in and :
where and depends on , , , and .
Proof.

bound. We proceed as indicated at the beginning of Sec. 2.4 and use the stability properties (3.4) of the splitting operator (see Proposition 3.3 in Sec. 3.2):
(3.1a) (3.1b) By the regularity result for the splitting operator , see Lemma 3.8 in Sec. 3.4, we can ensure the existence of the constant .