Convergence of a Strang splitting finite element discretization for the Schrödinger-Poisson equation
Operator splitting methods combined with finite element spatial discretizations are studied for time-dependent 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 second-order Strang splitting method and conforming polynomial finite element discretizations. For sufficiently regular solutions the classical orders of convergence are retained, that is, second-order 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, Convergence
1991 Mathematics Subject Classification:65J15, 65L05, 65M60 65M12, 65M15
1. Introduction and overview
We consider full discretization methods for the time-dependent Schrödinger–Poisson equation, which typically arises in models of quantum transport [10, 20]. Our approach relies on a second-order 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 solution-adapted non-uniform spatial grid, which can be updated in the course of the time integration .
Our main objective is to provide an error analysis for this full discretization, showing the expected second-order 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 right-hand side are of higher order than the discretization error and therefore will not be taken into account.
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 time-dependent Gross–Pitaevskii equations is given in , which summarizes most of the studies conducted in this field. The Crank-Nicholson 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. Semi-implicit 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. Semi-implicit finite difference schemes lose most of the desired properties. For regular solutions, time-splitting methods in conjunction with Fourier- or Sine-spectral 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 time-stepping schemes. For non-smooth 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 ).
Recently, full discretization of the Schrödinger–Poisson equation by splitting methods in conjunction with spectral space discretization has been investigated in , where the long-range interaction is approximated efficiently by nonuniform fast Fourier transform (NUFFT). The authors conclude superior accuracy and performance of their approach in particular over the Sine-spectral method.
The stability and error behavior of operator splitting methods for the Schrödinger–Poisson equation have first been analyzed in . For the structurally similar equations associated with the multi-configuration time-dependent Hartree–Fock method, a complete convergence analysis of high-order splitting methods has been given in . An error analysis of splitting methods applied to the Schrödinger–Poisson equation in the semiclassical regime is provided in .
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 time-dependent Schrödinger equations see for example [37, 21] and the more recent contribution , and for atomic and molecular systems see  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 . The investigation in the context of the Schrödinger–Poisson equation remains an open question.
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 higher-order 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
We consider the time-dependent Schrödinger-Poisson equation for ,
|where , , is a bounded domain with smooth boundary.|
We impose homogeneous Dirichlet boundary conditions and an initial condition
For the subsequent analysis we will assume that the initial state satisfies111For simplicity of notation we write instead of , etc. . The nonlocal nonlinear term describing the electrostatic self-interaction is the solution of the Poisson equation under homogeneous Dirichlet boundary conditions,
The evolution operator associated with problem (2.1) will be denoted by , i.e.,
Introducing the operator notation
|we employ a compact formulation of problem (2.1) as an abstract evolution equation|
2.2. Semidiscretization in time by the Strang splitting method
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 real-valued 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
such that . Clearly,
Strang splitting method
are determined by the recurrence
|For notational simplicity we shall employ a formal notation for the -fold composition,|
Weak formulation of the subproblems
where we require .
|where is the solution of the Poisson equation in weak formulation,|
In the following we use the standard denotation for the Sobolev semi-norms, i.e., for , and for .
2.3. Conforming finite element discretization of the subproblems
Finite element space
For the space discretization of the subproblems, we choose a tessellation over subdomains , with
which are affine-equivalent 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 quasi-uniform. As common we choose a linear indexing of the basis functions, . The subspace spanned by these functions is denoted by
Finite element interpolation and projection
By we denote the nodal interpolation operator,
The Rayleigh-Ritz projection is defined implicitly by the Galerkin orthogonality relation
|By the Poincaré inequality , this also implies|
|with a constant depending on .|
For a sufficiently smooth boundary, the regularity estimate
holds, see [9, Sec. 5.5].
Fully discrete solution and computational representation
The full discretization of the Schrödinger-Poisson 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,
The above computations imply that the unknown coefficients satisfy systems of linear ordinary differential equations. For a more compact formulation we introduce the vectors
|and the invertible symmetric matrices|
In this notation, the system (2.14b) reads
System (2.15c) takes the form
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
In particular, (2.21b) means that for the solution of we have
in the sense of Remark 2.1.
Moreover, in analogy to (2.2a) we set
In this notation,
For representing the solution of a system of the type
we will also employ an analogous notation as for problem (2.5),
|Then, analogously as in (2.6),|
For the resulting fully discrete Strang splitting solution we again write
determined by the recurrence
|and we again employ a formal notation for the -fold composition,|
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
The first term represents the error attributable to the space discretization and the second term is the splitting error at the semi-discrete 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 finite-element interpolants (see Theorem C.4) and the Rayleigh-Ritz-projection (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.
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 ),
3. Convergence analysis
3.1. Global error bound
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  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.
Let for , , and . Then for , the and bounds of the semi-discrete error can be bounded in and :
where and depends on , , , and .