A Solution representation by variation-of-constant formulas

# Convergence of a Strang splitting finite element discretization for the Schrödinger-Poisson equation

## Abstract

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.

Nonlinear Schrödinger equations, Operator splitting methods, Finite element discretization, Stability, Local error, Convergence

2 \subjclass65J15, 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 [[35]].

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.

#### 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 time-dependent Gross–Pitaevskii equations is given in [[2]], 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 [[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 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.

#### 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 multi-configuration time-dependent Hartree–Fock method, a complete convergence analysis of high-order 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 time-dependent 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 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

#### Schrödinger-Poisson equation

We consider the time-dependent Schrödinger-Poisson equation for ,

 i∂tψ(x,t)=−12Δψ(x,t)+Δ−1(|ψ(x,t)|2)ψ(x,t), (2.1a) where Ω⊂Rd, d∈{2,3}, is a bounded domain with smooth boundary. We impose homogeneous Dirichlet boundary conditions and an initial condition ψ(x,t)∣∣x∈∂Ω=0,ψ(x,0)=ψ0(x). (2.1b) For the subsequent analysis we will assume that the initial state satisfies2 ψ0∈H2=H2(Ω). The nonlocal nonlinear term Δ−1(|ψ|2) describing the electrostatic self-interaction is the solution Θ of the Poisson equation under homogeneous Dirichlet boundary conditions, ΔΘ(x,t)=|ψ(x,t)|2,Θ(x,t)∣∣x∈∂Ω=0. (2.1c)

The evolution operator associated with problem (2.1) will be denoted by , i.e.,

 ψ(⋅,t)=φ\tiny{SP}(t,ψ0).

#### Abstract formulation

Introducing the operator notation

 A:H2∩H10→L2:\leavevmode\nobreak \leavevmode\nobreak u↦12iΔu,ˆB:H10→H2∩H10:\leavevmode\nobreak \leavevmode\nobreak w↦−iΔ−1(|w|2),B:H10→H10:\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak u↦−iΔ−1(|u|2)u, (2.2a) we employ a compact formulation of problem (2.1) as an abstract evolution equation {\leavevmode\nobreak ∂tψ=Aψ+B(ψ)=Aψ+ˆB(ψ)ψ,\leavevmode\nobreak ψ∣∣t=0=ψ0. (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

 {\leavevmode\nobreak ∂tψ=Aψ,\leavevmode\nobreak ψ∣∣t=0=u, (2.3a) is denoted by φA(t), such that ψ(⋅,t)=φA(t)u. (2.3b)
• The evolution operator associated with the nonlinear initial value problem

 {\leavevmode\nobreak ∂tψ=B(ψ),\leavevmode\nobreak ψ∣∣t=0=u, (2.4a) is denoted by φB(t,⋅), such that ψ(⋅,t)=φB(t,u). (2.4b) Due to the fact that Δ−1(|ψ(⋅,t)|2) defines a real-valued function and thus ∂t|ψ(⋅,t)|2=2R(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ψ(⋅,t)∂tψ(⋅,t))=2R(ˆB(ψ(⋅,t))|ψ(⋅,t)|2)=0, the nonlinear equation (2.4a) reduces to the linear equation {\leavevmode\nobreak ∂tψ=ˆB(u)ψ,\leavevmode\nobreak ψ∣∣t=0=u. (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

 {\leavevmode\nobreak ∂tψ=ˆB(w)ψ,\leavevmode\nobreak ψ∣∣t=0=u, (2.5)

such that . Clearly,

 φB(t,u)=EB(t,u)u. (2.6)

#### Strang splitting method

Our main focus is on the symmetric second-order Strang splitting method applied to the splitting according to (2.3), (2.4). That is, for a time increment , the time-discrete solution values

 ψn≈ψ(nτ),n=0,1,2,…

are determined by the recurrence

 ψn=S(τ,ψn−1)=φA(12τ)φB(τ,φA(12τ)ψn−1). (2.7a) For notational simplicity we shall employ a formal notation for the n-fold composition, ψn=Snψ0:=S(τ,⋅)∘⋯∘S(τ,ψ0)n times. (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),

 {\leavevmode\nobreak (∂tψ,ϕ)L2=−12i(∇ψ,∇ϕ)L2for all\leavevmode\nobreak \leavevmode\nobreak ϕ∈H10,\leavevmode\nobreak ψ∣∣t=0=u, (2.8)

where we require .

 For (2.4c), {\leavevmode\nobreak (∂tψ,ϕ)L2=−i(Θψ,ϕ)L2for all\leavevmode\nobreak \leavevmode\nobreak ϕ∈H10,\leavevmode\nobreak ψ∣∣t=0=u, (2.9a) where Θ is the solution of the Poisson equation in weak formulation, (∇Θ,∇χ)L2=−(|u|2,χ)L2for all\leavevmode\nobreak \leavevmode\nobreak χ∈H10, (2.9b)

requiring .

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

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

 Ω=K⋃k=1Ωk,h=maxk∈{1,…,K}diamΩk,

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

 Vh=span{v1,…,vJ}⊂H10. (2.10)

#### Finite element interpolation and projection

By we denote the nodal interpolation operator,

 Ih(f)=J∑j=1f(xj)vj. (2.11)

The Rayleigh-Ritz projection is defined implicitly by the Galerkin orthogonality relation

 (∇(u−Phu),∇vh)L2=0for all\leavevmode\nobreak \leavevmode\nobreak vh∈Vh, (2.12a) satisfying |Phu|H1≤|u|H1for all\leavevmode\nobreak \leavevmode\nobreak u∈H10. (2.12b) By the Poincaré inequality ∥v∥L2≤C|v|H1, this also implies ∥Phu∥H1≤C∥u∥H1for all\leavevmode\nobreak \leavevmode\nobreak u∈H10, (2.12c) with a constant C depending on Ω.
###### Remark 2.1.

The Rayleigh-Ritz projection is connected to the finite element approximation in the following way. Consider a Poisson problem with homogeneous Dirichlet boundary conditions in weak formulation (see (2.9b)),

 (∇u,∇v)L2=−(f,v)L2for all\leavevmode\nobreak \leavevmode\nobreak v∈H10,

and its FEM discretization by the Galerkin equations (see (2.15b) below),

 (∇uh,∇vh)L2=−(f,vh)L2for all\leavevmode\nobreak \leavevmode\nobreak vh∈Vh.

Then,

 (∇(u−uh),∇vh)L2=0for all\leavevmode\nobreak \leavevmode\nobreak vh∈Vh,

i.e.,

 uh=Phu.

For a sufficiently smooth boundary, the regularity estimate

 |Δ−1(u−Phu)|H2≤C∥u−Phu∥L2 (2.13)

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,

 Ih(ψ0)=J∑j=1cj,0vj∈Vh.

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.

•  For the first subproblem (2.8), ψh is determined from {\leavevmode\nobreak (∂tψh,ϕh)L2=−12i(∇ψh,∇ϕh)L2% for all\leavevmode\nobreak \leavevmode\nobreak ϕh∈Vh,\leavevmode\nobreak ψh∣∣t=0=uh. (2.14a) With the ansatz in terms of the basis (2.10), ψh(t)=J∑j=1cj(t)vj, (2.14a) yields the Galerkin equations for the coefficients cj(t) in the form J∑j=1∂tcj(t)(vj,vi)L2=−12iJ∑j=1cj(t)(∇vj,∇vi)L2,i=1…J. (2.14b)
• For the second subproblem (2.9), is determined such that

 {\leavevmode\nobreak (∂tψh,ϕh)L2=−i(Θhψh,ϕh)L2for all\leavevmode\nobreak \leavevmode\nobreak ϕh∈Vh,\leavevmode\nobreak ψh∣∣t=0=uh, (2.15a) where Θh is the solution of the discretized Poisson problem (∇Θh,∇χh)L2=−(|uh|2,χh)L2for all\leavevmode\nobreak \leavevmode\nobreak χh∈Vh. (2.15b) With the ansatz in terms of the basis (2.10), ψh(t)=J∑j=1cj(t)vj∈Vh,Θh=J∑j=1djvj∈Vh, we obtain the Galerkin equations (2.15b) in the form J∑j=1∂tcj(t)(vj,vi)L2=−iJ∑j=1cj(t)J∑k=1dk(vkvj,vi)L2,i=1…J,J∑j=1dj(∇vj,∇vi)L2=−J∑j,k=1cj(0)¯¯¯¯¯¯¯¯¯¯¯¯ck(0)(vjvk,vi)L2,i=1…J. (2.15c)

The above computations imply that the unknown coefficients satisfy systems of linear ordinary differential equations. For a more compact formulation we introduce the vectors

 c(t)=(cj(t))Jj=1∈CJ,d=(dj)Jj=1∈RJ,F(c(0))=(J∑j,k=1cj(0)¯¯¯¯¯¯¯¯¯¯¯¯ck(0)(vjvk,vi)L2)Ji=1∈RJ, (2.16a) and the invertible symmetric matrices M =(Mij)Ji,j=1∈RJ×J,with\leavevmode\nobreak \leavevmode\nobreak Mij=(vi,vj)L2, (2.16b) K =(Kij)Ji,j=1∈RJ×J,with\leavevmode\nobreak \leavevmode\nobreak Kij=(∇vi,∇vj)L2, Φ(d) =(Φij(d))Ji,j=1∈RJ×J,with\leavevmode\nobreak \leavevmode\nobreak Φij(d)=∑kdk(vkvi,vj)L2.

In this notation, the system (2.14b) reads

 {\leavevmode\nobreak M∂tc(t)=−12iKc(t),\leavevmode\nobreak c(0)given, (2.17)

with solution

 c(t)=e−12itM−1Kc(0). (2.18)

System (2.15c) takes the form

 ⎧⎪⎨⎪⎩\leavevmode\nobreak Kd=−F(c(0)),\leavevmode\nobreak M∂tc(t)=−iΦ(d)c(t),\leavevmode\nobreak c(0)given, (2.19)

with solution

 c(t)=e−it(M−1Φ(d))c(0), where\leavevmode\nobreak \leavevmode\nobreak d=−K−1F(c(0)). (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

 (Δhuh,vh)L2=−(∇uh,∇vh)L2for all\leavevmode\nobreak \leavevmode\nobreak vh∈Vh, (2.21a) (∇Δ−1hf,∇vh)L2=−(f,vh)L2\leavevmode\nobreak for all\leavevmode\nobreak \leavevmode\nobreak vh∈Vh. (2.21b)

In particular, (2.21b) means that for the solution of we have

 Δ−1hf=Phu=uh (2.22)

in the sense of Remark 2.1.

Moreover, in analogy to (2.2a) we set

 Ah:Vh→H−1:\leavevmode\nobreak \leavevmode\nobreak uh↦12iΔhuh,ˆBh:Vh→Vh:\leavevmode\nobreak \leavevmode\nobreak wh↦−iΔ−1h(|wh|2),Bh:Vh→Vh:\leavevmode\nobreak \leavevmode\nobreak uh↦−iΔ−1h(|uh|2)uh. (2.23)

In this notation,

• is associated with subproblem (2.14a),

• is associated with subproblem (2.15a).

For representing the solution of a system of the type

 {\leavevmode\nobreak (∂tψh,ϕh)L2=(ˆBh(wh)ψh,ϕh)L2for all\leavevmode\nobreak \leavevmode\nobreak ϕh∈Vh,\leavevmode\nobreak ψh∣∣t=0=uh, (2.24)

we will also employ an analogous notation as for problem (2.5),

 ψh=EBh(t,wh)uh. (2.25a) Then, analogously as in (2.6), φBh(t,uh)=EBh(t,uh)uh. (2.25b)

For the resulting fully discrete Strang splitting solution we again write

 ψn≈ψ(nτ),n=0,1,2,…,

determined by the recurrence

 Missing or unrecognized delimiter for \big (2.26a) and we again employ a formal notation for the n-fold composition, ψn=Snhψ0:=Sh(τ,⋅)∘⋯∘Sh(τ,ψ0)n times. (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

 ψn−ψ(tn)=(SnhIhψ0−Snψ0)+(Snψ0−φ\tiny{SP}(tn,ψ0)). (2.27)

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)
• The second term in (2.27) can similarly be recast as

 Snψ0−φ\tiny{SP}(tn,ψ0)=n∑j=1Sn−j(Sφ\tiny{SP}(tj−1,ψ0)−φ\tiny{SP}(τ,φ\tiny{SP}(tj−1,ψ0)). (2.28b)
• Here apply a standard argument for estimating the splitting error at the semi-discrete level combining the stability of the splitting operator (see Sec. 3.2) with an estimate for the local splitting error (see Theorem 3.2 or [[27]]).

This leads to the following global error bound for the full discretization.

###### Theorem 2.2.

Suppose that , and that is such that (2.13) holds. Consider the fully discretized method from (2.26a) based on the Strang splitting scheme and conforming finite elements of degree , then

 Missing or unrecognized delimiter for \big (2.29)

where and . Here, depends on , , the - and the -norms of .

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 ,

 ∥ψn−ψ(tn)∥L2≤Ctn(τ2+hp+1),∥ψn−ψ(tn)∥H1≤Ctn(τ+hp).
• 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 ),

 ∥ψn−ψ(tn)∥L2≤Ctn(τ2+hsτ),∥ψn−ψ(tn)∥H1≤Ctn(τ+hs−1τ),

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.23.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 semi-discrete error can be bounded in and :

 ∥SnhIhψ0−Snψ0∥L2 ≤Ctnhs(1+1τ)β, ∥SnhIhψ0−Snψ0∥H1 ≤Ctnhs−1(1+1τ)β,

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):

 ∥SnhIhψ0−Snψ0∥L2≤∥Snh(Id−Ih)ψ0∥L2+∥Snh(Id−Ph)ψ0∥L2 +n∑j=1∥Sn−jhShPhSj−1ψ0−Sn−jhPhSSj−1ψ0∥L2+∥(Id−Ph)Snψ0∥L2 ≤eCtna2Sh(∥(Id−Ih)ψ0∥L2+∥(Id−Ph)ψ0∥L2) Extra open brace or missing close brace ≤eCtna2Sh(∥(Id−Ih)ψ0∥L2+∥(Id−Ph)ψ0∥L2) (3.1a) Extra open brace or missing close brace (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 .

The expressions in (3.1a) can be bounded using Theorems C.4 and C.5,

 ∥(Id−Ph)u∥L2 ≤Chs∥u∥Hs, ∥(Id−Ih)u∥L2 ≤Chs∥u∥Hs,

and by the bound (C.12) from Proposition C.12 we obtain

 ∥Snψ0∥Hs≤eLstn∥ψ0∥Hs.

It remains to bound . By Theorem 3.7 from Sec. 3.3, it follows that

 ∥Sh(τ,Phu)−PhS(τ,u)∥L2≤Cτ(1+1τ)βhs,

for , where depends on , , , and . Altogether, we obtain

 ∥SnhIhψ0−Snψ0∥L2 ≤eCtna2Shhs∥ψ0∥Hs+etnLshs∥ψ0∥Hs+neCtna2ShCτ(1+1τ)βhs ≤Chs(eCtna2Sh∥ψ0∥Hs+tn(1+1τ)βeCtna2Sh) ≤Ctn(1+1τ)βhs,

which concludes the proof for the bound.

• -bound. Analogously as for the -bound, we use Theorems C.4 and C.5 and Proposition C.12 and Theorem 3.7. Hence we obtain

 ∥SnhIhψ0−Snψ0∥H1 ≤∥SnhIhψ0−Snhψ0∥H1+∥Snhψ0−SnhPhψ0∥H1 +n∑j=1∥∥Sn−j</