A Priori Error Analysis of Linear and Nonlinear Periodic Schrödinger Equations with Analytic Potentials

Abstract

This paper is concerned with the numerical analysis of linear and nonlinear Schrödinger equations with periodic analytic potentials. We prove that, for linear equations, when the potential is analytic in a strip of width A of the complex plane, the solution is analytic in the same strip, ensuring an exponential convergence of the planewave discretization of the equation with rate A. On the other hand, for nonlinear equations, we find that the solution may be analytic only in a strip of width smaller than A. This behavior is illustrated by two examples using a combination of numerical and analytical arguments.

Appendices

Appendix A: Extension to the Multidimensional Case with Application to Kohn–Sham Models

The goal of this section is to extend the previous results to the multidimensional case and apply them to the linear version of the Kohn–Sham equations (1). To this end, consider a Bravais lattice $\mathbb {L}=\mathbb {Z}{\varvec{a}}_1+ \cdots + \mathbb {Z}{\varvec{a}}_d$ where ${\varvec{a}}_1,\dots , {\varvec{a}}_d$ are linearly independent vectors of $\mathbb {R}^d$ ($d=3$ for KS-DFT). Up to an affine change of variables (which preserves analyticity), we can take without loss of generality $\mathbb {L}= 2\pi \mathbb {Z}^{d}$. We denote by $\Omega =[0,2\pi )^{d}$ a unit cell, by $\mathbb {L}^* = \mathbb {Z}^{d}$ the reciprocal lattice, by $e_{{\varvec{G}}}({\varvec{x}})=|\Omega |^{-1/2}{\textrm{e}}^{i{\varvec{G}} \cdot {\varvec{x}}}$ the Fourier mode with wavevector ${\varvec{G}} \in \mathbb {L}^*$, and by

$$\begin{aligned} H^s_{{\textrm{per}},\mathbb {L}} {:}{=}\left\{ u= \sum _{{\varvec{G}} \in \mathbb {L}^*} \widehat{u}_{{\varvec{G}}} e_{{\varvec{G}}} \, \bigg \vert \, \sum _{{\varvec{G}} \in \mathbb {L}^*} (1+|{\varvec{G}}|^2)^s |{\widehat{u}}_{{\varvec{G}}}|^2 < \infty \right\} \end{aligned}$$

the $\mathbb {L}$-periodic Sobolev spaces endowed with their usual inner products. All the arguments in Sects. 3.1-3.3 can be extended to the multidimensional case by introducing the Hilbert spaces

$$\begin{aligned} \mathcal {H}_{A,\mathbb {L}} {:}{=}\left\{ u \in L^2_{{\textrm{per}},\mathbb {L}} \,\bigg \vert \, \sum _{{\varvec{G}}\in \mathbb {L}} {w}_{A,\mathbb {L}} ({\varvec{G}}) \left| \widehat{u}_{{\varvec{G}}} \right| ^2 < \infty \right\} ,\quad (u,v)_{A,\mathbb {L}} {:}{=}\sum _{{\varvec{G}}\in \mathbb {L}^*} {w}_{A,\mathbb {L}}({\varvec{G}}) \overline{\widehat{u}_{{\varvec{G}}}} \widehat{v}_{{\varvec{G}}}, \end{aligned}$$

where $\displaystyle {{w}_{A,\mathbb {L}}({\varvec{G}}) = \sum _{n=1}^d w_A(G_{n})}$. Each $u \in \mathcal {H}_{A,\mathbb {L}}$ can be extended to an analytic function $u(z_1,\dots ,z_d)$ of d complex variables defined on a neighborhood on $\mathbb {R}^d$, and it holds

$$\begin{aligned} \sum _{{\varvec{G}}\in \mathbb {L}^*} w_{A,\mathbb {L}}({\varvec{G}}) \left| \widehat{u}_{{\varvec{G}}} \right| ^2 = \frac{1}{2} \sum _{j=1}^d \int _{\Omega } \left| {u}({\varvec{x}}+ \textrm{i} e_{j}) \right| ^2 + \left| {u}({\varvec{x}}- \textrm{i} e_{j}) \right| ^2 \textrm{d}{\textbf{x}}. \end{aligned}$$

with $e_{1},\dots , e_d$ the canonical basis vectors. The approximation space $X_{N,\mathbb {L}}$ is then defined as

$$\begin{aligned} X_{N,\mathbb {L}} {:}{=}\text{ Span } (e_{{\varvec{G}}}, \, {\varvec{G}} \in \mathbb {L}^*, \, |{{\varvec{G}}}| \leqslant N ), \end{aligned}$$

and the inverse $T_{22,\mathbb {L}}^{-1}$ of the restriction $T_{22,\mathbb {L}}$ of the operator $-\Delta$ on $L^2_{{\textrm{per}},\mathbb {L}}$ to the invariant subspace $X_{N,\mathbb {L}}^\perp = \text{ Span } (e_{{\varvec{G}}}, \, {\varvec{G}} \in \mathbb {L}^*, \, |{{\varvec{G}}}| > N )$ satisfies

$$\begin{aligned} \Vert T^{-1}_{22,\mathbb {L}} \Vert _{\mathcal {L}(X_{N,\mathbb {L}}^\perp )} = \Vert T^{-1}_{22,\mathbb {L}} \Vert _{\mathcal {L}(\mathcal {H}_{A,\mathbb {L}} \cap X_{N,\mathbb {L}}^\perp )} \leqslant N^{-2}. \end{aligned}$$

The proofs of Theorems 1, 2 and 3 can thus be straightforwardly adapted to the multidimensional case.

Lastly, if $V \in \mathcal {H}_{B,\mathbb {L}}$ for some $B > 0$, the Schrödinger operator $H = -\Delta + V$ considered this time as a Schrödinger operator on $L^2(\mathbb {R}^d,\mathbb {C})$ with an $\mathbb {L}$-periodic potential, can be decomposed by the Bloch transform [22, Section XIII.16] and its Bloch fibers are the self-adjoint operators on $L^2_{{\textrm{per}},\mathbb {L}}$ with domain $H^2_{{\textrm{per}},\mathbb {L}}$ and form domain $H^1_{{\textrm{per}},\mathbb {L}}$ defined as $H_{{\varvec{k}}} = (-\i \nabla + {\varvec{k}})^2 + V$. The following result is concerned with the Bloch eigenmodes of the $H_{{\varvec{k}}}$’s.

Theorem 6

Let $B > 0$ and $V \in \mathcal {H}_{B,\mathbb {L}}$. For each ${\varvec{k}} \in \mathbb {R}^d$, the eigenfunctions of the Bloch fibers $H_{{\varvec{k}}} = (- i\nabla +{\varvec{k}})^2 + V$ of the periodic Schrödinger operator $H=-\Delta + V$ are in $\mathcal {H}_{A,\mathbb {L}}$ for any $0<A<B$. Let $\lambda _{1,{\varvec{k}}} \leqslant \lambda _{2,{\varvec{k}}} \leqslant \cdots$ be the eigenvalues of $H_{{\varvec{k}}}$ counted with multiplicities and ranked in non-decreasing order, and $\lambda _{1,{\varvec{k}},N} \leqslant \lambda _{2,{\varvec{k}},N} \leqslant \cdots \leqslant \lambda _{d_{\mathbb {L},N},{\varvec{k}},N}$ the eigenvalues of the variational approximation of $H_{{\varvec{k}}}$ in the $d_{\mathbb {L},N}$-dimensional space

$$\begin{aligned} X_{\mathbb {L},{\varvec{k}}, N}{:}{=}\,{\textrm{Span}}(e_{\varvec{G}}, \, {\varvec{G}} \in \mathbb {L}^*, \, |{\varvec{G}} + {\varvec{k}}| \leqslant N). \end{aligned}$$

Then, for each $0< A < B$ and $i \in \mathbb {N}^*$, there exists a constant $C_{i,A} \in \mathbb {R}_+$ such that

$$\begin{aligned} 0 \leqslant \max _{{\varvec{k}} \in \Omega ^*} (\lambda _{i,{\varvec{k}},N} - \lambda _{i,{\varvec{k}}}) \leqslant C_{i,A}{\textrm{e}}^{-2AN}, \end{aligned}$$

(19)

where $\Omega ^*$ is the first Brillouin zone (i.e. the Voronoi cell of the lattice $\mathbb {L}$ of $\mathbb {R}^d$ containing the origin).

Proof

It suffices to replace in the proofs of Theorems 2 and 3$X_N$ with $X_{\mathbb {L},{\varvec{k}}, N}$ and $T_{22}$ with the restriction $T_{22,\mathbb {L},{\varvec{k}}}$ of the operator $(- i\nabla + {\varvec{k}})^2$ to the invariant space $X_{\mathbb {L},{\varvec{k}}, N}^\perp$. The latter is invertible and such that $\Vert T_{22,\mathbb {L},{\varvec{k}}}^{-1} \Vert _{\mathcal {L}(X_{\mathbb {L},{\varvec{k}}, N}^\perp )} \leqslant N^{-2}$ and $\Vert T_{22,\mathbb {L},{\varvec{k}}}^{-1} \Vert _{\mathcal {L}(\mathcal {H}_{A,\mathbb {L}} \cap X_{\mathbb {L},{\varvec{k}}, N}^\perp )} \leqslant N^{-2}$. The result then follows similarly. $\square$

Remark 3

As a conclusion, let us mention that, under appropriate assumptions, we can extend these results to the KS-DFT equations (1) with GTH pseudopotentials in the case where an analytic parametrization of the exchange-correlation functional is used. This is relevant for instance when studying condensed-phase systems, for which the periodic setting is well suited, and where the commonly used approximations of $V_{\textrm{xc},\rho }$ are analytic on the positive density values that are of interest in such cases [19]. Then we can rewrite (1) as a system of elliptic PDEs:

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \left( - \frac{1}{2} \Delta + V_{\textrm{nl}} + V_{\textrm{loc}} + V_{\textrm{H},\rho } + V_{{\textrm{xc}},\rho } \right) \varphi _i - \lambda _i \varphi _i = 0, \quad (\varphi _i,\varphi _{j})_{L^2_{{\textrm{per}},\mathbb {L}}} = \delta _{ij}, \\ \displaystyle \rho ({\varvec{x}}) = \sum _{i=1}^{N_{\textrm{p}}} |\varphi _i({\varvec{x}})|^2, \\ \displaystyle - \Delta V_{\textrm{H},\rho }({\varvec{x}}) = 4\pi \left( \rho ({\varvec{x}}) - \frac{1}{|\Omega |}\int _\Omega \rho \right) , \quad \int _\Omega V_{\textrm{H},\rho } = 0, \end{array}\right. } \end{aligned}$$

We can then use known results for elliptic systems of PDEs [16, 17] to conclude, in a way similar to what we proved in Theorem 4, that the orbitals $\varphi _i$ belong to $\mathcal {H}_{A,\mathbb {L}}$ for some $A>0$. In particular, this leads to the exponential convergence of planewave approximations, justifying the use of GTH pseudopotentials.

B Proof of Lemma 1

We start by rewriting (18) as a first-order ODE on $\Psi _\varepsilon (y) {:}{=}\begin{bmatrix} \psi _\varepsilon (y) \\ \psi '_\varepsilon (y) \end{bmatrix} \in \mathbb {R}^2$, starting at $y_\eta \geqslant B_0$:

$$\begin{aligned} \Psi '_\varepsilon (y) = \begin{bmatrix} \Psi _{\varepsilon ,2}(y) \\ \varepsilon ^{-1} \left( \mu \sinh (y) - \Psi _{\varepsilon ,1}(y) + \Psi _{\varepsilon ,1}^3(y) \right) \end{bmatrix},\quad \Psi _\varepsilon (y_\eta ) {:}{=}\begin{bmatrix} \psi _\varepsilon (y_\eta ) \\ \psi '_\varepsilon (y_\eta ) \end{bmatrix}= \begin{bmatrix} 1+\eta \\ \psi '_\varepsilon (y_\eta ) \end{bmatrix}.\qquad \end{aligned}$$

(20)

We then use the following simplified version of more general comparison results on systems of differential inequalities [27, 28]. In the sequel, the inequality $a \geqslant b$ for two vectors $a,b \in \mathbb {R}^d$ means that $a_i \geqslant b_i$ for all $1 \leqslant \textrm{i} \leqslant d$.

Theorem 7

(See e.g. [27, p. 112]). Let $d\geqslant 1$ and $G: \mathbb {R}^d \rightarrow \mathbb {R}^d$ be locally Lipschitz and quasimonotone in the sense that for all $X,Z \in {\mathbb {R}}^d$,

$$\begin{aligned} \left( Z_i=X_i \hbox { and } X_{j} \leqslant Y_{j} \text{ for }\ j\ne i \right) \quad \Rightarrow \quad \left( G(X) \leqslant G(Z)\right) . \end{aligned}$$

Let $0 \leqslant y_0 < y_{\textrm{M}} \leqslant \infty$ and $\Phi \in C^1 ([y_0,y_{\textrm{M}}),\mathbb {R}^d)$ and $\Psi \in C^1 ([y_0,y_{\textrm{M}}),\mathbb {R}^d)$ satisfying respectively the ODE

$$\begin{aligned} \Phi '(y) = G(\Phi (y)),\qquad \Phi (y_0) \in \mathbb {R}^d, \end{aligned}$$

and the differential inequality

$$\begin{aligned} \Psi '(y) \geqslant G(\Psi (y)),\qquad \Psi (y_0) = \Phi (y_0). \end{aligned}$$

Then we have

$$\begin{aligned} \Psi (y) \geqslant \Phi (y). \end{aligned}$$

on $[y_0,y_\textrm{M})$.

To apply this result to (20), we introduce the function $G_\varepsilon : \mathbb {R}^2 \rightarrow \mathbb {R}^2$ defined by

$$\begin{aligned} \forall \ X = \begin{bmatrix} X_1 \\ X_2 \end{bmatrix} \in \mathbb {R}^2, \quad G_\varepsilon (X) = \begin{bmatrix} X_2 \\ \varepsilon ^{-1}(-X_1+X_1^3) \end{bmatrix}, \end{aligned}$$

and the maximal solution $\Phi _\varepsilon$ to

$$\begin{aligned} \Phi '_\varepsilon (y) = G_\varepsilon ( \Phi _\varepsilon (y)), \qquad \Phi _\varepsilon (y_\eta ) = \begin{bmatrix} 1+\eta \\ \psi _\varepsilon '(y_\eta ) \end{bmatrix}. \end{aligned}$$

(21)

As $\sinh (y) \geqslant 0$ for all $y \geqslant 0$, we have the following differential inequality:

$$\begin{aligned} \Psi '_\varepsilon (y) \geqslant G_\varepsilon (\Psi _\varepsilon (y)), \qquad \Psi _\varepsilon (y_\eta ) = \Phi _\varepsilon (y_\eta ). \end{aligned}$$

Note that, by convexity, $\psi _\varepsilon (y) \geqslant \psi _\varepsilon (y_\eta ) > 1$ for $y\geqslant y_\eta$: $G_\varepsilon$ is indeed quasimonotone on the domain of interest. We are now ready to compute an upper bound of $Y_\varepsilon$ with the use of Theorem 7, which yields

$$\begin{aligned} \forall \ y \geqslant y_\eta ,\ \Psi _\varepsilon (y) \geqslant \Phi _\varepsilon (y). \end{aligned}$$

This leads us to the study of the ODE

$$\begin{aligned} {\left\{ \begin{array}{ll} \varepsilon \phi ''_\varepsilon = -\phi _\varepsilon + \phi _\varepsilon ^3 \\ \phi _\varepsilon (y_\eta ) = 1+\eta ,\ \phi '_\varepsilon (y_\eta ) = \psi '_\varepsilon (y_\eta )>0. \end{array}\right. } \end{aligned}$$

We have

$$\begin{aligned} \frac{\varepsilon }{2}\frac{\textrm{d}}{\textrm{d}y} \left( \phi '_\varepsilon \right) ^2 = \frac{\textrm{d}}{\textrm{d}y}\left( -\frac{1}{2} \phi _\varepsilon ^2 + \frac{1}{4} \phi _\varepsilon ^4 \right) , \end{aligned}$$

from which we deduce

$$\begin{aligned} \frac{\varepsilon }{2}(\phi '_\varepsilon )^2 = \frac{1}{4}\phi _\varepsilon ^4 - \frac{1}{2}\phi _\varepsilon ^2 + C(\eta ), \end{aligned}$$

with

$$\begin{aligned}\begin{aligned} C(\eta )&= -\frac{1}{4}(1+\eta )^4 +\frac{1}{2}(1+\eta )^2 + \frac{\varepsilon }{2}(\psi '_\varepsilon (y_\eta ))^2\\&\geqslant -\frac{1}{4}(1+\eta )^4 +\frac{1}{2}(1+\eta )^2 + \frac{\varepsilon }{2}(\psi '_\varepsilon (y_0))^2,\\ \end{aligned} \end{aligned}$$

Fig. 7

Description of $X_\mu$ for $\mu =0.5$, along with the plot of $\psi _\varepsilon$ and the lower bound $\xi _{\varepsilon ,\eta }$ for $\varepsilon =0.1$, $\eta =0.5$. While $y < B_0$, $\psi _\varepsilon$ can possibly oscillate around $\psi _0$, but as soon as $y\geqslant B_0$, $\psi _\varepsilon$ is strictly convex and has no other choice than to explode in finite time $Y_\varepsilon \leqslant Y_{\varepsilon ,\eta }$, where $Y_{\varepsilon ,\eta }$ is the explosion time of the lower bound $\xi _{\varepsilon ,\eta }$

where $C(\eta )$ is computed from the initial conditions at $y=y_\eta$ and $B_0< y_0 < y_\eta$ is such that $\psi _\varepsilon (y_0)=1$. Note that $\eta \mapsto -\frac{1}{4}(1+\eta )^4 +\frac{1}{2}(1+\eta )^2 + \frac{\varepsilon }{2}(\psi '_\varepsilon (y_0))^2$ is decreasing on $\mathbb {R}_*^+$ and takes the value $\frac{1}{4} + \frac{\varepsilon }{2}(\psi '_\varepsilon (y_0))^2 > \frac{1}{4}$ at $\eta = 0$. Thus, for $\eta >0$ small enough, $C(\eta ) \geqslant \frac{1}{4}$ and we have

$$\begin{aligned} \begin{aligned} \frac{\varepsilon }{2}(\phi _\varepsilon ')^2&\geqslant \frac{1}{4}\phi _\varepsilon ^4 - \frac{1}{2}\phi _\varepsilon ^2 + \frac{1}{4} = \frac{1}{4}\left( \phi _\varepsilon ^2-1 \right) ^2, \end{aligned} \end{aligned}$$

hence

$$\begin{aligned} \phi _\varepsilon ' \geqslant \frac{1}{\sqrt{2\varepsilon }}\left( \phi _\varepsilon ^2-1 \right) \quad \text{ on } [y_\eta ,Y_\varepsilon ). \end{aligned}$$

Finally, we consider the ODE

$$\begin{aligned} \xi '_{\varepsilon ,\eta } = \frac{1}{\sqrt{2\varepsilon }}(\xi _{\varepsilon ,\eta }^2 -1),\qquad \xi _{\varepsilon ,\eta }(y_\eta ) = 1+\eta , \end{aligned}$$

whose solution is

$$\begin{aligned} \xi _{\varepsilon ,\eta }(y) = \frac{1+\frac{2}{\eta } + \exp \left( \frac{y-y_\eta }{\sqrt{\varepsilon /2}} \right) }{1+\frac{2}{\eta } - \exp \left( \frac{y-y_\eta }{\sqrt{\varepsilon /2}} \right) }, \end{aligned}$$

which is defined only up to $Y_{\varepsilon ,\eta } {:}{=}\sqrt{\frac{\varepsilon }{2}}\log \left( 1 + \frac{2}{\eta } \right) + y_\eta$. Applying again Theorem 7, we have that $\phi _\varepsilon (y) \geqslant \xi _{\varepsilon ,\eta }(y)$ for any $y \geqslant y_\eta$ such that both functions are still finite. Putting everything together, we obtain that $\psi _\varepsilon$ is only defined up to some $Y_\varepsilon$ with $B_0< Y_\varepsilon \leqslant Y_{\varepsilon ,\eta }<\infty$ and that

$$\begin{aligned} \forall \ y\in [y_\eta ,Y_\varepsilon ),\quad \psi _\varepsilon (y)\geqslant \xi _{\varepsilon ,\eta }(y). \end{aligned}$$

(22)

These results are illustrated on Fig. 7, where we plotted the lower bound $\xi _{\varepsilon ,\eta }$ for $\varepsilon =0.1$ and $\eta =0.5$ as well as a numerical approximation of $\psi _\varepsilon$.