# A priori error analysis of linear and nonlinear periodic Schrödinger equations with analytic potentials

Eric Cancès, Gaspard Kemlin, Antoine Levitt

<figure><img src="image_1.png" /></figure>

## To cite this version:

Eric Cancès, Gaspard Kemlin, Antoine Levitt. A priori error analysis of linear and nonlinear periodic Schrödinger equations with analytic potentials. Journal of Scientific Computing, 2023, 98 (1), pp.25. 10.1007/s10915-023-02421-0. hal-03692851v3

HAL Id: hal-03692851

https://inria.hal.science/hal-03692851v3

Submitted on 13 Nov 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A PRIORI ERROR ANALYSIS OF LINEAR AND NONLINEAR PERIODIC SCHRÖDINGER EQUATIONS WITH ANALYTIC POTENTIALS

ERIC CANCÈS<fnref n="1" />, GASPARD KEMLIN<fnref n="2" />, AND ANTOINE LEVITT<fnref n="3" />

**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.

## 1. INTRODUCTION

In this paper, we study models inspired by Kohn–Sham Density Functional Theory (KS-DFT). KS-DFT is currently the most popular model in quantum chemistry and materials science as it offers a good compromise between accuracy and computational efficiency. It aims at computing, for a given configuration of the nuclei of the molecular system or material of interest, the electronic ground-state energy and density. From the latter, it is possible to compute the effective forces acting on the nuclei in this configuration, and thus to identify the (meta)stable equilibrium configurations of the system, or to simulate the dynamics of the molecular system in various thermodynamic conditions. In materials science applications, computations are commonly done in a periodic simulation cell, which can be either the unit cell of a crystal (for the special case of perfect crystals), or a supercell (for all the other cases: crystals with defects, disordered alloys, glassy materials, liquids...).

We denote by $\mathbb{L} = \mathbb{Z}\mathbf{a}_1 + \mathbb{Z}\mathbf{a}_2 + \mathbb{Z}\mathbf{a}_3$ the periodic lattice, where $(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3)$ is a non-necessarily orthonormal basis of $\mathbb{R}^3$, and by $\Omega = [0, 1]\mathbf{a}_1 + [0, 1]\mathbf{a}_2 + [0, 1]\mathbf{a}_3$ the simulation cell. Let us denote by

$$L_{\text{per},\mathbb{L}}^2 := \{u \in L_{\text{loc}}^2(\mathbb{R}^3, \mathbb{C}) \mid u \text{ is } \mathbb{L}\text{-periodic}\}$$

the Hilbert space of complex-valued $\mathbb{L}$-periodic locally square integrable functions on $\mathbb{R}^3$, endowed with its usual inner product. The KS-DFT equations read

$$H_{\rho}\varphi_i = \lambda_i\varphi_i, \quad \langle\varphi_i, \varphi_j\rangle_{L_{\text{per},\mathbb{L}}^2} = \delta_{ij}, \quad \rho(\mathbf{x}) = 2\sum_{i=1}^{N_{\rho}}|\varphi_i(\mathbf{x})|^2, \qquad (1)$$

where $H_\rho$ is the Kohn–Sham Hamiltonian, a self-adjoint operator on $L_{\text{per},\mathbb{L}}^2$ bounded below and with compact resolvent. The $\varphi_i$'s are the Kohn–Sham orbitals, and the $\lambda_i$'s their energies. Since $H_\rho$ depends on $\rho$, which in turn depends on the eigenfunctions $\varphi_i$, $1 \le i \le N_\rho$, (1) is a nonlinear eigenproblem. The parameter $N_\rho$ represents physically the number of valence electron pairs per simulation cell and $\rho$ the ground-state electronic density. We assume here, as is the case for most physical systems, that $\lambda_1 \le \lambda_2 \le \dots \le \lambda_{N_\rho}$ are the lowest $N_\rho$ eigenvalues of $H_\rho$ (Aufbau principle). The Kohn–Sham Hamiltonian with pseudopotentials reads

$$H_{\rho} = -\frac{1}{2}\Delta + V_{\text{nl}} + V_{\text{loc},\rho}^{\text{Hxc}}$$

where $V_{\text{nl}}$ is a finite-rank self-adjoint operator (the nonlocal part of the pseudopotential), and

$$V_{\text{loc},\rho}^{\text{Hxc}}(\mathbf{x}) = V_{\text{loc}}(\mathbf{x}) + V_{\text{H},\rho}(\mathbf{x}) + V_{\text{xc},\rho}(\mathbf{x})$$

is a periodic real-valued function depending (nonlocally) on $\rho$. The function $V_{\text{loc}}$ is the local component of the pseudopotential, the Hartree potential $V_{\text{H},\rho}$ is the unique solution with zero mean to the periodic Poisson equation

$$-\Delta V_{\text{H},\rho}(\mathbf{x}) = 4\pi \left( \rho(\mathbf{x}) - \frac{1}{|\Omega|} \int_{\Omega} \rho \right), \quad \int_{\Omega} V_{\text{H},\rho} = 0,$$

<footnote n="1">¹CERMICS, ÉCOLE DES PONTS AND INRIA PARIS, 6 & 8 AVENUE BLAISE PASCAL, 77455 MARNE-LA-VALLÉE, FRANCE</footnote>

<footnote n="2">²CERMICS, ÉCOLE DES PONTS AND INRIA PARIS, 6 & 8 AVENUE BLAISE PASCAL, 77455 MARNE-LA-VALLÉE, FRANCE, NOW AT LAMFA, UNIVERSITÉ DE PICARDIE JULES VERNE UFR DES SCIENCES, 33 RUE SAINT-LEU, 80039 AMIENS, FRANCE</footnote>

<footnote n="3">³LMO, UNIVERSITÉ PARIS-SACLAY, 307 RUE MICHEL MAGAT, 91400 ORSAY, FRANCE</footnote>

<footnote>EMAILS: ERIC.CANCES@ENPC.FR, GASPARD.KEMLIN@U-PICARDIE.FR, ANTOINE.LEVITT@UNIVERSITE-PARIS-SACLAY.FR</footnote>
and the function $V_{xc,\rho}$, called the exchange-correlation potential, depends on the chosen approximation of the exchange-correlation energy functional. In the simple $X\alpha$ model [24], $V_{xc,\rho}(\mathbf{x}) = -C_D\rho(\mathbf{x})^{1/3}$, where $C_D > 0$ is the Dirac constant.

It is not mandatory to use pseudopotentials in KS-DFT calculations. Some softwares allow for all-electrons calculations in which the total pseudopotential operator $V_{nl}+V_{loc}$ is replaced with a local potential with Coulomb singularities at the positions of the nuclei. However, most calculations are done with pseudopotentials, or use the formally similar Projector Augmented Wave (PAW) method [5], for three reasons: (i) core electrons are barely affected by the chemical environment and can usually be considered to occupy “frozen states”, (ii) in heavy atoms, core electrons must be dealt with relativistic quantum models which makes the simulation more expensive from a computational viewpoint, (iii) due to the Coulomb singularities, all-electron Kohn–Sham orbitals have cusps at the positions of the nuclei and are therefore only Lipschitz continuous, while the Kohn–Sham orbitals computed with pseudopotentials are much more regular and can be well approximated with Fourier spectral methods (usually called planewave discretization methods in the field).

Several methods for constructing pseudopotentials have been proposed in the literature, leading to local and nonlocal functions of different regularity. As expected, the rate of convergence of the planewave discretization method is directly linked to the regularity of these functions. The a priori error analysis of this problem was performed in [7] for pseudopotential with Sobolev regularity. It was proved in particular, for the simple $X\alpha$ exchange-correlation functional, but also for the much more popular local density approximation (LDA) exchange-correlation functional, that if the local and nonlocal parts of the pseudopotential are in the periodic Sobolev space of order $s > 3/2$, then the Kohn–Sham orbitals $\phi_i$ and the density $\rho$ are in the periodic Sobolev space of order $s+2$, and (optimal) polynomial convergence rates were obtained in any Sobolev spaces of order $r$ with $-s < r < s+2$. In addition, as for linear second-order elliptic eigenproblems, the error on the eigenvalues converges to zero as the square of the error on the eigenfunctions evaluated in $H^1$-norm. The analysis in [7] covers for example the case of Troullier–Martins pseudopotentials [25], for which $s = \frac{7}{2} - \epsilon$ for all $\epsilon > 0$. On the other hand, these estimates are not sharp in the case of Goedecker–Teter–Hutter (GTH) pseudopotentials [10, 12], for which the local and nonlocal contributions are periodic sums of Gaussian-polynomial functions, and therefore have entire continuations to the whole complex plane. Such pseudopotentials are implemented in different DFT softwares, such as BigDFT [21], Quantum Espresso [9] or Abinit [11, 23], as well as DFTK, a recent electronic structure package in the Julia language [14].

The purpose of this paper is to investigate this case. While it has been known for a long time (see e.g. [3, 8, 16, 17, 20] and references therein for historical insight or [4, 13] for more recent developments) that the solutions to elliptic equations on $\mathbb{R}^d$ with real-analytic data have an analytic continuation in a complex neighborhood of $\mathbb{R}^d$, the size of this neighborhood is a priori unknown. In the periodic case we are considering, the latter directly impact the decay rate of the Fourier coefficients of the solution, hence the convergence rate of the planewave discretization method. For pedagogical reasons, we will work most of the time with one dimensional linear or nonlinear Schrödinger equations, because (i) it is easier to visualize analytic or entire continuations of functions originally defined on the real space $\mathbb{R}^d$ when $d = 1$, and (ii) exponential convergence rates of planewave discretization methods are easier to spot in 1D. However, most of our arguments extend to the multidimensional case. In Section 2, we introduce a hierarchy of spaces $(\mathcal{H}_A)_{A>0}$ of complex-valued $2\pi$-periodic functions on the real line having analytic continuations to the strip $\mathbb{R}+i(-A, A)$ of the complex plane. We then pick a real-valued function $V \in \mathcal{H}_B$ for some $B > 0$ and consider the one-dimensional Schrödinger operator $H = -\Delta + V$. A low vs high-frequency decomposition of the periodic $L^2$ space allows to prove that for all $0 < A < B$, the solution $u$ to the linear equation $Hu = f$ is in $\mathcal{H}_A$ whenever $f \in \mathcal{H}_A$ (see Section 3.1), and that the eigenfunctions of $H$ are in $\mathcal{H}_A$ (see Section 3.2). We rely on this result to prove in Section 3.3 that the planewave discretization method converges exponentially fast in this case and we provide a numerical illustration of these results. We turn in Section 4 to the nonlinear setting, where we present two examples for which we show, using a combination of numerical and analytical tools, that the analyticity strip of the solution may be much smaller than the one of the data. Finally, we consider in the Appendix the multidimensional case, which is an immediate extension, and its application to Kohn–Sham models.

## 2. SPACES OF ANALYTIC FUNCTIONS

Let us first introduce some notation. We denote by $L_{\text{per}}^2(\mathbb{R}, \mathbb{C})$ the space of locally square integrable complex-valued $2\pi$-periodic functions on $\mathbb{R}$, endowed with its natural inner product

$$ (u, v)_{L_{\text{per}}^2} := \int_0^{2\pi} \overline{u(x)} v(x) dx, $$
and by $S'_{\text{per}}(\mathbb{R}, \mathbb{C})$ the space of tempered complex-valued $2\pi$-periodic distributions on $\mathbb{R}$. For each $u \in S'_{\text{per}}(\mathbb{R}, \mathbb{C})$, we denote by $(\hat{u}_k)_{k \in \mathbb{Z}}$ the Fourier coefficients of $u$ with the following normalization convention:

$$ \forall u \in L_{\text{per}}^2(\mathbb{R}, \mathbb{C}), \quad \forall k \in \mathbb{Z}, \quad \hat{u}_k := (e_k, u)_{L_{\text{per}}^2} = \frac{1}{\sqrt{2\pi}} \int_0^{2\pi} u(x)e^{-ikx}dx, $$

where $e_k(x) := \frac{1}{\sqrt{2\pi}}e^{ikx}$ is the $L_{\text{per}}^2$-normalized Fourier mode with wavevector $k \in \mathbb{Z}$. Recall that the $2\pi$-periodic Sobolev spaces are the Hilbert spaces $H_{\text{per}}^s(\mathbb{R}, \mathbb{C})$, $s \in \mathbb{R}$, defined by

$$ H_{\text{per}}^s(\mathbb{R}, \mathbb{C}) := \left\{ u \in L_{\text{per}}^2(\mathbb{R}, \mathbb{C}) \;\middle|\; \sum_{k \in \mathbb{Z}} (1 + |k|^2)^s |\hat{u}_k|^2 < \infty \right\}, \\ (u, v)_{H_{\text{per}}^s} := \sum_{k \in \mathbb{Z}} (1 + |k|^2)^s \overline{\hat{u}_k} \hat{v}_k. $$

We will also use the self-explanatory notation $C_{\text{per}}^k(\mathbb{R}, \mathbb{R})$, $C_{\text{per}}^k(\mathbb{R}, \mathbb{C})$, $L_{\text{per}}^p(\mathbb{R}, \mathbb{R})$, $L_{\text{per}}^p(\mathbb{R}, \mathbb{C})$ for $k \in \mathbb{N} \cup \{\infty\}$ and $1 \le p \le \infty$, all these spaces being endowed with their natural norms or topologies. We now introduce, for any $A > 0$, the (Hardy-like) space

$$ \mathcal{H}_A := \left\{ u \in L_{\text{per}}^2(\mathbb{R}, \mathbb{C}) \;\middle|\; \sum_{k \in \mathbb{Z}} w_A(k) |\hat{u}_k|^2 < \infty \right\} $$

where

$$ w_A(k) := \cosh(2Ak) = \frac{1}{2}(e^{2Ak} + e^{-2Ak}), $$

endowed with the inner product

$$ (u, v)_A := \sum_{k \in \mathbb{Z}} w_A(k) \overline{\hat{u}_k} \hat{v}_k. $$

Note that $\mathcal{H}_A$ can be canonically identified with the space of analytic functions

$$ \tilde{\mathcal{H}}_A := \left\{ u : \Omega_A \to \mathbb{C} \text{ analytic} \;\middle|\; \begin{array}{l} \left[-A, A\right] \ni y \mapsto u(\cdot + iy) \in L_{\text{per}}^2(\mathbb{R}, \mathbb{C}) \text{ continuous,} \\ \displaystyle \int_0^{2\pi} \left( |u(x + iA)|^2 + |u(x - iA)|^2 \right) dx < \infty \end{array} \right\}, $$

where $\Omega_A := \mathbb{R} + i(-A, A) \subset \mathbb{C}$ is the horizontal strip of width $2A$ of the complex plane centered on the real axis, endowed with the inner product

$$ (u, v)_{\tilde{\mathcal{H}}_A} = \frac{1}{2} \left( (u(\cdot + iA), v(\cdot + iA))_{L_{\text{per}}^2} + (u(\cdot - iA), v(\cdot - iA))_{L_{\text{per}}^2} \right). $$

The canonical unitary mapping $\mathcal{H}_A$ onto $\tilde{\mathcal{H}}_A$ is the analytic continuation: any function $u \in \mathcal{H}_A$ has a unique analytic continuation $u : \Omega_A \to \mathbb{C}$ given by

$$ \forall z = x + iy \in \Omega_A, \quad u(z) = \sum_{k \in \mathbb{Z}} \hat{u}_k \frac{e^{ikz}}{\sqrt{2\pi}} = \sum_{k \in \mathbb{Z}} \hat{u}_k e^{-ky} e_k(x). $$

It can be easily seen that the Fourier coefficients of $u(\cdot \pm iA)$ are the Fourier coefficients of $u$ rescaled by a factor $e^{\mp kA}$ and that the function $(-A, A) \ni y \mapsto u(\cdot + iy) = \sum_{k \in \mathbb{Z}} \hat{u}_k e^{-ky} e_k(\cdot) \in L_{\text{per}}^2(\mathbb{R}, \mathbb{C})$ has a unique continuation to $[-A, A]$. Therefore,

$$ \begin{aligned} \|u\|_{\tilde{\mathcal{H}}_A}^2 &= \frac{1}{2} \left( \|u(\cdot + iA)\|_{L_{\text{per}}^2}^2 + \|u(\cdot - iA)\|_{L_{\text{per}}^2}^2 \right) \\ &= \frac{1}{2} \left( \sum_{k \in \mathbb{Z}} |\hat{u}_k e^{-kA}|^2 + \sum_{k \in \mathbb{Z}} |\hat{u}_k e^{+kA}|^2 \right) \\ &= \sum_{k \in \mathbb{Z}} w_A(k) |\hat{u}_k|^2 = \|u\|_A^2. \end{aligned} $$

We record for future use the

**Proposition 1.** Let $B > 0$. Then, for all $0 < A < B$, the multiplication by a function $V_B \in \mathcal{H}_B$ defines a bounded operator on $\mathcal{H}_A$.

*Proof.* The proof is immediate from the analyticity of $V$ in $\Omega_B$, which implies that $\sup_{z, |\text{Im}(z)| \le A} |V(z)| < \infty$. □
### 3. THE LINEAR CASE

**3.1. The linear elliptic problem.** We consider in a first stage the one-dimensional linear elliptic problem

$$ \text{seek } u \in H_{\text{per}}^2(\mathbb{R}, \mathbb{C}) \text{ such that } -\Delta u + Vu = f \text{ on } \mathbb{R}, \quad (2) $$

where $V \in L_{\text{per}}^2(\mathbb{R}, \mathbb{R})$ and $f \in L_{\text{per}}^2(\mathbb{R}, \mathbb{C})$ are given $2\pi$-periodic functions. We assume in this section that $V \ge 1$. It is then standard that (2) has a unique solution $u$ satisfying the a priori bounds

$$ \|u\|_{L_{\text{per}}^2} \le \frac{\|f\|_{L_{\text{per}}^2}}{\alpha} \quad \text{and} \quad \|u\|_{H_{\text{per}}^1} \le \|f\|_{H_{\text{per}}^{-1}}, \quad (3) $$

where $\alpha := \lambda_1(-\Delta + V) \ge 1$ is the smallest eigenvalue of the self-adjoint operator $H = -\Delta + V$ on $L_{\text{per}}^2(\mathbb{R}, \mathbb{C})$. By bootstrap arguments, $u \in H_{\text{per}}^{s+2}(\mathbb{R}, \mathbb{C})$ whenever $V$ and $f$ are in $H_{\text{per}}^s$, for any $s \ge 0$. The following result deals with the case of real-analytic potentials $V$ and right-hand sides $f$.

**Theorem 1.** Let $B > 0$ and $V \in \mathcal{H}_B$ be real-valued and such that $V \ge 1$ on $\mathbb{R}$. Then, for all $0 < A < B$ and $f \in \mathcal{H}_A$, the unique solution $u$ of (2) is in $\mathcal{H}_A$. Moreover, we have the following estimate

$$ \exists C > 0 \text{ independent of } f \text{ such that } \|u\|_A \le C \|f\|_A. \quad (4) $$

As a consequence, if $V$ and $f$ are entire, then so is $u$.

**Proof.** For $N > 0$, we consider the decomposition $L_{\text{per}}^2(\mathbb{R}, \mathbb{C}) = X_N \oplus X_N^{\perp}$ where

$$ X_N := \text{Span}(e_k, |k| \le N) = \{u \in L_{\text{per}}^2(\mathbb{R}, \mathbb{C}) \mid \hat{u}_k = 0, \forall |k| > N\}. \quad (5) $$

Let $\Pi_N$ be the orthogonal projector on $X_N$ and $\Pi_N^{\perp} := 1 - \Pi_N$ the orthogonal projector on $X_N^{\perp}$. Note that the restriction of $\Pi_N$ to the Sobolev space $H_{\text{per}}^s(\mathbb{R}, \mathbb{C})$, $s > 0$, is also the orthogonal projector on $X_N$ for the $H_{\text{per}}^s$ inner product, and that the same property holds for the Hilbert spaces $\mathcal{H}_A$.

For a fixed $N$, we decompose $u$ as $u = u_1 + u_2$ with $u_1 \in X_N$ and $u_2 \in X_N^{\perp}$. As $u_1$ has compact Fourier support, it obviously belongs to $\mathcal{H}_A$ and we have the estimate

$$ \|u_1\|_A \le \|u\|_{L_{\text{per}}^2} \sqrt{w_A(N)} \le \frac{\|f\|_{L_{\text{per}}^2}}{\alpha} \sqrt{w_A(N)}. \quad (6) $$

Let us show that, for $N$ large enough, $u_2$ also belongs to $\mathcal{H}_A$. Projecting $-\Delta u + Vu = f$ onto $X_N^{\perp}$, we get

$$ T_{22}u_2 + V_{22}u_2 = f_2 - V_{21}u_1, \quad (7) $$

where $T_{22}$ is the restriction to the invariant subspace $\mathcal{H}_A \cap X_N^{\perp}$ of the self-adjoint operator $-\Delta$ on $L_{\text{per}}^2(\mathbb{R}, \mathbb{C})$, and where, in view of Proposition 1, $V_{22} := \Pi_N^{\perp}V\Pi_N^{\perp} \in \mathcal{L}(\mathcal{H}_A \cap X_N^{\perp})$, $V_{21} := \Pi_N^{\perp}V\Pi_N \in \mathcal{L}(\mathcal{H}_A \cap X_N)$ and $f_2 := \Pi_N^{\perp}f \in \mathcal{H}_A \cap X_N^{\perp}$. The operator $T_{22}$ on $\mathcal{H}_A \cap X_N^{\perp}$ is bounded from below by $N^2$ and is therefore invertible with inverse $T_{22}^{-1}$ bounded by $N^{-2}$. As $f_2, V_{21}u_1 \in \mathcal{H}_A \cap X_N^{\perp}$, we can therefore rewrite (7) as

$$ (1 + T_{22}^{-1}V_{22})u_2 = T_{22}^{-1}(f_2 - V_{21}u_1). \quad (8) $$

Since $\|T_{22}^{-1}\|_{\mathcal{L}(\mathcal{H}_A \cap X_N^{\perp})} \le N^{-2}$ and $\|V_{22}\|_{\mathcal{L}(\mathcal{H}_A \cap X_N^{\perp})} \le \|V\|_{\mathcal{L}(\mathcal{H}_A)}$, we can choose $N$ large enough so that the operator $(1 + T_{22}^{-1}V_{22})$ is invertible on $\mathcal{H}_A \cap X_N^{\perp}$. It then holds

$$ u_2 = (1 + T_{22}^{-1}V_{22})^{-1}T_{22}^{-1}(f_2 - V_{21}u_1) \quad (9) $$

and the result follows from $\|V_{21}\|_{\mathcal{L}(\mathcal{H}_A)} \le \|V\|_{\mathcal{L}(\mathcal{H}_A)}$ and the bound (6) on $u_1$. $\square$

**Remark 1** (Additional regularity). We mention here some cases where we can obtain additional regularity on the solution $u$ to (2). First, note that, if we additionally require that $V(\cdot \pm iB) \in L_{\text{per}}^{\infty}(\mathbb{R}, \mathbb{C})$, then the same reasoning leads to $u \in \mathcal{H}_B$ if $f \in \mathcal{H}_B$. Next, one might wonder whether there exist $0 < B < B_0 < A$, $f \in \mathcal{H}_A$ and $V \in \mathcal{H}_B$ such that $V \notin \mathcal{H}_C$ for all $B_0 \le C \le A$ but nevertheless $u \in \mathcal{H}_A$. The following example shows that this can happen. Consider the $2\pi$-periodic real-valued functions $V$ and $u_0$ on $\mathbb{R}$ defined by

$$ \forall x \in \mathbb{R}, \quad V(x) = 1 + \frac{1}{\gamma + \sin^2(x)} \quad \text{and} \quad u_0(x) = \frac{\gamma + \sin^2(x)}{2\gamma + \sin^2(x)}. $$

We have that (i) $V \in \mathcal{H}_B$ if and only if $B < B_0 := \arcsin h(\sqrt{\gamma})$ and (ii) $u_0$ and $Vu_0$ both belong to $\mathcal{H}_A$ for any $A < A_0 := \arcsin h(\sqrt{2\gamma})$ (in particular for any $A$ such that $B_0 \le A < A_0$). Therefore, if we set $f := -\Delta u_0 + Vu_0$, then $f \in \mathcal{H}_A$ for any $A$ such that $B_0 < A < A_0$ and the unique solution to (2) with data $V$ and $f$ is $u = u_0 \in \mathcal{H}_A$ despite $V \notin \mathcal{H}_C$ for all $B_0 \le C \le A$.
## 3.2. The linear eigenvalue problem.
We now focus on the linear eigenvalue problem,

$$ -\Delta u + Vu = \lambda u, \quad \|u\|_{L_{\text{per}}^2(\mathbb{R}, \mathbb{C})} = 1, \qquad (10) $$

where $V \in \mathcal{H}_B$ for some $B > 0$. Using the same technique as for the proof of Theorem 1, we get the

**Theorem 2.** Let $B > 0$, $V \in \mathcal{H}_B$ be real-valued, and $(u, \lambda) \in H_{\text{per}}^2(\mathbb{R}, \mathbb{C}) \times \mathbb{R}$ be a normalized eigenmode of $H = -\Delta + V$, with isolated eigenvalue (i.e. a solution to (10)). Then, $u$ is in $\mathcal{H}_A$ for all $0 < A < B$. As a consequence, if $V$ is entire, then so is $u$.

*Proof.* The proof proceeds along the same lines as that of Theorem 1. Let $(u, \lambda) \in H_{\text{per}}^2(\mathbb{R}, \mathbb{C}) \times \mathbb{R}$ be a solution to (10). Using the same notation as in the proof of Theorem 1, we decompose $u$ as $u = u_1 + u_2$ with $u_1 \in X_N$ and $u_2 \in X_N^\perp$, and observe that for $N$ large enough,

$$ u_2 = -(1 + T_{22}^{-1}(V_{22} - \lambda))^{-1} T_{22}^{-1} V_{21} u_1, $$

with

$$ \|T_{22}^{-1}(V_{22} - \lambda)\|_{\mathcal{L}(\mathcal{H}_A \cap X_N^\perp)} \le \frac{\|V\|_{\mathcal{L}(\mathcal{H}_A)} + |\lambda|}{N^2}. $$

Therefore, choosing $N$ large enough, we have

$$ \|(1 + T_{22}^{-1}(V_{22} - \lambda))^{-1} T_{22}^{-1}\|_{\mathcal{L}(\mathcal{H}_A \cap X_N^\perp)} < 1, $$

from which we deduce that $u_2 \in \mathcal{H}_A$ and the result follows. $\square$

## 3.3. Planewave approximation of the linear Schrödinger equation.
Using $X_N = \text{Span}(e_k, |k| \le N) \subset H_{\text{per}}^1(\mathbb{R}, \mathbb{C})$ as a variational approximation space for (10), we obtain the finite-dimensional problem

$$ \left\{ \begin{array}{ll} \text{seek } (u_N, \lambda_N) \in X_N \times \mathbb{R} \text{ such that } \|u_N\|_{L_{\text{per}}^2(\mathbb{R}, \mathbb{C})} = 1 \text{ and} & \\ \forall v_N \in X_N, \quad \int_0^{2\pi} \overline{\nabla u_N} \cdot \nabla v_N + \int_0^{2\pi} V \overline{u_N} v_N = \lambda_N \int_0^{2\pi} \overline{u_N} v_N, & (11) \end{array} \right. $$

which is equivalent to seeking the eigenpairs of the Hermitian matrix $H_N$ with entries

$$ [H_N]_{kk'} := |k|^2 \delta_{kk'} + \hat{V}_{k-k'}, \quad k, k' \in \mathbb{Z}, \quad |k| \le N, \quad |k'| \le N. $$

The following theorem states that if $V \in \mathcal{H}_B$ for some $B > 0$, the planewave discretization method has an exponential convergence rate. Note that a similar result holds for the planewave approximation of the linear problem $-\Delta u + Vu = f$, whenever $f \in \mathcal{H}_A$.

**Theorem 3.** Let $B > 0$, $V \in \mathcal{H}_B$ be real-valued, $i \in \mathbb{N}^*$ and $0 < A < B$. Let $\lambda_i$ be the $i^{\text{th}}$ lowest eigenvalue of the self-adjoint operator $H = -\Delta + V$ on $L_{\text{per}}^2(\mathbb{R}, \mathbb{C})$ counting multiplicities, and $\mathcal{E}_i = \text{Ker}(H - \lambda_i)$ the corresponding eigenspace. For $N$ large enough, we denote by $\lambda_{i,N}$ the $i^{\text{th}}$ lowest eigenvalue of (11), and by $u_{i,N}$ an associated normalized eigenvector. Then, there is a constant $C_{i,A} \in \mathbb{R}_+$ such that, for $N$ large enough,

$$ d_{H_{\text{per}}^1}(u_{i,N}, \mathcal{E}_i) \le C_{i,A} \exp(-AN) \quad \text{and} \quad 0 \le \lambda_{i,N} - \lambda_i \le C_{i,A} \exp(-2AN). $$

*Proof.* First, note that $-\Delta + V$ has compact resolvent so that its eigenvalues $\lambda_i$ are isolated. Let $0 < A < B$ and $A' = \frac{A+B}{2}$. We have

$$ \forall v \in \mathcal{H}_{A'}, \quad \|v - \Pi_N v\|_{H_{\text{per}}^1} \le C_{A,B} \|v\|_{A'} e^{-AN} $$

with

$$ C_{A,B} := \left( \max_{k \in \mathbb{Z}} \frac{(1 + |k|^2) e^{2A|k|}}{w_{A'}(k)} \right)^{1/2} < \infty. $$

The operator $H$ is self-adjoint on $L_{\text{per}}^2(\mathbb{R}, \mathbb{C})$ with form domain $H_{\text{per}}^1(\mathbb{R}, \mathbb{C})$ and we infer from Theorem 2 that all the eigenfunctions of the operator $H$ are in $\mathcal{H}_{A'}$. Theorem 3 then follows from classical arguments on the variational approximations of the eigenmodes of bounded below self-adjoint operators with compact resolvent (see e.g. [1, Theorems 8.1 and 8.2]). $\square$

## 3.4. Numerical results in the linear case.

### 3.4.1. Computational framework.
In this paper, all the numerical tests are realized with the DFTK software [14]. This Julia package uses a planewave basis $X_N$, as defined in Section 3.3, through a discretization parameter $E_{\text{cut}} := N^2/2$. Then, the numerical strategy depends on the nature of the problem:

* linear eigenproblems are solved with a LOBPCG solver (see e.g. [18]);

* nonlinear eigenproblems (see Section 4.2.1) and nonlinear elliptic problems (see Section 4.2.2) with source terms are solved by direct minimization of the corresponding energy, through the Julia implementation of the LBFGS minimization algorithm [15].
### 3.4.2. Results in the linear case.
In this section, we provide some numerical experiments that illustrate Theorems 2 and 3. To this end, we consider the potential V defined by

$$ \forall x \in [0, 2\pi], \quad V(x) = \frac{1}{\gamma + \sin^2(x)}, $$

with $\gamma = 1/500$. The analytic continuation of V has branching points at $\pi\mathbb{Z} \pm iB_0$ with

$$ B_0 = \arcsinh(\sqrt{\gamma}) \approx 0.0447, $$

and a direct calculation shows that $V \in \mathcal{H}_B$ for any $B < B_0$. From the previous results, we therefore expect the eigenfunctions of $-\Delta + V$ to belong to $\mathcal{H}_B$ for any $B < B_0$ and the planewave approximation to converge with rate proportional to $\exp(-2B_0N)$ for the eigenvalues and proportional to $\exp(-B_0N)$ for the eigenfunctions in the $H^1$ norm.

In Figure 1, we compute numerically the ground state $u_1$ of $-\Delta + V$ with discretization parameter $E_{cut} = 5 \cdot 10^5$. In order to be able to properly spot the convergence of the Fourier coefficients, we use quadruple precision and the tolerance of the linear solver is set to $10^{-20}$. As expected, the Fourier coefficients of (the numerical approximation of) $u_1$ decrease with rate of order $\sqrt{w_{B_0}(k)}$, which is confirmed by looking at the convergence of two successive nonzero Fourier coefficients, since

$$ \lim_{k \to \infty} \log \left( \sqrt{\frac{\cosh(2B_0k)}{\cosh(2B_0(k+1))}} \right) = -B_0. $$

In Figure 2, we solved the same problem for various $N$'s, and computed the error on the eigenvalues and the eigenvectors in the $H^1$ norm with respect to a reference solution obtained with $N = \sqrt{2E_{cut}} = 1000$. This time, the tolerance of the linear solver is set to $10^{-13}$ and we use standard double precision. The convergence rates predicted by Theorem 3 are observed, as expected. We can also remark that the error on the eigenvalues seems to decrease slightly faster (with an additional $1/N^2$ factor), the asymptotic rate being still given by the exponential factor $\exp(-2B_0N)$.

<figure><img src="image_3.png" /><figcaption>FIGURE 1 – [Linear eigenvalue problem] (Left) Fourier coefficients of (the numerical approximation of) $u_1$: they seem to decrease like $1/(|k|^\alpha \sqrt{\cosh(2B_0k)})$ with $\alpha = 4$, the prefactor $2B_0$ in the cosh being confirmed by the plot on the right. (Right) Logarithm of the ratio of two successive nonzero Fourier coefficients.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 2 – [Linear eigenvalue problem] (Left) Discretization error for the first eigenvalue and the associated eigenvector in the $H^1$ norm. Both errors seem to converge with rates given by Theorem 3, the prefactors in the cosh being confirmed by the plot on the right. (Right) Logarithm of the ratio between discretization errors for various N's.</figcaption></figure>

## 4. THE NONLINEAR CASE

In the nonlinear case, the analyticity strip of the solution can be much smaller than the one of the data. First, we report numerical simulations on a nonlinear periodic Gross-Pitaevskii eigenvalue problem with a cubic nonlinearity and a potential admitting an entire continuation, illustrating the fact that the eigenfunction belongs to $H_A$ for some $A > 0$, but does not seem to be entire. Next, we study a nonlinear Schrödinger equation with entire potential and source term and cubic nonlinearity. We show that the solution is not entire and we provide an upper bound of the width of the horizontal analyticity strip of the solution.

### 4.1. Analyticity on a strip.

$$ - \Delta u + V u + |u|^2 u = \lambda u, \quad \|u\|_{L_{\text{per}}^2(\mathbb{R}, \mathbb{C})} = 1, \qquad (12) $$

where $V \in H_B$ for any $B > 0$ is real-valued. Similarly to the linear case, we consider the planewave variational approximation: using again

$$ X_N = \operatorname{Span}(e_k, \quad |k| \le N), $$

we obtain the finite-dimensional problem

$$ \left\{ \begin{array}{l} \text{seek } (u_N, \lambda_N) \in X_N \times \mathbb{R} \text{ such that } \|u_N\|_{L_{\text{per}}^2(\mathbb{R}, \mathbb{C})} = 1 \text{ and} \\ \forall v_N \in X_N, \quad \int_0^{2\pi} \overline{\nabla u_N} \cdot \nabla v_N + \int_0^{2\pi} V \overline{u_N} v_N + \int_0^{2\pi} |u_N|^2 \overline{u_N} v_N = \lambda_N \int_0^{2\pi} \overline{u_N} v_N. \end{array} \right. \qquad (13) $$

Then, the following theorem holds, which is a simple corollary of known results on the analyticity of elliptic partial differential equations. Note that this result only yields the existence of a finite strip of analyticity in the complex plane. However, even if the data is entire, there is a priori no reason for the solution to be entire too, as shown by the counter-examples that follow.

**Theorem 4.** Let $V$ real-valued be in $H_B$ for any $B > 0$ and $(u, \lambda) \in H_{\text{per}}^2(\mathbb{R}, \mathbb{R}) \times \mathbb{R}$ be the ground-state of (12), uniquely defined under the assumption that $u > 0$. Then, there exists $A > 0$ such that $u \in H_A$. Moreover, if $(u_N, \lambda_N)$ is the variational approximation of $(u, \lambda)$ in $X_N \times \mathbb{R}$, then

$$ \exists C_A > 0, \quad \|u - u_N\|_{H_{\text{per}}^1} \le C_A \exp(-AN). $$

*Proof.* Existence and uniqueness of $u$ is a classical result, proved for instance in [6, Appendix]. In addition, as $V \in H_{\text{per}}^s(\mathbb{R}, \mathbb{R})$ for any $s > 0$, we know from [6, Theorem 2] and by an immediate bootstrap argument that $u$ also belongs to $H_{\text{per}}^s(\mathbb{R}, \mathbb{R})$ for any $s > 0$. Therefore, $u \in C_{\text{per}}^{\infty}(\mathbb{R}, \mathbb{R})$.
Knowing that $u \in C_{\text{per}}^{\infty}(\mathbb{R}, \mathbb{R})$ yields that $u$ is real-analytic on a neighborhood of every point $x_0$ of the periodicity cell $[0, 2\pi]$, according to well-known results on the analyticity of solutions to nonlinear elliptic equations, which we recall here for the reader's convenience (see for instance [4, Theorem 1.1], [13, Theorem 1] or [8, 16, 17]).

**Theorem 5.** Let $\Omega$ be an open subset of $\mathbb{R}^n$, $D$ an open subset of $\mathbb{R}^n \times \mathbb{R} \times \mathbb{R}^n \times \mathbb{R}_{\text{sym}}^{n \times n}$, $F: D \ni (x, v, g, M) \mapsto F(x, v, g, M) \in \mathbb{R}$ a real analytic function, and $u \in C^{\infty}(\Omega)$ a real-valued function such that for all $x \in \Omega$, $(x, u(x), \nabla u(x), \nabla^2 u(x)) \in D$. Assume that $u$ solves $F(x, u(x), \nabla u(x), \nabla^2 u(x)) = 0$ and that this equation is elliptic in the sense that

$$ \forall (x, \xi) \in \Omega \times \mathbb{R}^n, \xi \neq 0, \quad \xi^T \nabla_M F(x, u(x), \nabla u(x), \nabla^2 u(x)) \xi \neq 0. $$

Then $u$ is real-analytic in $\Omega$.

By a compactness argument, we therefore have that $u$ is analytic on a strip of size $A > 0$ around the real axis, for some $A > 0$.

We now consider the variational approximation $(u_N, \lambda_N)$ of $(u, \lambda)$ in $X_N \times \mathbb{R}$. We know from [6, Theorem 1] that there exists $C > 0$ such that

$$ \|u - u_N\|_{H_{\text{per}}^1} \le C \|u - \Pi_N u\|_{H_{\text{per}}^1}. $$

The result then follows similarly to the proof of Theorem 3. $\square$

**Remark 2.** In order to establish the convergence rate of the eigenvalues of the ground-state of (12), one can follow the proof of [6, Theorem 2]. Using in particular equations (50) and (53) from this reference, with the negative Sobolev norms replaced by the dual norms of $H_A$, one can obtain that there exists $C_A > 0$ such that $|\lambda - \lambda_N| \le C_A \exp(-2AN)$.

### 4.2. Counter-examples.
We analyze in this section two counter-examples of Theorems 1 and 2 based on the nonlinear Gross-Pitaevskii equation with entire potentials and source terms, for which the solutions are not entire. We also present numerical tests to support our analysis.

#### 4.2.1. A nonlinear Gross-Pitaevskii eigenvalue problem.
We study in this section the following nonlinear eigenvalue problem:

$$ \begin{cases} -\varepsilon \Delta u_{\varepsilon} + \gamma \cos(x) u_{\varepsilon} + |u_{\varepsilon}|^2 u_{\varepsilon} = \lambda_{\varepsilon} u_{\varepsilon} & \text{in } H_{\text{per}}^{1}(\mathbb{R}, \mathbb{C}), \\ \|u_{\varepsilon}\|_{L_{\text{per}}^2} = 1, & u_{\varepsilon} > 0, \end{cases} \quad (14) $$

with entire data. In view of the results we obtained in the linear case, we could expect the solution $u_{\varepsilon}$ to (14) (which we know to exist) to be entire. However, our numerical experiments suggest on the contrary that $u_{\varepsilon}$ is analytic only on a band of finite size.

The potential $V$ being real-valued, $u_{\varepsilon}$ is actually real-valued too and the singular limit $\varepsilon \to 0$ yields the algebraic equation $\gamma \cos(x)u_0(x) + u_0(x)^3 = \lambda_0 u_0(x)$. If we assume that $u_0$ does not vanish on the unit cell $[0, 2\pi]$, we can divide by $u_0(x)$ and then integrate over $[0, 2\pi]$: using the normalization condition, we obtain $\lambda_0 = 1/2\pi$. This yields

$$ u_0(x) = \sqrt{\frac{1}{2\pi} - \gamma \cos(x)}, \quad (15) $$

which is indeed bounded away from zero as soon as $\beta_0 := 1/(2\pi\gamma) > 1$, which we take in the following. The analytic continuation of $u_0$, still denoted by $u_0$, satisfies

$$ \gamma \cos(z)u_0(z) + u_0(z)^3 = \frac{1}{2\pi}u_0(z), $$

with $\sqrt{\cdot}$ in (15) being the continuation of the square root with branch cut $\mathbb{R}_-$. The maximal horizontal strip of the complex plane on which $u_0$ is analytic is $\mathbb{R} + i(-B_0, B_0)$, where $B_0$ is such that $\cos(\pm iB_0) = \beta_0$, that is to say $B_0 = \log(\beta_0 + \sqrt{\beta_0^2 - 1})$. More precisely, the function $z \mapsto u_0(z)$ admits two branching points at $z_{\pm} = \pm iB_0$, for which $u_0(z_{\pm}) = 0$ and

$$ \left| \frac{du_0}{dt}(tz_{\pm}) \right| = \left| \frac{\gamma \sin(tz_{\pm})}{2\sqrt{\frac{1}{2\pi} - \gamma \cos(tz_{\pm})}} \right| \to \infty \quad \text{when } \mathbb{R} \ni t \to 1_{-}, $$

see Figure 3.

When $\varepsilon > 0$, we can approximate numerically the solution of (14) using the planewave approximation presented before, with discretization parameter $E_{\text{cut}} = 10^6$ and solver tolerance $10^{-13}$. The Fourier coefficients of the numerical solution for various $\varepsilon$ are displayed in Figure 4. Although it is not possible to rule out a transition to a different asymptotic regime in the limit of extremely small values of $\varepsilon$ or very large values of
<figure><img src="image_5.png" /><figcaption>FIGURE 3 – [Nonlinear eigenvalue problem] Analytic continuation of $u_0$ for $\beta_0 = 1.00001$, for which $B_0 \approx 0.0045$. (Left) Imaginary and real parts of $y \mapsto u_0(iy)$. A discontinuity of the derivative appears at the expected positions. (Right) Phase of $z \mapsto u_0(z)$ where $z = x + iy$. Branching points, in red, appear at the expected positions.</figcaption></figure>

$k$ because of finite numerical accuracy, our numerical results are compatible with the decay of the Fourier coefficients of $u_\epsilon$ as $1/(|k|^{3/2}\sqrt{\cosh(2B_0k)})$. This leads to the conclusion that the solution $u_\epsilon$ is analytic on a band strip of size close to $B_0$, which is in contradiction with the results from the linear case. In the next section, we study a similar nonlinear elliptic equation where we reach similar conclusions and where, in addition, we can provide an upper bound of the width of the horizontal strip on which the solutions are analytic.

<figure><img src="image_6.png" /><figcaption>FIGURE 4 – [Nonlinear eigenvalue problem] (Left) Fourier coefficients of (the numerical approximation of) $u_\epsilon$: they seem to decrease like $1/(|k|^\alpha \sqrt{\cosh(2B_0k)})$ with $\alpha = 3/2$, the prefactor $2B_0$ in the cosh being confirmed by the plot on the right. (Right) Logarithm of the ratio of two successive nonzero Fourier coefficients.</figcaption></figure>
### 4.2.2. A nonlinear elliptic equation.
Still in the perspective of studying nonlinear elliptic problems with analytic data, we now consider the nonlinear periodic elliptic equation with cubic nonlinearity

$$ -\varepsilon\Delta u_{\varepsilon} + u_{\varepsilon} + |u_{\varepsilon}|^2 u_{\varepsilon} = f \quad \text{in } H_{\text{per}}^{1}(\mathbb{R}, \mathbb{C}), \qquad (16) $$

where $\varepsilon > 0$, and $f: \mathbb{R} \to \mathbb{R}$ is a real-analytic $2\pi$-periodic function admitting an entire continuation, still denoted by $f$, to the complex plane. We will show that, in this particular case, the same kind of results we showed for the linear case are not true any more and we provide an estimation of the width of the horizontal analyticity strip of the solution, which is finite even though the source term $f$ is entire.

The right-hand-side $f$ being real-valued, $u_{\varepsilon}$ is also real-valued and the singular limit $\varepsilon = 0$ gives rise to the algebraic equation $u_0(x) + u_0(x)^3 = f(x)$, which has a unique real solution for each $x \in \mathbb{R}$. The latter can be computed by Cardano's formula: the discriminant of the cubic equation is

$$ R(x) = -(4 + 27f^2(x)) < 0, $$

so that

$$ u_0(x) = \sqrt[3]{\frac{1}{2}\left(f(x) + \sqrt{\frac{-R(x)}{27}}\right)} + \sqrt[3]{\frac{1}{2}\left(f(x) - \sqrt{\frac{-R(x)}{27}}\right)}, \qquad (17) $$

the other two roots being complex conjugates with nonzero imaginary parts. The function $u_0$ can be analytically continued from the real axis upwards in the complex plane as long as $R(z)$ does not touch zero; it has branching points (with an exploding first derivative) at the points where $R(z) = 0$. In the rest of this section, we will use the function $f(x) = \mu \sin(x)$, for a fixed $\mu > 0$. The analytic continuation of $u_0$ has branching points at $2\pi\mathbb{Z} \pm iB_0$, with

$$ B_0 = \operatorname{arcsinh}\sqrt{\frac{4}{27\mu^2}} > 0. $$

In particular, although the source term $f$ has an entire continuation, the solution $u_0$ given by (17) does not.

When $\varepsilon > 0$, we can approximate numerically the solution to (16) with the planewave approximation introduced before, with discretization parameter $E_{\text{cut}} = 10^6$ and solver tolerance $10^{-13}$. The Fourier coefficients of the numerical solution are then displayed in Figure 6 for various $\varepsilon$: the plots suggest that increasing $\varepsilon$ increases the width $B_{\varepsilon}$ of the horizontal analyticity strip of $u_{\varepsilon}(z)$, but does not make it entire.

<figure><img src="image_7.png" /><figcaption>FIGURE 5 – [Nonlinear elliptic equation] Analytic continuation of $u_0$ for $\mu = 10$, for which $B_0 \approx 0.0385$. (Left) Imaginary part of $y \mapsto u_0(iy)$. Discontinuities also appear at the expected positions $\pm B_0$. (Right) Phase of $z \mapsto u_0(z)$ where $z = x + iy$. Branching points, in red, appear at the expected positions $\pm iB_0$.</figcaption></figure>

In this particular case, we are able to obtain an upper bound of the value of $B_{\varepsilon}$ by an ODE technique. Let $\varphi_{\varepsilon}(y) := u_{\varepsilon}(iy)$ be the value of the solution along the imaginary axis. Using the Cauchy–Riemann equations,
<figure><img src="image_8.png" /><figcaption>FIGURE 6 – [Nonlinear elliptic equation] (Left) Fourier coefficients of (the numerical approximation of) $u_\epsilon$. (Right) Logarithm of the ratio of two successive nonzero Fourier coefficients.</figcaption></figure>

we see that $\varphi_\epsilon$ satisfies the following second-order ODE:

$$ \begin{cases} \varepsilon\varphi_{\varepsilon}''(y) + \varphi_{\varepsilon}(y) + \varphi_{\varepsilon}^3(y) = i\mu\sinh(y), \\ \varphi_{\varepsilon}(0) = u_{\varepsilon}(0) = 0, \quad \varphi_{\varepsilon}'(0) = iu_{\varepsilon}'(0), \end{cases} $$

with $u'_\epsilon(0) \in \mathbb{R}$. Decomposing $\varphi_\epsilon$ in its real part $\theta_\epsilon$ and imaginary part $\psi_\epsilon$, we see that $\theta_\epsilon$ satisfies the equation

$$ \begin{cases} \varepsilon\theta_{\varepsilon}'' + \theta_{\varepsilon} + \theta_{\varepsilon}^3 - 3\theta_{\varepsilon}\psi_{\varepsilon}^2 = 0, \\ \theta_{\varepsilon}(0) = 0, \quad \theta_{\varepsilon}'(0) = 0, \end{cases} $$

and therefore vanishes. As a consequence, $u_\epsilon$ is purely imaginary along the imaginary axis (as could have been anticipated by the symmetries $u_\epsilon(-z) = -u(z)$, $u_\epsilon(\bar{z}) = \overline{u_\epsilon(z)}$). We are then left with studying the imaginary part $\psi_\epsilon$, which satisfies the ODE

$$ \begin{cases} \varepsilon\psi_{\varepsilon}'' + \psi_{\varepsilon} - \psi_{\varepsilon}^3 = \mu\sinh, \\ \psi_{\varepsilon}(0) = 0, \quad \psi_{\varepsilon}'(0) = u_{\varepsilon}'(0). \end{cases} \quad (18) $$

If we can prove that $\psi_\epsilon$ becomes non-analytic at a finite $0 < Y_\epsilon < \infty$, this will imply that the width $B_\epsilon$ of the horizontal analyticity strip of $u_\epsilon$ satisfies $B_\epsilon \le Y_\epsilon < \infty$ and is therefore finite even though the data in (16) is entire. To this end, we prove in the Appendix B the following lemma, which gives a sufficient condition for the non-analyticity of $\psi_\epsilon$ in finite time.

**Lemma 1.** Let $Y_\epsilon \in \mathbb{R}_+ \cup \{\infty\}$ be the maximum time of definition of $\psi_\epsilon$. If, for some $\eta > 0$ small enough, there exists $y_\eta \ge B_0$ such that $\psi_\epsilon(y_\eta) = 1 + \eta > 1$ and $\psi'_\epsilon(y_\eta) > 0$, then $Y_\epsilon < \infty$ and $\psi_\epsilon(y) \to \infty$ when $y \to Y_\epsilon$.

To prove that the sufficient condition from Lemma 1 is satisfied, we introduce the set

$$ X_{\mu} = \{(y, v) \in \mathbb{R}^2, \mu \sinh(y) - v + v^3 \ge 0\}. $$

This set is such that, if $(y, \psi_\epsilon(y))$ lies strictly inside $X_\mu$, then $\psi_\epsilon$ is locally strictly convex. For $y < B_0$, $\psi_\epsilon$ might oscillate on both sides of the boundary of $X_\mu$ [26]. We investigated numerically the behavior of this function for the set of parameters used in Figure 7 ($\epsilon = 0.1$, $\mu = 0.5$), and observed that $0 < \psi_\epsilon(B_0) < \frac{1}{\sqrt{3}}$ and $\psi'_\epsilon(B_0) > 0$. This numerical observation can be trusted as the ODE satisfied by $\psi'_\epsilon$ on the interval $[0, B_0]$ for $\epsilon = 0.1$ and $\mu = 0.5$ is not stiff. It is therefore easy to solve it numerically with high accuracy with a posteriori error estimates guaranteeing that $\psi_\epsilon(B_0)$ is indeed strictly between 0 and $\frac{1}{\sqrt{3}}$, and $\psi'_\epsilon(B_0)$ is positive. Therefore, $\psi_\epsilon$ is locally strictly convex at $y = B_0$. Given the shape of $X_\mu$ (see Figure 7 in appendix), $\psi_\epsilon$ lies, for any $B_0 < y < Y_\epsilon$, above its tangent at $y = B_0$, whose slope is $\psi'_\epsilon(B_0) > 0$. Strict convexity allows one to
conclude that for any $\eta > 0$, there exists $y_\eta > B_0$ such that

$$ \begin{cases} \psi_{\varepsilon}(y_{\eta}) = 1 + \eta > 1, \\ \psi'_{\varepsilon}(y_{\eta}) > \psi'_{\varepsilon}(B_0) > 0. \end{cases} $$

Any $\eta > 0$ is therefore suitable to apply Lemma 1 and conclude that $\psi_{\varepsilon}$ blows up in finite time $Y_{\varepsilon}$. We can deduce from these investigations that $u_{\varepsilon}$ is analytic only on a horizontal strip of finite width of the complex plane although the source term $f$ is an entire function: our results in the linear case are therefore no longer valid in general in the nonlinear case.

## CONFLICT OF INTEREST

All the authors declare that they have no conflicts of interest.

## DATA AVAILABILITY

All the scripts used to generate the plots of Figures 1, 2, 3, 4, 5 and 6 are available online at

https://github.com/gkemlin/analytic potenti als.

## ACKNOWLEDGEMENTS

The authors would like to thank Geneviève Dusson and Michael F. Herbst for fruitful discussions. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 810367). The authors would also like to thank the anonymous reviewers for their comments and suggestions.

## 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 $L = Za_1 + \dots + Za_d$ where $a_1, \dots, 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 $L = 2\pi\mathbb{Z}^d$. We denote by $\Omega = [0, 2\pi)^d$ a unit cell, by $L^* = \mathbb{Z}^d$ the reciprocal lattice, by $e_G(x) = |\Omega|^{-1/2}e^{iG \cdot x}$ the Fourier mode with wavevector $G \in L^*$, and by

$$ H_{\text{per},L}^s := \left\{ u = \sum_{G \in L^*} \hat{u}_G e_G \left| \sum_{G \in L^*} (1 + |G|^2)^s |\hat{u}_G|^2 < \infty \right. \right\} $$

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

$$ \mathcal{H}_{A,L} := \left\{ u \in L_{\text{per},L}^2 \left| \sum_{G \in L} w_{A,L}(G) |\hat{u}_G|^2 < \infty \right. \right\}, \quad (u,v)_{A,L} := \sum_{G \in L^*} w_{A,L}(G) \overline{\hat{u}_G} \hat{v}_G, $$

where $w_{A,L}(G) = \sum_{n=1}^d w_A(G_n)$. Each $u \in \mathcal{H}_{A,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

$$ \sum_{G \in L^*} w_{A,L}(G) |\hat{u}_G|^2 = \frac{1}{2} \sum_{j=1}^d \int_{\Omega} |u(x + ie_j)|^2 + |u(x - ie_j)|^2 dx. $$

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

$$ X_{N,L} := \text{Span}(e_G, G \in L^*, |G| \le N), $$

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

$$ \|T_{22,L}^{-1}\|_{\mathcal{L}(X_{N,L}^{\perp})} = \|T_{22,L}^{-1}\|_{\mathcal{L}(\mathcal{H}_{A,L} \cap X_{N,L}^{\perp})} \le N^{-2}. $$

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

Lastly, if $V \in \mathcal{H}_{B,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 $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_{\text{per},L}^2$ with domain $H_{\text{per},L}^2$ and form domain $H_{\text{per},L}^1$ defined as $H_k = (-i\nabla + k)^2 + V$. The following result is concerned with the Bloch eigenmodes of the $H_k$'s.
**Theorem 6.** Let $B > 0$ and $V \in \mathcal{H}_{B,\mathbb{L}}$. For each $k \in \mathbb{R}^d$, the eigenfunctions of the Bloch fibers $H_k = (-i\nabla + 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,k} \le \lambda_{2,k} \le \dots$ be the eigenvalues of $H_k$ counted with multiplicities and ranked in non-decreasing order, and $\lambda_{1,k,N} \le \lambda_{2,k,N} \le \dots \le \lambda_{d_{\mathbb{L},N},k,N}$ the eigenvalues of the variational approximation of $H_k$ in the $d_{\mathbb{L},N}$-dimensional space

$$X_{\mathbb{L},k,N} := \operatorname{Span}(e_G, G \in \mathbb{L}^*, |G + k| \le N).$$

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

$$0 \le \max_{k \in \Omega^*} (\lambda_{i,k,N} - \lambda_{i,k}) \le C_{i,A}e^{-2AN}, \qquad (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 Theorem 2 and Theorem 3 $X_N$ with $X_{\mathbb{L},k,N}$ and $T_{22}$ with the restriction $T_{22,\mathbb{L},k}$ of the operator $(-i\nabla + k)^2$ to the invariant space $X_{\mathbb{L},k,N}^{\perp}$. The latter is invertible and such that $\|T_{22,\mathbb{L},k}^{-1}\|_{\mathcal{L}(X_{\mathbb{L},k,N}^{\perp})} \le N^{-2}$ and $\|T_{22,\mathbb{L},k}^{-1}\|_{\mathcal{L}(\mathcal{H}_{A,\mathbb{L}} \cap X_{\mathbb{L},k,N}^{\perp})} \le 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_{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:

$$\left\{
\begin{array}{l}
\left(-\frac{1}{2}\Delta + V_{nl} + V_{loc} + V_{H,\rho} + V_{xc,\rho}\right) \varphi_i - \lambda_i \varphi_i = 0, \quad (\varphi_i, \varphi_j)_{L_{\text{per},\mathbb{L}}^2} = \delta_{ij}, \\[1.5ex]
\rho(\mathbf{x}) = \sum_{i=1}^{N_p} |\varphi_i(\mathbf{x})|^2, \\[1.5ex]
-\Delta V_{H,\rho}(\mathbf{x}) = 4\pi \left(\rho(\mathbf{x}) - \frac{1}{|\Omega|} \int_{\Omega} \rho\right), \quad \int_{\Omega} V_{H,\rho} = 0,
\end{array}
\right.
$$

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.

## APPENDIX 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} \ge B_0$:

$$\Psi'_{\varepsilon}(y) = \begin{bmatrix} \Psi_{\varepsilon,2}(y) \\ \varepsilon^{-1}(\mu \sinh(y) - \Psi_{\varepsilon,1}(y) + \Psi_{\varepsilon,1}^3(y)) \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 (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 \ge b$ for two vectors $a, b \in \mathbb{R}^d$ means that $a_i \ge b_i$ for all $1 \le i \le d$.

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

$$(Z_i = X_i \text{ and } X_j \le Y_j \text{ for } j \ne i) \Rightarrow (G(X) \le G(Z)).$$

Let $0 \le y_0 < y_M \le \infty$ and $\Phi \in C^1([y_0, y_M], \mathbb{R}^d)$ and $\Psi \in C^1([y_0, y_M], \mathbb{R}^d)$ satisfying respectively the ODE

$$\Phi'(y) = G(\Phi(y)), \quad \Phi(y_0) \in \mathbb{R}^d,$$

and the differential inequality

$$\Psi'(y) \ge G(\Psi(y)), \quad \Psi(y_0) = \Phi(y_0).$$

Then we have

$$\Psi(y) \ge \Phi(y).$$

on $[y_0, y_M)$.

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

$$\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},$$
and the maximal solution $\Phi_{\varepsilon}$ to

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

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

$$ \Psi'_{\varepsilon}(y) \ge G_{\varepsilon}(\Psi_{\varepsilon}(y)), \quad \Psi_{\varepsilon}(y_{\eta}) = \Phi_{\varepsilon}(y_{\eta}). $$

Note that, by convexity, $\psi_{\varepsilon}(y) \ge \psi_{\varepsilon}(y_{\eta}) > 1$ for $y \ge 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

$$ \forall y \ge y_{\eta}, \Psi_{\varepsilon}(y) \ge \Phi_{\varepsilon}(y). $$

This leads us to the study of the ODE

$$ \begin{cases} \varepsilon\phi''_{\varepsilon} = -\phi_{\varepsilon} + \phi_{\varepsilon}^3 \\ \phi_{\varepsilon}(y_{\eta}) = 1 + \eta, \quad \phi'_{\varepsilon}(y_{\eta}) = \psi'_{\varepsilon}(y_{\eta}) > 0. \end{cases} $$

We have

$$ \frac{\varepsilon}{2} \frac{d}{dy} (\phi'_{\varepsilon})^2 = \frac{d}{dy} \left( -\frac{1}{2}\phi_{\varepsilon}^2 + \frac{1}{4}\phi_{\varepsilon}^4 \right), $$

from which we deduce

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

with

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

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) \ge \frac{1}{4}$ and we have

$$ \frac{\varepsilon}{2}(\phi'_{\varepsilon})^2 \geq \frac{1}{4}\phi_{\varepsilon}^4 - \frac{1}{2}\phi_{\varepsilon}^2 + \frac{1}{4} = \frac{1}{4}(\phi_{\varepsilon}^2 - 1)^2, $$

hence

$$ \phi'_{\varepsilon} \ge \frac{1}{\sqrt{2\varepsilon}} (\phi_{\varepsilon}^2 - 1) \quad \text{on } [y_{\eta}, Y_{\varepsilon}). $$

Finally, we consider the ODE

$$ \xi'_{\varepsilon, \eta} = \frac{1}{\sqrt{2\varepsilon}}(\xi_{\varepsilon, \eta}^2 - 1), \quad \xi_{\varepsilon, \eta}(y_{\eta}) = 1 + \eta, $$

whose solution is

$$ \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)}, $$

which is defined only up to $Y_{\varepsilon, \eta} := \sqrt{\frac{\varepsilon}{2}}\log(1 + \frac{2}{\eta}) + y_{\eta}$. Applying again Theorem 7, we have that $\phi_{\varepsilon}(y) \ge \xi_{\varepsilon, \eta}(y)$ for any $y \ge 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} \le Y_{\varepsilon, \eta} < \infty$ and that

$$ \forall y \in [y_{\eta}, Y_{\varepsilon}), \quad \psi_{\varepsilon}(y) \ge \xi_{\varepsilon, \eta}(y). \qquad (22) $$

These results are illustrated on Figure 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}$.
$$ \epsilon = 0.1, \mu = 0.5, \eta = 0.5 $$

<figure><img src="image_9.png" /><figcaption>FIGURE 7 – Description of $X_μ$ for $μ = 0.5$, along with the plot of $\psi_ε$ and the lower bound $\xi_{ε,η}$ for $ε = 0.1, η = 0.5$. While $y < B_0$, $\psi_ε$ can possibly oscillate around $\psi_0$, but as soon as $y \ge B_0$, $\psi_ε$ is strictly convex and has no other choice than to explode in finite time $Y_ε \le Y_{ε,η}$, where $Y_{ε,η}$ is the explosion time of the lower bound $\xi_{ε,η}$.</figcaption></figure>

REFERENCES

[1] I. Babuška and J. Osborn. Eigenvalue problems. In *Handbook of Numerical Analysis*, volume 2 of *Finite Element Methods (Part 1)*, pages 641–787. Elsevier, 1991.

[2] C. M. Bender and S. A. Orszag. *Advanced Mathematical Methods for Scientists and Engineers. 1: Asymptotic Methods and Perturbation Theory*. Springer, New York Heidelberg, 1999.

[3] S. Bernstein. Sur la nature analytique des solutions des équations aux dérivées partielles du second ordre. *Mathematische Annalen*, 59(1-2):20–76, 1904.

[4] S. Blatt. On the analyticity of solutions to non-linear elliptic partial differential systems. *arXiv:2009.08762 [math.AP]*, 2020.

[5] P. E. Blöchl. Projector augmented-wave method. *Physical Review B*, 50(24):17953–17979, 1994.

[6] E. Cancès, R. Chakir, and Y. Maday. Numerical Analysis of Nonlinear Eigenvalue Problems. *Journal of Scientific Computing*, 45(1):90–117, 2010.

[7] E. Cancès, R. Chakir, and Y. Maday. Numerical analysis of the planewave discretization of some orbital-free and Kohn-Sham models. *ESAIM: Mathematical Modelling and Numerical Analysis*, 46(2):341–388, 2012.

[8] A. Friedman. On the Regularity of the Solutions of Non-Linear Elliptic and Parabolic Systems of Partial Differential Equations. *Indiana University Mathematics Journal*, 7(1):43–59, 1958.

[9] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch. QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials. *Journal of Physics. Condensed Matter: An Institute of Physics Journal*, 21(39):395502, 2009.

[10] S. Goedecker, M. Teter, and J. Hutter. Separable dual-space Gaussian pseudopotentials. *Physical Review B*, 54(3):1703, 1996.

[11] X. Gonze, B. Amadon, G. Antonius, F. Arnardi, L. Baguet, J.-M. Beuken, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, N. Brouwer, F. Bruneval, G. Brunin, T. Cavignac, J.-B. Charraud, W. Chen, M. Côté, S. Cottenier, J. Denier, G. Geneste, P. Ghosez, M. Giantomassi, Y. Gillet, O. Gingras, D. R. Hamann, G. Hautier, X. He, N. Helbig, N. Holzwarth, Y. Jia, F. Jollet,
W. Lafargue-Dit-Hauret, K. Lejaeghere, M. A. L. Marques, A. Martin, C. Martins, H. P. C. Miranda, F. Naccarato, K. Persson, G. Petretto, V. Planes, Y. Pouillon, S. Prokhorenko, F. Ricci, G.-M. Rignanese, A. H. Romero, M. M. Schmitt, M. Torrent, M. J. van Setten, B. Van Troeye, M. J. Verstraete, G. Zérah, and J. W. Zwanziger. The Abinit project: Impact, environment and recent developments. *Computer Physics Communications*, 248:107042, 2020.

[12] C. Hartwigsen, S. Goedecker, and J. Hutter. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. *Physical Review B*, 58(7):3641-3662, 1998.

[13] Y. Hashimoto. A Remark on the Analyticity of the Solutions for Non-Linear Elliptic Partial Differential Equations. *Tokyo Journal of Mathematics*, 29(2):271-281, 2006.

[14] M. F. Herbst, A. Levitt, and E. Cancès. DFTK: A Julian approach for simulating electrons in solids. *Proceedings of the JuliaCon Conferences*, 3(26):69, 2021.

[15] D. C. Liu and J. Nocedal. On the limited memory BFGS method for large scale optimization. *Mathematical Programming*, 45(1-3):503-528, Aug. 1989.

[16] C. B. Morrey. On the Analyticity of the Solutions of Analytic Non-Linear Elliptic Systems of Partial Differential Equations: Part I. Analyticity in the Interior. *American Journal of Mathematics*, 80(1):198, 1958.

[17] C. B. Morrey. On the Analyticity of the Solutions of Analytic Non-Linear Elliptic Systems of Partial Differential Equations: Part II. Analyticity at the Boundary. *American Journal of Mathematics*, 80(1):219, 1958.

[18] T. Nottoli, I. Gianni, A. Levitt, and F. Lipparini. A robust, open-source implementation of the locally optimal block preconditioned conjugate gradient for large eigenvalue problems in quantum chemistry. *Theoretical Chemistry Accounts*, 142(8):69, Aug. 2023.

[19] J. P. Perdew and Y. Wang. Accurate and simple analytic representation of the electron-gas correlation energy. *Physical Review B*, 45(23):13244-13249, 1992.

[20] I. G. Petrovskii. Sur l'analyticité des solutions des systèmes d'équations différentielles. *Matematiceskij sbornik*, 47(1):3-70, 1939.

[21] L. E. Ratcliff, W. Dawson, G. Fisicaro, D. Caliste, S. Mohr, A. Degomme, B. Videau, V. Cristiglio, M. Stella, M. D'Alessandro, S. Goedecker, T. Nakajima, T. Deutsch, and L. Genovese. Flexibilities of wavelets as a computational basis set for large-scale electronic structure calculations. *Journal of Chemical Physics*, 152(19):194110, 2020.

[22] M. Reed and B. Simon. *Analysis of Operators*. Number 4 in Methods of Modern Mathematical Physics. Academic Press, 1978.

[23] A. H. Romero, D. C. Allan, B. Amadon, G. Antonius, T. Applencourt, L. Baguet, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, F. Bruneval, G. Brunin, D. Caliste, M. Côté, J. Denier, C. Dreyer, P. Ghosez, M. Giantomassi, Y. Gillet, O. Gingras, D. R. Hamann, G. Hautier, F. Jollet, G. Jomard, A. Martin, H. P. C. Miranda, F. Naccarato, G. Petretto, N. A. Pike, V. Planes, S. Prokhorenko, T. Rangel, F. Ricci, G.-M. Rignanese, M. Royo, M. Stengel, M. Torrent, M. J. van Setten, B. Van Troeye, M. J. Verstraete, J. Wiktor, J. W. Zwanziger, and X. Gonze. ABINIT: Overview and focus on selected capabilities. *Journal of Chemical Physics*, 152(12):124102, 2020.

[24] J. C. Slater. A Simplification of the Hartree-Fock Method. *Physical Review*, 81(3):385-390, 1951.

[25] N. Troullier and J. L. Martins. Efficient pseudopotentials for plane-wave calculations. *Physical Review B*, 43(3):1993-2006, 1991.

[26] R. Vrabel, P. Tanuska, P. Vazan, P. Schreiber, and V. Liska. Duffing-Type Oscillator with a Bounded from above Potential in the Presence of Saddle-Center Bifurcation and Singular Perturbation: Frequency Control. *Abstract and Applied Analysis*, 2013:1-7, 2013.

[27] W. Walter. *Ordinary Differential Equations*. Number 182 in Graduate Texts in Mathematics; Readings in Mathematics. Springer, New York, 1998.

[28] T. Wazewski. Systèmes des équations et des inégalités différentielles ordinaires aux deuxièmes membres monotones et leurs applications. *Annales de la Société polonaise de mathématique*, 23:112-166, 1950.