# Variational projector-augmented wave method: a full-potential approach for
electronic structure calculations in solid-state physics

Mi-Song Dupuy Zentrum Mathematik, Technische Universität München,
Boltzmannstraße 3, 85747 Garching, Germany
(Email<EMAIL_ADDRESS>
###### Abstract

In solid-state physics, energies of crystals are usually computed with a
plane-wave discretization of Kohn-Sham equations. However the presence of
Coulomb singularities requires the use of large plane-wave cut-offs to produce
accurate numerical results. In this paper, an analysis of the plane-wave
convergence of the eigenvalues of periodic linear Hamiltonians with Coulomb
potentials using the variational projector-augmented wave (VPAW) method is
presented. In the VPAW method, an invertible transformation is applied to the
original eigenvalue problem, acting locally in balls centered at the
singularities. In this setting, a generalized eigenvalue problem needs to be
solved using plane-waves. We show that cusps of the eigenfunctions of the VPAW
eigenvalue problem at the positions of the nuclei are significantly reduced.
These eigenfunctions have however a higher-order derivative discontinuity at
the spheres centered at the nuclei. By balancing both sources of error, we
show that the VPAW method can drastically improve the plane-wave convergence
of the eigenvalues with a minor additional computational cost. Numerical tests
are provided confirming the efficiency of the method to treat Coulomb
singularities.

Keywords. Eigenvalue problems, Spectral method, Error analysis, Weighted
Sobolev space.

## Introduction

In solid-state physics, the electronic properties of a material are given by
the bottom of the spectrum of a Kohn-Sham Hamiltonian. The spectrum of this
Hamiltonian can be approximated by solving eigenvalue problems arising from
the Bloch transform. To take advantage of the periodicity of the system, a
plane-wave discretization is a natural choice. However the Coulomb potentials
in the Hamiltonian considerably impede the convergence rate of Fourier methods
because of the cusps [18] of the eigenfunctions located at each nucleus.

Two approaches to deal with the Coulomb singularities can be distinguished. In
the pseudopotential approach, the potential is regularized in a neighborhood
of each nucleus, using a process which preserves the bottom of the spectrum of
specific Hamiltonians (see [6, 10] for more details on the generation of
pseudopotentials). In doing so, a modelling error is introduced which is
rarely estimated. Another way to handle the singularities is to use a
different basis functions. Since the Coulomb potentials are point
singularities, one can modify the plane-wave basis locally around the nuclei
and use functions that are less sensitive to those singularities. This is the
main idea of the augmented plane-wave (APW) method [30] or of the linearised
augmented plane-wave (LAPW) method [19]. As opposed to the pseudopotential
approach, these methods dealing directly with Coulomb singularities are called
_full-potential_ approach. The APW and LAPW methods can be reinterpreted as
discontinuous Galerkin methods and variants of these methods have been studied
in [9, 25].

The purpose of this paper is to provide an analysis of the variational
projector-augmented wave (VPAW) method for the plane-wave discretization of
periodic linear Hamiltonians with Coulomb potentials. The method has been
introduced in [2]. Its efficiency for the plane-wave discretization of one-
dimensional toy model with Dirac potentials has been shown [3]. The VPAW
method is based on the projector-augmented wave (PAW) method [4] which has
become a state-of-the-art method in computational solid-state physics and
implemented in several popular electronic structure simulation codes (Abinit
[32], VASP [21]). The idea of the PAW method is to modify the original
eigenvalue problem by an invertible transformation. This transformation acts
locally in the neighborhood of each nucleus and maps atomic wave functions to
smooth functions called pseudo wave functions. Moreover the PAW treatment
allows for the use of pseudopotentials without any approximation. In practice
however, the PAW method yields equations with infinite expansions that have to
be truncated. This induces an error that has been analyzed in [11] for the
same one-dimensional model mentioned previously. In the VPAW method, the
invertible transformation is a finite sum of operators, hence no truncation is
needed and no approximation is made. Although Coulomb singularities are still
present in the equations that are discretized, an acceleration of the plane-
wave convergence can be proved.

The paper is organized as follows. The VPAW method applied to a periodic
linear Hamiltonian is presented in Section 1. Estimates on the eigenvalues of
the plane-wave discretization of the VPAW equations are given in Section 2.
Numerical tests confirming the efficiency of the VPAW method on a simple model
can be found in Section 3. Proofs of the acceleration of convergence are
gathered in Section 5.

##  1 The VPAW method

###  1.1 The eigenvalue problem

For simplicity, we restrict ourselves to the linear model. A quick overview of
the spectral theory of periodic Hamiltonians can be found in [15]. More
thoroughful expositions of this theory are presented in [12, 23]. For
extensions to nonlinear equations, the interested reader is referred to [8]
for the Hartree model and to [7] for the Hartree-Fock model.

The crystal is modelled as an infinite periodic motif of $N_{at}$ point
charges at positions $\mathbf{R}_{I}$ in the unit cell

$\Gamma=\left\\{\alpha_{1}\mathbf{a}_{1}+\alpha_{2}\mathbf{a}_{2}+\alpha_{3}\mathbf{a}_{3},(\alpha_{1},\alpha_{2},\alpha_{3})\in[-1/2,1/2)^{3}\right\\}.$

and repeated over a periodic lattice

$\mathcal{R}=\mathbb{Z}\mathbf{a}_{1}+\mathbb{Z}\mathbf{a}_{2}+\mathbb{Z}\mathbf{a}_{3}.$

where $\mathbf{a}_{1},\mathbf{a}_{2},\mathbf{a}_{3}$ are linearly independent
vectors of $\mathbb{R}^{3}$.

The electronic properties of the crystal are determined by the spectrum of the
periodic Hamiltonian $H_{\mathrm{per}}$ acting on $L^{2}(\mathbb{R}^{3})$:

$H_{\mathrm{per}}=-\frac{1}{2}\Delta+V_{\mathrm{per}}+W_{\mathrm{per}},$ (1.1)

where $V_{\mathrm{per}}$ is an $\mathcal{R}$-periodic potential defined by

$\begin{cases}-\Delta
V_{\mathrm{per}}=4\pi\left(\sum\limits_{\mathbf{T}\in\mathcal{R}}\sum\limits_{I=1}^{N_{at}}Z_{I}\left(\delta_{\mathbf{R}_{I}}(\cdot+\mathbf{T})-\frac{1}{|\Gamma|}\right)\right)\\\
V_{\mathrm{per}}\text{ is }\mathcal{R}\text{-periodic}.\end{cases}$ (1.2)

In this paper, $W_{\mathrm{per}}$ is a smooth $\mathcal{R}\text{-periodic}$
potential so that Equation (1.2) has a solution. In practice,
$W_{\mathrm{per}}$ is a nonlinear potential depending on the model chosen to
describe the electronic self-interaction (typically a Kohn-Sham potential).

The standard way to study the spectrum of $H_{\mathrm{per}}$ is through Bloch
theory which will be outlined in the next few lines. Let $\mathcal{R}^{*}$ be
the dual lattice

$\mathcal{R}^{*}=\mathbb{Z}\mathbf{a}_{1}^{*}+\mathbb{Z}\mathbf{a}_{2}^{*}+\mathbb{Z}\mathbf{a}_{3}^{*},$

where $(\mathbf{a}_{1}^{*},\mathbf{a}_{2}^{*},\mathbf{a}_{3}^{*})$ satisfies
$\mathbf{a}_{i}\cdot\mathbf{a}_{j}^{*}=2\pi\delta_{ij}$. The reciprocal unit
cell is defined by

$\Gamma^{*}=\left\\{\alpha_{1}\mathbf{a}_{1}^{*}+\alpha_{2}\mathbf{a}_{2}^{*}+\alpha_{3}\mathbf{a}_{3}^{*},(\alpha_{1},\alpha_{2},\alpha_{3})\in[-1/2,1/2)^{3}\right\\}.$

As $H_{\mathrm{per}}$ commutes with $\mathcal{R}$-translations,
$H_{\mathrm{per}}$ admits a Bloch decomposition in operators $H_{\mathbf{q}}$
acting on

$L^{2}_{\mathrm{per}}(\Gamma)=\\{f\in L^{2}_{\mathrm{loc}}(\mathbb{R}^{3})\ |\
f\ \text{is }\mathcal{R}\text{-periodic}\\},$

with domain

$H^{2}_{\mathrm{per}}(\Gamma)=\\{f\in H^{2}_{\mathrm{loc}}(\mathbb{R}^{3})\ |\
f\ \text{is }\mathcal{R}\text{-periodic}\\}.$

The operator $H_{\mathbf{q}}$ is given by:

$H_{\mathbf{q}}=\frac{1}{2}|-i\nabla+\mathbf{q}|^{2}+V_{\mathrm{per}}+W_{\mathrm{per}},\quad\mathbf{q}\in\Gamma^{*}.$

For each $\mathbf{q}\in\Gamma^{*}$, the operator $H_{\mathbf{q}}$ is self-
adjoint, bounded below and with compact resolvent. Thus it has a discrete
spectrum of infinite eigenvalues $E_{1,\mathbf{q}}\leq
E_{2,\mathbf{q}}\leq\dots\leq E_{n,\mathbf{q}}\to\infty$, counted with
multiplicities, and the associated eigenfunctions
$(\psi_{n,\mathbf{q}})_{n\in\mathbb{N}^{*}}$ form an orthonormal basis of
$L^{2}_{\mathrm{per}}(\Gamma)$:

$H_{\mathbf{q}}\psi_{n,\mathbf{q}}=E_{n,\mathbf{q}}\psi_{n,\mathbf{q}}.$ (1.3)

By Bloch theorem [29, Chapter XIII], the spectrum of $H_{\mathrm{per}}$ is
given by the union of the discrete spectra of an infinite number of eigenvalue
problems parameterized by the vector $\mathbf{q}$ belonging to the reciprocal
unit cell $\Gamma^{*}$:

$\sigma(H_{\mathrm{per}};L^{2}(\mathbb{R}^{3}))=\bigcup\limits_{\mathbf{q}\in\Gamma^{*}}\sigma(H_{\mathbf{q}};L^{2}_{\mathrm{per}}(\Gamma)).$
(1.4)

The VPAW method aims to ease the resolution of the eigenvalue problem (1.3).
For clarity, we will only present the case $\mathbf{q}=0$ and denote $H_{0}$
by $H$ as this example contains all the main difficulties encountered in the
numerical resolution of Equation (1.3).

###  1.2 The VPAW method for solids

Following the idea of the PAW method, an invertible transformation
$(\mathrm{Id}+T)$ is applied to the eigenvalue problem (1.3), where $T$ is the
sum of operators $T_{I}$ acting locally around each nucleus. For each operator
$T_{I}$, two parameters $N_{\mathrm{paw}}$ and ${r_{c}}$ need to be fixed
($r_{c}$ and $N_{\mathrm{paw}}$ may depend on the atomic site $I$):

  1. 1.

$N_{\mathrm{paw}}$ is the number of PAW functions used to build $T_{I}$,

  2. 2.

${r_{c}}$ is a cut-off radius which will set the acting domain of $T_{I}$,
more precisely:

     * •

for all $f\in L_{\mathrm{per}}^{2}(\Gamma)$,
$\mathrm{supp}(T_{I}f)\subset\bigcup\limits_{\mathbf{T}\in\mathcal{R}}B(\mathbf{R}_{I}+\mathbf{T},{r_{c}})$,
where $B(\mathbf{R},r)$ is the closed ball of $\mathbb{R}^{3}$ with center
$\mathbf{R}$ and radius $r$,

     * •

if
$\mathrm{supp}(f)\bigcap\bigcup\limits_{\mathbf{T}\in\mathcal{R}}B(\mathbf{R}_{I}+\mathbf{T},{r_{c}})=\emptyset$,
then $T_{I}f=0$.

The operator $T_{I}$ is given by:

$T_{I}=\sum\limits_{\mathbf{T}\in\mathcal{R}}\sum\limits_{k=1}^{N_{\mathrm{paw}}}(\phi^{I}_{k}(\mathbf{r}-\mathbf{R}_{I})-\widetilde{\phi}^{I}_{k}(\mathbf{r}-\mathbf{R}_{I}))\langle\widetilde{p}^{I}_{k}(\cdot-\mathbf{R}_{I})\,,\bm{\cdot}\rangle,$
(1.5)

where $\langle\bm{\cdot}\,,\bm{\cdot}\rangle$ is the $L^{2}$-scalar product on
the unit cell $\Gamma$ and the functions $\phi^{I}_{k}$,
$\widetilde{\phi}^{I}_{k}$ and $\widetilde{p}^{I}_{k}$ are functions in
$L^{2}_{\mathrm{per}}(\Gamma)$. The PAW functions $(\phi_{k}^{I})_{1\leq k\leq
N_{\mathrm{paw}}}$, $(\widetilde{\phi}_{k}^{I})_{1\leq k\leq
N_{\mathrm{paw}}}$ and $(\widetilde{p}_{k}^{I})_{1\leq k\leq
N_{\mathrm{paw}}}$ must satisfy the following essential properties:

  1. 1.

$\mathrm{supp}\,(\phi_{k}^{I}-\widetilde{\phi}_{k}^{I})\subset\bigcup\limits_{\mathbf{T}\in\mathcal{R}}B(\mathbf{T},{r_{c}})$;

  2. 2.

$\widetilde{\phi}_{k}^{I}$ restricted to $B(0,{r_{c}})$ is smooth;

  3. 3.

$\widetilde{p}_{k}^{I}$ are supported in
$\bigcup\limits_{\mathbf{T}\in\mathcal{R}}B(\mathbf{T},{r_{c}})$ and
$\forall\,1\leq i,j\leq N_{\mathrm{paw}},\
\langle\widetilde{p}^{I}_{i}\,,\widetilde{\phi}^{I}_{j}\rangle=\delta_{ij}$
(_i.e._ $(\widetilde{p}^{I}_{j})_{1\leq j\leq N_{\mathrm{paw}}}$ is dual to
$(\widetilde{\phi}^{I}_{j})_{1\leq j\leq N_{\mathrm{paw}}}$).

The operators $T_{I}$ act locally in
$\bigcup\limits_{\mathbf{T}\in\mathcal{R}}B(\mathbf{R}_{I}+\mathbf{T},{r_{c}})$.

Several schemes exist in the literature to generate the PAW functions. In this
paper, the PAW functions are close to the Vanderbilt scheme [22] where only
the projector functions differ from ours. The Blöchl scheme [4] is another
popular way to generate PAW functions although the first seems to be preferred
[17]. See [3, 17] for more details on the generation of the PAW functions.

#### Atomic wave function

Let $(\varphi^{I}_{k})_{1\leq k\leq
N_{\mathrm{paw}}}\in(L^{2}(\mathbb{R}^{3}))^{N_{\mathrm{paw}}}$ be
eigenfunctions of an atomic _non-periodic_ Hamiltonian

$H_{I}\varphi^{I}_{k}=\epsilon_{k}\varphi^{I}_{k},\quad\epsilon^{I}_{1}\leq\epsilon^{I}_{2}\leq\epsilon^{I}_{3}\leq\dots,\quad\int_{\mathbb{R}^{3}}\varphi^{I}_{k}\varphi^{I}_{k^{\prime}}=\delta_{kk^{\prime}},$

with $H_{I}$ defined by

$H_{I}=-\frac{1}{2}\Delta-\frac{Z_{I}}{|\mathbf{r}|}+W_{\mathrm{at}}(|\mathbf{r}|),$
(1.6)

where $W_{\mathrm{at}}$ is a smooth bounded potential. The operator $H_{I}$ is
self-adjoint on $L^{2}(\mathbb{R}^{3})$ with domain $H^{2}(\mathbb{R}^{3})$.
Again, in practice, $W_{\mathrm{at}}$ is a radial nonlinear potential
belonging to the same family of models as $W_{\mathrm{per}}$ in Equation
(1.2). Since the atomic Hamiltonian is rotationnaly invariant, $H_{I}$ is
block-diagonal in the decomposition of $L^{2}(\mathbb{R}^{3})$ associated with
the eigenspaces of the operator $\mathbf{L}^{2}$ (the square of the angular
momentum $\mathbf{L}=\mathbf{r}\times\mathbf{p}=\mathbf{r}\times(-i\nabla)$).
The eigenfunctions $\varphi_{k}^{I}$ can be decomposed into a radial function
and a spherical harmonics (see [29, Chapter XIII.3.B] for further details):

$\varphi_{k}^{I}(\mathbf{r})=r^{\ell}R_{n\ell}(r)Y_{\ell
m}(\hat{\mathbf{r}}),$ (1.7)

where $Y_{\ell m}$ is the real spherical harmonics and $R_{n\ell}$ is a
continuous function such that $\lim\limits_{r\to 0}|R_{n\ell}(r)|<\infty$. For
$\mathbf{r}\in\mathbb{R}^{3}$, we define
$\hat{\mathbf{r}}:=\frac{\mathbf{r}}{|\mathbf{r}|}$ and when there is no
ambiguity we will denote by $r$ the euclidean norm of $\mathbf{r}$. The
decomposition (1.7) also holds for some nonlinear models, see [31, 5]. The
functions $R_{n\ell}$ satisfies the following radial Schrödinger equation

$\mathfrak{h}_{\ell}R_{n\ell}(r)=-\frac{1}{2}R_{n\ell}^{\prime\prime}(r)-\frac{\ell+1}{r}R_{n\ell}^{\prime}(r)-\frac{Z_{I}}{r}R_{n\ell}(r)+W_{\mathrm{at}}(r)R_{n\ell}(r)=\epsilon_{n\ell}R_{n\ell}(r).$
(1.8)

The eigenvalues of $\mathfrak{h}_{\ell}$, if they exist, are all simple. The
discrete spectrum of $H_{I}$ is then the collection of all the eigenvalues of
the operators $\mathfrak{h}_{\ell}$, $\ell\geq 0$.

The PAW atomic wave functions $(\phi^{I}_{k})_{1\leq k\leq
N_{\mathrm{paw}}}\in(L^{2}_{\mathrm{per}}(\Gamma))^{N_{\mathrm{paw}}}$ are
then defined by

  * •

for $1\leq k\leq N_{\mathrm{paw}}$ and $\mathbf{r}\in\Gamma$,
$\phi^{I}_{k}(\mathbf{r})=\varphi^{I}_{k}(\mathbf{r})$,

  * •

$\phi^{I}_{k}$ is $\mathcal{R}$-periodic.

If $W_{\mathrm{at}}\not=0$, there is a minimal angular momentum
$\ell_{\mathrm{adm}}$ for which $\mathfrak{h}_{\ell}$ for all
$\ell\geq\ell_{\mathrm{adm}}$ has no eigenvalue (see [29, Theorem XIII.8] for
the linear case and [31] for the reduced Hartree-Fock extension). As an
immediate consequence, PAW functions can only be selected for a finite range
of angular momentum $\ell\leq\ell_{\mathrm{adm}}$.

We denote by $(n_{0},n_{1},\dots,n_{\ell_{\mathrm{max}}})$ the number of PAW
functions for each admissible angular momentum, _i.e._ there are $n_{0}$ PAW
functions for the angular momentum $\ell=0,m=0$, $n_{1}$ PAW functions for
$\ell=1$, $|m|\leq 1$, …The total number of PAW functions for one atomic site
is thus given by
$N_{\mathrm{paw}}=\sum\limits_{\ell=0}^{\ell_{\mathrm{max}}}(2\ell+1)n_{\ell}$.

#### Pseudo wave function

The pseudo wave functions $\widetilde{\phi}_{k}^{I}$ are the
$\mathcal{R}$-periodic functions given in the unit cell $\Gamma$ by:

$\forall\,\mathbf{r}\in\Gamma,\ \widetilde{\phi}_{k}^{I}({\bf
r})=r^{\ell}\widetilde{R}_{n\ell}(r)\,Y_{\ell m}(\hat{\mathbf{r}}).$ (1.9)

where $k$ stands for the multiple index $(n,\ell,m)$. The radial functions
$\widetilde{R}_{n\ell}$, $1\leq n\leq
n_{\ell},0\leq\ell\leq\ell_{\mathrm{max}}$ are polynomial inside the
augmentation region $B(0,{r_{c}})$:

$\widetilde{R}_{n\ell}(r)=\begin{cases}\sum\limits_{k=0}^{d}c_{2k}r^{2k}&\text{for
}0\leq r\leq{r_{c}}\\\ R_{n\ell}(r)&\text{for }r>{r_{c}}\end{cases}$ (1.10)

and the coefficients are chosen to match $R_{n\ell}$ and its first $(d-1)$
derivatives of $R_{n\ell}$ at ${r_{c}}$.

#### Projector functions

The projector functions $(\widetilde{p}_{k}^{I})_{1\leq k\leq
N_{\mathrm{paw}}}$ chosen here are the $\mathcal{R}$-periodic functions given
in the unit cell $\Gamma$ by:

$\forall\,\mathbf{r}\in\Gamma,\ \widetilde{p}^{I}_{n\ell
m}(\mathbf{r})=r^{\ell}p_{n\ell}(r)Y_{\ell m}(\hat{\mathbf{r}}).$ (1.11)

The functions $p_{n\ell}$ for $0\leq\ell\leq\ell_{\mathrm{max}},1\leq n\leq
n_{\ell}$ are defined by

$p_{n\ell}(r)=\sum\limits_{n^{\prime}=1}^{n_{\ell}}\left(B_{\ell}^{-1}\right)_{nn^{\prime}}\chi(r)\widetilde{R}_{n^{\prime}\ell}(r),$
(1.12)

with $\chi$ a smooth positive cut-off function supported in $(0,{r_{c}})$ and

$B_{\ell}=\left(\int_{0}^{{r_{c}}}\chi(r)\widetilde{R}_{n\ell}(r)\widetilde{R}_{n^{\prime}\ell}(r)r^{2+2\ell}\,\mathrm{d}r\right)_{1\leq
n,n^{\prime}\leq n_{\ell}}.$ (1.13)

By definition, the projector functions $(\widetilde{p}^{I}_{k})_{1\leq k\leq
N_{\mathrm{paw}}}$ are supported in
$\bigcup\limits_{\mathbf{T}\in\mathcal{R}}B(\mathbf{T},{r_{c}})$ and form a
dual family to the pseudo wave functions $(\widetilde{\phi}^{I}_{k})_{1\leq
k\leq N_{\mathrm{paw}}}$:
$\langle\widetilde{p}^{I}_{k}\,,\widetilde{\phi}^{I}_{k^{\prime}}\rangle=\delta_{kk^{\prime}}$.

The VPAW equations to solve are then:

$H^{\mathrm{VPAW}}\widetilde{\psi}=ES^{\mathrm{VPAW}}\widetilde{\psi},$ (1.14)

where

$H^{\mathrm{VPAW}}=(\mathrm{Id}+T)^{*}H(\mathrm{Id}+T),\quad
S^{\mathrm{VPAW}}=(\mathrm{Id}+T)^{*}(\mathrm{Id}+T),$ (1.15)

and

$T=\sum\limits_{I=1}^{N_{at}}T_{I}.$

Thus if $(\mathrm{Id}+T)$ is invertible, the eigenfunctions of $H$ can be
recovered by the formula

$\psi=(\mathrm{Id}+T)\widetilde{\psi},$ (1.16)

and the eigenvalues are identical to the original eigenvalue problem (1.3).

By construction, the operator $(\mathrm{Id}+T_{I})$ maps the pseudo wave
functions $\widetilde{\phi}^{I}$ to the atomic eigenfunctions $\phi^{I}$:

$(\mathrm{Id}+T_{I})\widetilde{\phi}^{I}_{k}(\cdot-\mathbf{R}_{I})=\phi^{I}_{k}(\cdot-\mathbf{R}_{I}),$

so if locally around each nucleus, the function $\psi$ “behaves” like the
atomic wave functions $\phi^{I}_{k}$, we can hope that the cusp behavior of
$\psi$ is captured by the operator $T$. The VPAW eigenfunction
$\widetilde{\psi}$ would therefore be smoother than $\psi$ and the plane-wave
expansion of $\widetilde{\psi}$ would converge faster than the expansion of
$\psi$.

###  1.3 Well-posedness of the VPAW method

To be well-posed the VPAW method requires

  1. 1.

for each $0\leq\ell\leq\ell_{\mathrm{max}}$, the family of pseudo wave
functions $(\widetilde{R}_{n\ell})_{1\leq n\leq n_{\ell}}$ to be linearly
independent in $[0,r_{c}]$, so that the projector functions
$(p_{n\ell})_{1\leq n\leq n_{\ell}}$ are well defined;

  2. 2.

$(\mathrm{Id}+T)$ to be invertible.

To fulfill the first condition, the following assertion is assumed.

######  Assumption 1.

For all $0<r_{c}<r_{\mathrm{min}}$ and each
$0\leq\ell\leq\ell_{\mathrm{max}}$, $({\mathcal{R}}^{(k)}(r_{c}))_{0\leq k\leq
n_{\ell}-1}$ is a linearly independent family, where ${\mathcal{R}}$ is the
vector of the functions $(R_{1\ell},\dots,R_{n_{\ell}\ell})$.

This condition ensures that the family of pseudo-wave functions
$\left(\widetilde{R}_{n\ell}\right)_{1\leq n\leq n_{\ell}}$ is linearly
independent. This assumption holds in the particular case of the hydrogenoid
atom (see Lemma A.4 in the appendix).

It can be shown that the second condition is equivalent to the invertibility
of the matrix $\left(\langle
p_{j\ell}\,,r^{\ell}R_{k\ell}\rangle\right)_{1\leq j,k\leq n_{\ell}}$ for each
$0\leq\ell\leq\ell_{\mathrm{max}}$. Since the proof of this statement is very
close to the proof of Proposition 2.3 in [3], we will not reproduce it here.
For the rest of the paper, we make the following assumption.

######  Assumption 2.

For all $0<r_{c}<r_{\mathrm{min}}$ and any $0\leq\ell\leq\ell_{\mathrm{max}}$,
the matrix $\left(\langle p_{n\ell}\,,R_{n^{\prime}\ell}\rangle\right)_{1\leq
n,n^{\prime}\leq n_{\ell}}$ is invertible.

Finally, it appears in the proof of Theorem 2.1 that we need the following
assumption. It is shown in the appendix that this assumption is satisfied for
the hydrogenoid eigenfunctions.

######  Assumption 3.

For every $0\leq\ell\leq\ell_{\mathrm{max}}$ and $0\leq n\leq n_{\ell}$, the
radial functions $R_{n\ell}$ defined in (1.7) satisfies $R_{n\ell}(0)\not=0$.

There is a natural interpretation to this condition. The poor convergence of
plane-wave expansions in the computation of the eigenvalues of (1.3) is due to
the cusps of the molecular wave functions. Hence we want to use atomic wave
functions with cusps to reduce them. By the Kato cusp condition (see Theorem
1.6 below), this is equivalent to $R_{n\ell}(0)\not=0$.

###  1.4 Computational cost of the VPAW method

A detailed analysis of the computational cost of the PAW method can be found
in [24]: the cost scales like
${\mathcal{O}}(N_{\mathrm{at}}N_{\mathrm{paw}}M+N_{\mathrm{at}}N_{\mathrm{paw}}^{2}+M\log
M)$ where $N_{\mathrm{at}}$ is the number of nuclei, $N_{\mathrm{paw}}$ is the
number of PAW functions per atomic site and $M$ the number of plane-waves.
Usually, $N_{\mathrm{paw}}$ is chosen relatively small, but $M$ may be large,
so it is important to avoid a computational cost of order $M^{2}$.

In practice, we are interested in the cost of the computation of
$H^{\mathrm{VPAW}}\widetilde{\psi}$ and $S^{\mathrm{VPAW}}\widetilde{\psi}$
where $\widetilde{\psi}$ is expanded in $M$ plane-waves as the generalized
eigenvalue problem is solved by a conjugate gradient algorithm. We will only
focus on $H^{\mathrm{VPAW}}\widetilde{\psi}$ since the analysis
$S^{\mathrm{VPAW}}\widetilde{\psi}$ is similar. Let us split
$H^{\mathrm{VPAW}}$ into four terms:

$H^{\mathrm{VPAW}}\widetilde{\psi}=H\widetilde{\psi}+PD_{H}P^{T}\widetilde{\psi}+H\left(\Phi-\widetilde{\Phi}\right)P^{T}\widetilde{\psi}+PH\left(\Phi-\widetilde{\Phi}\right)^{T}\widetilde{\psi},$

where $P$ is the $M\times N_{\mathrm{at}}N_{\mathrm{paw}}$ matrix of the
projector functions, $H(\Phi-\widetilde{\Phi})$ the $M\times
N_{\mathrm{at}}N_{\mathrm{paw}}$ matrix of the Fourier representation of the
$N_{\mathrm{at}}N_{\mathrm{paw}}$ functions
$H(\phi^{I}_{i}-\widetilde{\phi}^{I}_{i})$, and $D_{H}$ is the
$N_{\mathrm{at}}N_{\mathrm{paw}}\times N_{\mathrm{at}}N_{\mathrm{paw}}$ matrix
$\langle\phi^{I}_{i}-\widetilde{\phi}^{I}_{i}\,,H(\phi^{J}_{j}-\widetilde{\phi}^{J}_{j})\rangle$.
Note that $D_{H}$ is a block diagonal matrix, because balls of radius $r_{c}$
centered at different atomic site do not overlap.

The computational cost can be estimated as follows (the cost at each step is
given in brackets):

  1. 1.

$H\widetilde{\psi}$ is assembled in two steps. First,
$-\frac{1}{2}\Delta\widetilde{\psi}$ is computed in ${\mathcal{O}}(M)$ since
the operator $\frac{1}{2}\Delta$ is diagonal in Fourier representation. For
the potential $V$, apply an inverse FFT to $\widetilde{\psi}$ to have the real
space representation of $\widetilde{\psi}$, multiply pointwise by $V$ and
apply a FFT to the whole result (${\mathcal{O}}(M\log M)$);

  2. 2.

for $PD_{H}P^{T}\widetilde{\psi}$, compute the
$N_{\mathrm{at}}N_{\mathrm{paw}}$ projections $P^{T}\widetilde{\psi}$
(${\mathcal{O}}(MN_{\mathrm{at}}N_{\mathrm{paw}})$), then successively apply
the matrices $D_{H}$ (${\mathcal{O}}(N_{\mathrm{at}}N_{\mathrm{paw}}^{2})$
since $D_{H}$ is block-diagonal) and $P$
(${\mathcal{O}}(MN_{\mathrm{at}}N_{\mathrm{paw}})$);

  3. 3.

for $PH(\Phi-\widetilde{\Phi})^{T}\widetilde{\psi}$, similarly apply
successively $H(\Phi-\widetilde{\Phi})^{T}$ to $\widetilde{\psi}$
(${\mathcal{O}}(MN_{\mathrm{at}}N_{\mathrm{paw}})$) and $P$ to
$H(\Phi-\widetilde{\Phi})^{T}\widetilde{\psi}$
(${\mathcal{O}}(MN_{\mathrm{at}}N_{\mathrm{paw}})$);

  4. 4.

for $H(\Phi-\widetilde{\Phi})P^{T}\widetilde{\psi}$, we proceed as in step 3.

Thus, the total numerical cost is of order
${\mathcal{O}}(MN_{\mathrm{at}}N_{\mathrm{paw}}+N_{\mathrm{at}}N_{\mathrm{paw}}^{2}+M\log
M)$ which is the same as for the PAW method.

###  1.5 Singular expansion

It appears that the theory of weighted Sobolev spaces and the singular
expansion of eigenfunctions of Hamiltonians with Coulomb potentials provides a
nice framework to study the Fourier decay of the VPAW pseudo wave functions
$\widetilde{\psi}$. The singular expansion gives a generalization of the Kato
cusp condition [18] to any order. This theory is closely linked to the
$b$-calculus of pseudodifferential operators developed by Melrose [28]. It has
been applied successfully to characterize precisely the behaviour of the
electronic wave function close the nucleus [14, 16] and used in the analysis
of the muffin-tin and LAPW methods [9]. It is also a key element of the
analysis of the convergence of _hp_ -finite elements approximation for similar
models [26, 27]. The interested reader may refer to [20, 13] for a detailed
exposition of this theory.

For simplicity and without loss of generality, we assume that $\Gamma$ is the
cube $[-\frac{1}{2},\frac{1}{2}]^{3}$. We denote by
$\mathcal{S}\subset\mathbb{R}^{3}$ the set of the positions of the nuclei

$\mathcal{S}=\\{\mathbf{R}_{I}+\mathbf{T},\ I=1,\dots,N_{\mathrm{at}},\
\mathbf{T}\in\mathcal{R}\\}.$

Let $\chi$ be a $\mathcal{R}$-periodic continuous function such that
$\chi(\mathbf{R}_{I}+\mathbf{r})=r$ for small $r$, $\chi\in
C^{\infty}_{\mathrm{loc}}(\mathbb{R}^{3}\setminus\mathcal{S})$.

######  Definition 1.1.

Let $k\in\mathbb{N}$ and $\gamma\in\mathbb{R}$. We define the $k$-th weighted
Sobolev space with index $\gamma$ by

${\mathcal{K}}^{k,\gamma}(\Gamma)=\left\\{u\in
L^{2}_{\mathrm{per}}(\Gamma):\chi^{|\alpha|-\gamma}\partial^{\alpha}u\in
L^{2}_{\mathrm{per}}(\Gamma)\ \forall\ |\alpha|\leq k\right\\}.$ (1.17)

Consider a subspace of functions with the asymptotic expansions

$\forall
I=1,\dots,N_{\mathrm{at}},u(\mathbf{r}+\mathbf{R}_{I})\sim\sum\limits_{j\in\mathbb{N}}c^{I}_{j}(\hat{\mathbf{r}})r^{j}\
\text{as }r\to 0,$ (1.18)

where $c^{I}_{j}$ belongs to the finite dimensional subspace
$M_{j}=\mathrm{span}\\{Y_{\ell m},0\leq\ell\leq j,|m|\leq\ell\\}$.

######  Definition 1.2.

Let $k\in\mathbb{N}$ and $\gamma\in\mathbb{R}$. We define the weighted Sobolev
spaces with asymptotic type (1.18) by

$\begin{split}\mathscr{K}^{k,\gamma}(\Gamma)=\Bigg{\\{}u\in\mathcal{K}^{k,\gamma}(\Gamma)\,\bigg{|}\,\eta_{N}\in\mathcal{K}^{k,\gamma+N+1}(\Gamma)\text{
where }\eta_{N}\text{ is the }\Gamma\text{-periodic function defined in
}\Gamma\text{ by}\\\ \left.\forall N\in\mathbb{N},\
\forall\,\mathbf{r}\in\Gamma,\
\eta_{N}(\mathbf{r})=u(\mathbf{r})-\sum\limits_{I=1}^{N_{\mathrm{at}}}\omega(|\mathbf{r}-\mathbf{R}_{I}|)\sum\limits_{j=0}^{N}c_{j}^{I}(\widehat{\mathbf{r-R}_{I}})|\mathbf{r-R}_{I}|^{j}\right\\},\end{split}$
(1.19)

where $\omega$ is a smooth positive cutoff function, _i.e._ $\omega=1$ near 0
and $\omega=0$ outside some neighbourhood of $0$.

Definition 1.2 slightly differs from the definition of the weighted Sobolev
space given in [9] (Equation (2.6)). However, our definition is consistent
with the results that can be found in [16] (see Theorem I.1) and the original
paper [14] (see Proposition 1) from which the definition appearing in [9] is
taken.

The expansion (1.18) can be viewed as a “regularity expansion”. Let us suppose
that the functions $c_{j}$ in the singular expansion are constant. Then all
the even terms appearing in (1.19) are smooth since for any $k\in\mathbb{N}$,
$r\mapsto r^{2k}$ is smooth. For the odd terms in the expansion, the function
$r\mapsto r$ is continuous but not differentiable at the origin, the function
$r\mapsto r^{3}$ is $C^{2}$ but not $C^{3}$ and so on. Since the decay of the
Fourier coefficients depends on the regularity of the function, this expansion
enables one to characterize precisely this decay. A precise estimation of this
decay for all the terms appearing in (1.19) is given in Lemma 5.4 below.

######  Definition 1.3.

A function $u$ is asymptotically well-behaved if
$u\in\mathscr{K}^{\infty,\gamma}(\Gamma)$ for $\gamma<3/2$.

######  Remark 1.4.

Note that if $u$ is asymptotically well-behaved then by the definition of the
weighted Sobolev space with asymptotic type (1.18), the remainder
$\eta_{N}(\mathbf{r})=u(\mathbf{r})-\sum\limits_{I=1}^{N_{\mathrm{at}}}\omega(|\mathbf{r}-\mathbf{R}_{I}|)\sum\limits_{j=0}^{N}c^{I}_{j}(\widehat{\mathbf{r-R}_{I}})|\mathbf{r-R}_{I}|^{j}$
is in the classical Sobolev space
$H^{5/2+N-\varepsilon}_{\mathrm{per}}(\Gamma)$.

The following result, stated in [16, 9], gives the regularity of the
eigenfunction of (1.3) in terms of the previously defined weighted Sobolev
space.

######  Theorem 1.5 ([16, 9]).

Let $\psi$ be an eigenfunction of $H\psi=E\psi$ where $H$ is defined in (1.3).
Then $\psi$ is asymptotically well-behaved.

Theorem 1.5 enables to characterize precisely the singularity of the
Hamiltonian wave function and generalizes the Kato cusp condition for
eigenfunctions of 3D-Hamiltonians. Let $V^{I}$ be the smooth potential such
that in a neighborhood of $\mathbf{R}_{I}$,

$V_{\mathrm{per}}(\mathbf{r})+W_{\mathrm{per}}(\mathbf{r})=-\frac{Z}{|\mathbf{r}-\mathbf{R}_{I}|}+V^{I}(\mathbf{r}-\mathbf{R}_{I}),$
(1.20)

and denote by $(v^{I}_{k})_{k\geq\ell}$ the coefficients of the Taylor
expansion of

$V^{I}_{\ell m}(r)=\int_{S(0,1)}V^{I}(\mathbf{r})Y_{\ell
m}(\hat{\mathbf{r}})\,\mathrm{d}\hat{\mathbf{r}}.$ (1.21)

######  Theorem 1.6.

Let $\ell\in\mathbb{N}$ and $|m|\leq\ell$. Let $\psi$ be an eigenfunction of
$H\psi=E\psi$ where $H$ is defined in (1.3) and $(\psi_{j\ell m})_{j,\ell\leq
j,|m|\leq\ell}$ be the coefficients of the singular expansion of $\psi$,
_i.e._ for all $\varepsilon>0$,

$\psi(\mathbf{r})-\sum\limits_{I=1}^{N_{\mathrm{at}}}\omega(|\mathbf{r-R}_{I}|)\sum\limits_{j=0}^{N}\sum\limits_{|m|\leq\ell\leq
j}\psi^{I}_{j\ell m}|\mathbf{r-R}_{I}|^{j}Y_{\ell
m}(\widehat{|\mathbf{r-R}_{I}|})\in\mathscr{K}^{\infty,\frac{5}{2}+N-\varepsilon}(\Gamma).$

Let $(v^{I}_{k})_{k\geq\ell}$ be the coefficients of the Taylor expansion of
the function $V^{I}_{\ell m}$ defined in (1.21). Then the sequence
$(\psi^{I}_{j\ell m})_{j\geq\ell}$ satisfies

$\forall\,j\geq\ell,\frac{(j+1)(j+2+2\ell)}{2}\psi^{I}_{j+1,\ell
m}=-Z\psi^{I}_{j\ell m}+(v^{I}*\psi^{I})_{j-1}-E\psi^{I}_{j-1,\ell m},$ (1.22)

where $v^{I}*\psi^{I}$ denotes the convolution

$(v^{I}*\psi^{I})_{k}=\sum\limits_{j=\ell}^{k}v^{I}_{k-j}\psi^{I}_{k\ell m}.$

The proof of this theorem follows from the definition of the weighted Sobolev
space ${\mathcal{K}}^{\infty,\frac{5}{2}+N-\varepsilon}$ and the Coulomb
singularity of the potential. For $\ell=0$, the Kato cusp condition is
recovered since $\psi^{I}_{000}=\psi(\mathbf{R}_{I})$ and
$\psi^{I}_{100}=\left.\frac{\partial}{\partial
r}\right|_{r=0}\int_{S(\mathbf{R}_{I},1)}\psi(\mathbf{r})Y_{00}(\hat{\mathbf{r}})\,\mathrm{d}\hat{\mathbf{r}}$.

##  2 Main results

We focus on the analysis of the VPAW method restricted to a set of PAW
functions associated to the lowest angular momentum $\ell=m=0$ (_i.e._
$\ell_{\mathrm{max}}=0$).

From Definition 1.2 and the asymptotic expansion of the molecular wave
function $\psi$, it is possible to identify the cause of the slow decay of the
Fourier coefficients of $\psi$, which is the cusps at each nucleus. We can
show that the cusp of the pseudo wave function $\widetilde{\psi}$ is
significantly reduced by the VPAW method. More precisely (see Proposition
5.1), if $n_{0}$ PAW functions associated to the angular momentum $\ell=m=0$
are used, then there is a constant $C$ independent of $r_{c}$ such that for
all $I=1,\dots,N_{\mathrm{at}}$, for any $0<r_{c}\leq r_{\mathrm{min}}$ and
for all $\varepsilon>0$:

$\left|\frac{\partial}{\partial
r}\bigg{|}_{r=0}\int_{S(\mathbf{R}_{I},1)}\widetilde{\psi}(\mathbf{r})\,\mathrm{d}\hat{\mathbf{r}}\right|\leq
C{r_{c}}^{\min(2n_{0},5)-\varepsilon}.$

The blow-up of the $d$-th derivative introduced at each sphere
$S(\mathbf{R}_{I},r_{c})$ is controlled similarly. It is possible to show (see
Proposition 5.2) that there exists a constant $C$ independent of $r_{c}$ such
that $I=1,\dots,N_{\mathrm{at}}$, for any $0<r_{c}\leq r_{\mathrm{min}}$ and
for all $\varepsilon>0$:

$\left|\left[\int_{S(\mathbf{R}_{I},1)}\widetilde{\psi}^{(d)}(\mathbf{r})\,\mathrm{d}\hat{\mathbf{r}}\right]_{r_{c}}\right|\leq\frac{C}{{r_{c}}^{d-1}}.$

From both estimates, the following plane-wave convergence for the computation
of the eigenvalues with the VPAW method can be proved.

######  Theorem 2.1.

Let $E_{M}$ be an eigenvalue of the variational approximation of (1.15) in a
plane-wave basis with wavenumber $|\mathbf{K}|\leq M$, with $n_{0}$ PAW
functions associated to the angular momentum $\ell=0,m=0$ with smoothness
$d\geq n$ and cut-off radius $r_{c}$. Let $E$ be the corresponding exact
eigenvalue. Under Assumptions 1, 2 and 3, there exists a constant $C>0$
independent of ${r_{c}}$ and $M$ such that for $M$ sufficiently large, for all
$\varepsilon>0$, and for all $0<r_{c}<r_{\mathrm{min}}$, we have

$0<E_{M}-E\leq
C\left(\frac{{r_{c}}^{2\min(2n_{0},5)-2\varepsilon}}{M^{3}}+\frac{{r_{c}}^{\min(2n_{0},5)-\varepsilon}}{M^{4-\varepsilon}}+\frac{1}{{r_{c}}^{2d-2}}\frac{1}{M^{2d-1}}+o\left(\frac{1}{M^{5-\varepsilon}}\right)\right).$
(2.1)

The VPAW method does not erase the cusps appearing in the molecular wave
function $\psi$, hence in the asymptotic regime, the plane-wave convergence
rate is the same as the brute force discretization of the original eigenvalue
problem. The prefactor ${r_{c}}^{2\min(2n_{0},5)-2\varepsilon}$ can be
significantly reduced by taking a small cut-off radius $r_{c}$. However in
doing so, the second prefactor $\frac{1}{{r_{c}}^{2d-2}}$ can become dominant
in the eigenvalue error. Balancing both error terms gives an optimal cut-off
radius equal to $r_{\mathrm{opt}}=(M^{\frac{2d-4}{4n_{0}+2d-2}})^{-1}$. For
$n_{0}=2$ and $d=5$ (which are typical for PAW simulations) and
$r_{c}=r_{\mathrm{opt}}$, both error terms behave like $\frac{1}{M^{6}}$.

Theorem 2.1 holds for _any_ eigenvalue of $H$ given in (1.3), however
numerical tests provided in Section 3 are restricted to the ground-state
eigenvalue. They suggest that the VPAW method can be an efficient strategy to
solve accurately the eigenvalue problem (1.3) (see Figure 1).

######  Remark 2.2.

By incorporating $n_{1}$ functions for each angular momentum $\ell=1$ and
$m\in\\{-1,0,1\\}$, we can improve the convergence estimate (2.1) to

$\forall 0<r_{c}<r_{\mathrm{min}},\ 0<E_{M}-E\leq
C\left(\frac{{r_{c}}^{2\min(2n_{0},5)-2\varepsilon}}{M^{3}}+\frac{{r_{c}}^{\min(2n_{0},5)-\varepsilon}}{M^{4-\varepsilon}}\right.\\\
\left.+\frac{{r_{c}}^{2\min(2n_{1},5)-\varepsilon}}{M^{5}}+\frac{1}{{r_{c}}^{2d-2}}\frac{1}{M^{2d-1}}+o\left(\frac{1}{M^{7-\varepsilon}}\right)\right).$
(2.2)

The only difference between (2.1) and (2.2) is the prefactor of
$\frac{1}{M^{5}}$. In our example (Section 3), improvements for the
computation of ground-state are marginal and visible for large plane-wave
cutoffs (see Figure 1). However, introducing PAW functions for $\ell=1$ might
be beneficial for higher eigenvalues where in a pre-asymtotic regime, the
prefactor of $\frac{1}{M^{5}}$ may be preponderant.

It is interesting to compare the VPAW method convergence with another full-
potential approach like the augmented plane-wave (APW) method [30]. In the APW
method, instead of modifying the Hamiltonian, a different basis set is used
which is not sensitive to the cusps resulting from the Coulomb interaction
with the nuclei. The new basis set is defined by partitioning the unit cell
$\Gamma$ into two types of regions (the so-called _muffin-tin division_):

  1. i)

balls $B(\mathbf{R}_{I},r_{c})$, $I=1,\dots,M$;

  2. ii)

the remaining interstitial region $\mathcal{D}$.

The basis functions consist of augmentations of plane-waves:

$\omega_{\mathbf{K}}(\mathbf{r})=\begin{cases}e^{i\mathbf{K}\cdot\mathbf{r}}&\text{in
}\mathcal{D}\\\
\sum\limits_{\ell=0}^{\ell_{\mathrm{max}}}\sum\limits_{m=-\ell}^{\ell}\alpha_{\ell
m}^{\mathbf{K}}\chi_{\ell}(r_{I})Y_{\ell m}(\widehat{\mathbf{r}_{I}})&\text{in
each }B(\mathbf{R}_{I},r_{c})\end{cases}$ (2.3)

where $\mathbf{r}_{I}=\mathbf{r}-\mathbf{R}_{I}$. The coefficients
$(\alpha_{\ell m}^{\mathbf{K}})_{|m|\leq\ell\leq\ell_{\mathrm{max}}}$ are set
to match the spherical harmonics expansion of $e^{i\mathbf{K}\cdot\mathbf{r}}$
at the boundaries of the balls $B(\mathbf{R}_{I},r_{c})$. These basis
functions are however not continuous at the boundary of the balls
$B(\mathbf{R}_{I},r_{c})$, hence they do not belong to $H^{1}(\Gamma)$: the
APW method is a _nonconforming_ method. In [9], the convergence of the APW
method for a particular choice of $\chi_{\ell}$ is studied where each
$\chi_{\ell}$ is a polynomial of degree less than $N$. The authors showed that
the error on the eigenvalues of the problem (1.3) by the APW method is bounded
by

$\forall s>\frac{3}{2},\
|E^{\mathrm{APW}}_{\eta}-E|\leq\frac{C_{s}}{\eta^{s-\frac{3}{2}}},$

where $\eta=\min(M,\ell_{\mathrm{max}},N)$.

Although this bound holds for any $s>\frac{3}{2}$, the prefactor depends on
$s$ and this dependency is not explicit in the paper. Moreover, in most
situations, $\eta$ is equal to the maximal angular momentum
$\ell_{\mathrm{max}}$. Hence, increasing this parameter is more and more
costly since it introduces $(2\ell+1)(N+1)$ basis functions in the
nonconforming method. On the other hand, the convergence of the VPAW method is
already very good for $n_{0}\leq 2$ PAW functions for $\ell=0$ (see Figure 1).

##  3 Numerical results

In this section, we present some numerical results applied to the Hamiltonian
$H$ in $[-\tfrac{L}{2},\tfrac{L}{2}]^{3}$

$H=-\frac{1}{2}\Delta-\frac{Z}{\left|\mathbf{r}-\frac{\mathbf{R}}{2}\right|}-\frac{Z}{\left|\mathbf{r}+\frac{\mathbf{R}}{2}\right|},$
(3.1)

with periodic boundary conditions. The lowest eigenvalue is sought using
iterative schemes, hence we are interested in the cost of the matrix-vector
multiplication.

The problem is solved using plane-waves. The kinetic operator is diagonal in
the reciprocal space. The potential is discretized using a radial grid around
the nuclei and a uniform grid in the rest of the domain. For the VPAW method,
the following integrals are pre-computed:

  1. 1.

$\langle e_{\mathbf{K}}\,,\widetilde{p}\rangle$ : since
$\widetilde{p}(\mathbf{r})=p(r)Y_{\ell m}(\hat{\mathbf{r}})$ using (5.6),
$\langle e_{\mathbf{K}}\,,\widetilde{p}\rangle$ can be evaluated on a radial
grid.

  2. 2.

$\langle e_{\mathbf{K}}\,,\phi-\widetilde{\phi}\rangle$ : we proceed like for
$\langle e_{\mathbf{K}}\,,\widetilde{p}\rangle$ using a radial grid;

  3. 3.

$\langle e_{\mathbf{K}}\,,H(\phi-\widetilde{\phi})\rangle$ :
$H(\phi-\widetilde{\phi})$ is decomposed into a radial and a non-radial part.
The radial part is evaluated on a radial grid and the non-radial part on a
uniform grid. For non-linear approximations (Hartree-Fock and Kohn-Sham DFT),
this term can be critical since it may be necessary to re-compute these
integrals at each iteration. This is the main drawback of the VPAW method
compared to the PAW method where this term does not exist.

  4. 4.

$\langle\phi-\widetilde{\phi}\,,\phi-\widetilde{\phi}\rangle$,
$\langle\phi-\widetilde{\phi}\,,H(\phi-\widetilde{\phi})\rangle$ : these
integrals are computed using radial grids when possible or using 3D
integration schemes. For nonlinear models, the last integral needs to be
recomputed at each iteration, however, since there are $N_{\mathrm{paw}}^{2}$
of them, it is not too costly.

The numerical results using a Julia [1] homemade code are summarized in the
following figures with $Z=3$, $R=1$ and $L=5$. The atomic PAW function
$\phi_{k}$ are the eigenfunctions of the hydrogenoid atom. For the pseudo
atomic function $\widetilde{\phi}_{k}$, continuity of the function and of the
first four derivatives are enforced (_i.e._ $d=5$). The lowest eigenvalue is
computed using a conjugate-gradient algorithm stopped when the norm of the
residual is less than $10^{-5}$.

For all the plots presented in this section, _VPAW 1s_ denotes the VPAW method
with one PAW function for $\ell=0$ per atom, _VPAW 2s_ with two PAW functions
for $\ell=0$ per atom and _VPAW 2s1p_ with two PAW functions for $\ell=0$ and
one function for $\ell=1$, $|m|\leq 1$ per atom. The reference value for the
lowest eigenvalue is given by the VPAW method for 200 plane waves per
direction. Computation of the reference is out of reach by a direct plane-wave
method. Figures 1, 2 and 4 are log-log plots of the convergence of the lowest
eigenvalue of (3.1) with respect to the number of plane-waves per direction
for different choices of the PAW parameters. Figure 3 is a log-log plot of the
energy difference

In Figure 1, convergence rates of the energy for the VPAW method for a cutoff
radius $r_{c}=\tfrac{R}{2}=0.5$ are presented. We clearly notice that the VPAW
method converges faster than the direct method, gaining up to three orders of
magnitude compared to a brute force plane-wave discretization. The convergence
seems marginally faster when increasing the number of PAW functions per atom.
This can be explained by the relatively large cut-off radius chosen for this
example. We can distinctly see two regimes:

  * •

for $M\leq 20$, in the pre-asymptotic regime, the error on the eigenvalue is
dominated by the term $\frac{1}{{r_{c}}^{2d-2}M^{2d-1}}$ as the estimated
convergence rate suggests (numerically we observe a convergence rate
$M^{-7.8}$ and theoretically $M^{-9}$ is expected). Since the same regularity
of the atomic pseudo wave function $\widetilde{\phi}_{k}$ is used, it is not
surprising to witness a similar behavior in that regime.

  * •

for $M\geq 20$, the convergence rate is close to $M^{-4}$ which is the next
error term given by Theorem 2.1.

Figure 1: Error on the lowest eigenvalue with the VPAW method for different
choices of PAW functions.

The size of the VPAW acting region can significantly impact the convergence
rate in the pre-asymptotic regime (Figure 2). This plot suggests that there
are different phases in the convergence of the VPAW eigenvalue:

  1. 1.

for very low plane-wave cut-off ($M\leq 15$), the VPAW acting region is too
small to be seen by the Fourier grid, hence no improvement is observed;

  2. 2.

as the plane-wave cut-off grows, the VPAW eigenvalue converges very fast,
since for this regime, the prefactors kill the $\frac{1}{M^{3}}$ and
$\frac{1}{M^{4}}$ decay;

  3. 3.

for a larger plane-wave cut-off $M\geq 50$, the convergence slows down since
the prefactor for the $\frac{1}{M^{4}}$ decay is not negligible anymore. Note
that in this regime the error on the eigenvalue decreases as the VPAW cut-off
radius $r_{c}$ is small, in agreement with Theorem 2.1.

Figure 2: Error on the lowest eigenvalue with the VPAW method for different
choices of cut-off radius.

Figure 3 illustrates the dependence on the VPAW cut-off radius $r_{c}$ of the
error for the VPAW 1s method. The numerically observed behavior of the error
is approximately the one predicted by Theorem 2.1 ($r_{c}^{-8}$ for small $M$
and $r_{c}^{2}$ for large $M$)

(a) VPAW 1s with $M=9$

(b) VPAW 1s with $M=28$

Figure 3: Error on the lowest eigenvalue with respect to the cut-off radius
$r_{c}$ for a fixed plane-wave cut-off $M$

In Figure 4, a comparison between the original PAW method and the VPAW method
is provided. Notice that the convergence of the eigenvalue for the PAW method
is not monotone because the limit is below $E_{0}$. If very accurate results
are awaited on the lowest eigenvalue of (3.1), the VPAW method seems the
method of choice.

Figure 4: Error on the lowest eigenvalue with the PAW and VPAW methods.

##  4 Concluding remarks

This paper provides an analysis of the VPAW method for the plane-wave
discretization of eigenvalue problems of periodic linear Hamiltonians with
Coulomb potentials. It theoretically and numerically shows the efficiency of
the method to deal with Coulomb type singularities. Numerical tests to Kohn-
Sham Hamiltonians are in preparation to assess the method on models closer to
practice. This work also opens the way of an analysis of the PAW method for 3D
models. The PAW and VPAW methods are closely related. Since the VPAW method
gives the same eigenvalues as the original Hamiltonian, using the connection
between the PAW and VPAW methods should give an estimation of the error
introduced by the PAW method. This approach has already proven to be
successful in the analysis of this error for a one-dimensional model [11].

##  5 Proof of Theorem 2.1

The general idea of the proof is to isolate the main convergence difficulty
which is the cusp of the eigenfunction $\psi$ located at each nucleus and see
how the VPAW method reduces it. Although the VPAW method reduces the cusps of
the VPAW eigenfunction $\widetilde{\psi}$, it introduces a derivative jump on
the spheres $\partial B(\mathbf{R}_{I},r_{c})$ that blows up as the cut-off
radius shrinks. As in [3], we split the pseudo wave function
$\widetilde{\psi}$ into three parts using the singular expansion (1.19). Let
$\omega$ be a smooth nonnegative cut-off function such that
$\omega(r)=g(\frac{r}{r_{c}})$ where $g$ satisfies

  * •

$g$ is equal to 1 in $B(0,1/4)$,

  * •

$\mathrm{supp}(g)\subset B(0,1/2)$.

Therefore, $\omega$ satisfies $\mathrm{supp}(\omega)\subset B(0,r_{c}/2)$,
$\mathrm{supp}(1-\omega)^{c}\subset B(0,r_{c}/4)$ and
$\|\omega^{(k)}\|_{L^{\infty}(0,r_{c})}\leq\frac{C}{{r_{c}}^{k}}$. Let
$\eta\in{\mathcal{K}}^{\infty,\frac{5}{2}+N-\varepsilon}(\Gamma)$ be the
remainder of the singular expansion (1.19) applied to $\psi$:

$\eta(\mathbf{r})=\psi(\mathbf{r})-\sum\limits_{I=1}^{N_{\mathrm{at}}}\omega(|\mathbf{r}-\mathbf{R}_{I}|)\sum\limits_{j=0}^{N}c_{j}^{I}(\widehat{\mathbf{r-R}_{I}})|\mathbf{r-R}_{I}|^{j},$

where by Theorem 1.6 $c^{I}_{j}=\sum\limits_{|m|\leq\ell\leq j}\psi^{I}_{j\ell
m}Y_{\ell m}$. By Equation (1.16), we have:

$\displaystyle\widetilde{\psi}(\mathbf{r})$
$\displaystyle=\psi(\mathbf{r})-\sum\limits_{I=1}^{N_{\mathrm{at}}}\sum\limits_{k=1}^{N_{\mathrm{paw}}}\langle\widetilde{p}_{k}^{I}\,,\widetilde{\psi}\rangle(\phi^{I}_{k}(\mathbf{r-R}_{I})-\widetilde{\phi}^{I}_{k}(\mathbf{r-R}_{I})),$
$\displaystyle=\sum\limits_{I=1}^{N_{\mathrm{at}}}\omega(|\mathbf{r}-\mathbf{R}_{I}|)\sum\limits_{j=0}^{N}c_{j}^{I}(\widehat{\mathbf{r-R}_{I}})|\mathbf{r-R}_{I}|^{j}+\eta(\mathbf{r})-\sum\limits_{I=1}^{N_{\mathrm{at}}}\sum\limits_{k=1}^{N_{\mathrm{paw}}}\langle\widetilde{p}_{k}^{I}\,,\widetilde{\psi}\rangle(\phi^{I}_{k}(\mathbf{r-R}_{I})-\widetilde{\phi}^{I}_{k}(\mathbf{r-R}_{I})),$
$\displaystyle=\sum\limits_{I=1}^{N_{\mathrm{at}}}\omega(|\mathbf{r}-\mathbf{R}_{I}|)\left(\sum\limits_{j=0}^{N}c_{j}^{I}(\widehat{\mathbf{r-R}_{I}})|\mathbf{r-R}_{I}|^{j}-\sum\limits_{k=1}^{N_{\mathrm{paw}}}\langle\widetilde{p}_{k}^{I}\,,\widetilde{\psi}\rangle(\phi^{I}_{k}(\mathbf{r-R}_{I})-\widetilde{\phi}^{I}_{k}(\mathbf{r-R}_{I}))\right)$
$\displaystyle\qquad\qquad+\sum\limits_{I=1}^{N_{\mathrm{at}}}(1-\omega(|\mathbf{r}-\mathbf{R}_{I}|)\sum\limits_{k=1}^{N_{\mathrm{paw}}}\langle\widetilde{p}_{k}^{I}\,,\widetilde{\psi}\rangle(\phi^{I}_{k}(\mathbf{r-R}_{I})-\widetilde{\phi}^{I}_{k}(\mathbf{r-R}_{I}))+\eta(\mathbf{r}).$
(5.1)

The first part corresponds to the cusp of the pseudo wave function in a
neighborhood of a nucleus. The second part is the $d$-th derivative jump
caused by the lack of regularity at the sphere. The last part is the remainder
appearing in the singular expansion of the original wave function $\psi$. In
this section, we analyze the decay of the Fourier coefficients of the three
parts separately.

We study the plane-wave convergence of the error on the eigenvalue for the
VPAW method where $\ell_{\mathrm{max}}=0$. In this setting, for our purpose, a
singular expansion for $N=1$ is sufficient. Since only PAW functions for the
angular momentum $\ell=0$, $m=0$ are considered, the $00$ index in the PAW
functions is dropped.

######  Proposition 5.1.

Let $(c^{I}_{j})_{0\leq j\leq 1}$ be the functions of the singular expansion
(1.19) of $\psi$ . Let $n\geq 1$ be the number of PAW functions associated to
the angular momentum $\ell=0,m=0$. Then there exists a positive constant $C$
independent of $r_{c}$ and $K$ such that for all $I=1,\dots,N_{\mathrm{at}}$,
for all $\varepsilon>0$ and for $K$ sufficiently large,

$\forall\,0<{r_{c}}<r_{\mathrm{min}},\
\left|\int_{\Gamma}\omega(r)\left(\sum\limits_{j=0}^{1}c^{I}_{j}(\hat{\mathbf{r}})r^{j}-\sum\limits_{k=1}^{n}\langle\widetilde{p}_{k}^{I}\,,\widetilde{\psi}\rangle(\phi^{I}_{k}(\mathbf{r})-\widetilde{\phi}^{I}_{k}(\mathbf{r}))\right)e^{-i\mathbf{K}\cdot\mathbf{r}}\,\mathrm{d}\mathbf{r}\right|\leq\frac{C{r_{c}}^{\min(2n,5)-\varepsilon}}{K^{4}}.$
(5.2)

This proposition states that the cusp of the VPAW eigenfunction
$\widetilde{\psi}$ is reduced by a factor of order
${r_{c}}^{\min(2n,5)-\varepsilon}$. The proof can be found in Section 5.1.

######  Proposition 5.2.

Let $n\geq 1$ be the number of PAW functions associated to the angular
momentum $\ell=0,m=0$. There exists a positive constant $C$ independent of
$r_{c}$ such that for all $I=1,\dots,N_{\mathrm{at}}$ and for $K$ sufficiently
large we have

$\forall\,0<{r_{c}}<r_{\mathrm{min}},\
\left|\int_{\Gamma}(1-\omega(r))\sum\limits_{k=1}^{n}\langle\widetilde{p}_{k}^{I}\,,\widetilde{\psi}\rangle(\phi^{I}_{k}(\mathbf{r})-\widetilde{\phi}^{I}_{k}(\mathbf{r}))e^{-i\mathbf{K}\cdot\mathbf{r}}\,\mathrm{d}\mathbf{r}\right|\leq\frac{C}{{r_{c}}^{d-1}K^{d+2}}.$
(5.3)

By reducing the cusp at a nucleus, the VPAW method introduces a derivative
jump but for a higher order derivative. The proof can be found in Section 5.2.

######  Proposition 5.3.

Let $\eta$ be the remainder of the expansion (5.1) for $N=1$. Let $\eta_{M}$
be the truncation to the wavenumber $M$ of the plane-wave expansion of $\eta$.
Then for all $\varepsilon>0$, we have

$\|\eta_{M}-\eta\|_{H^{1}_{\mathrm{per}}}\leq\frac{1}{M^{5/2-\varepsilon}}\|\eta\|_{H^{7/2-\varepsilon}_{\mathrm{per}}}.$
(5.4)

This is a direct consequence of the regularity of the remainder of the
expansion (1.19) given by Theorem 1.5 (see Remark 1.4).

###  5.1 Proof of Proposition 5.1

Without loss of generality, we assume that there is a nucleus located at
$\mathbf{R}_{I}=0$. Since the estimate in Proposition 5.1 and 5.2 do not
depend on the atomic site, we will drop the upper index $I$.

The following notation is introduced

$\displaystyle p({r})$
$\displaystyle:=(p_{1}({r}),\dots,p_{n}({r}))^{T}\in\mathbb{R}^{n},$
$\displaystyle\widetilde{p}(\mathbf{r})$
$\displaystyle:=(\widetilde{p}_{1}(\mathbf{r}),\dots,\widetilde{p}_{n}(\mathbf{r}))^{T}\in\mathbb{R}^{n},$
$\displaystyle\langle\widetilde{p}\,,f\rangle$
$\displaystyle:=\left(\langle\widetilde{p}_{1}\,,f\rangle,\dots,\langle\widetilde{p}_{n}\,,f\rangle\right)^{T}\in\mathbb{R}^{n},\forall\,f\in
L^{2}_{\mathrm{per}}(\Gamma),$ $\displaystyle\Phi(\mathbf{r})$
$\displaystyle:=(\phi_{1}(\mathbf{r}),\dots,\phi_{n}(\mathbf{r}))^{T}\in\mathbb{R}^{n},$
$\displaystyle\widetilde{\Phi}(\mathbf{r})$
$\displaystyle:=(\widetilde{\phi}_{1}(\mathbf{r}),\dots,\widetilde{\phi}_{n}(\mathbf{r}))^{T}\in\mathbb{R}^{n},$
$\displaystyle{\mathcal{R}}(r)$
$\displaystyle:=(R_{1}(r),\dots,R_{n}(r))^{T}\in\mathbb{R}^{n},$
$\displaystyle\widetilde{{\mathcal{R}}}(r)$
$\displaystyle:=(\widetilde{R}_{1}(r),\dots,\widetilde{R}_{n}(r))^{T}\in\mathbb{R}^{n}.$

For a function $f\in L^{2}([-\frac{1}{2},\frac{1}{2}]^{3})$, we denote by
$f_{\ell m}$ the averaged function

$f_{\ell m}(r)=\int_{S(0,1)}f(\mathbf{r})Y_{\ell
m}(\hat{\mathbf{r}})\,\mathrm{d}\hat{\mathbf{r}}.$ (5.5)

We recall the following identity that will be extensively used in the rest of
the paper:

$e^{-i\mathbf{K}\cdot\mathbf{r}}=4\pi\sum\limits_{|m|\leq\ell}i^{\ell}Y_{\ell
m}(-\hat{\mathbf{K}})Y_{\ell m}(\hat{\mathbf{r}})j_{\ell}(Kr),$ (5.6)

where $j_{\ell}$ is the spherical Bessel function of the first kind.

To prove Proposition 5.1, we start with a lemma that caracterizes the main
difficulties in the plane-wave convergence of the molecular wave function.

######  Lemma 5.4.

Let $\ell$ and $j$ be integers such that $\ell\leq j$. Let $K>0$. Then if
$j+\ell$ is even, asymptotically as $K\to\infty$, we have for any positive
integer $n\geq j+3$,

$\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+2}j_{\ell}(Kr)\,\mathrm{d}r=\frac{\beta_{j,\ell}}{K^{j+3}}+o\left(\frac{1}{K^{n}}\right),$
(5.7)

where

$\beta_{j,\ell}=(-1)^{(j+\ell)/2}(j-\ell+1)!\prod_{k=0}^{\ell}(j-\ell+1+2k),$

and if $j+\ell$ is odd, for any positive integer $n\geq j+3$,

$\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+2}j_{\ell}(Kr)\,\mathrm{d}r=o\left(\frac{1}{K^{n}}\right).$
(5.8)

###### Proof.

We prove the lemma by induction on $\ell$. Let $a_{j+2,\ell}$ be defined by

$a_{j+2,\ell}=\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+2}j_{\ell}(Kr)\,\mathrm{d}r.$
(5.9)

#### Initialization

For $\ell=0$, we have

$j_{0}(x)=\frac{\sin(x)}{x},$

hence for any $j\in\mathbb{N}$,

$a_{j+2,0}=\frac{1}{K}\mathrm{Im}\left(\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+1}e^{iKr}\,\mathrm{d}r\right).$
(5.10)

By integration by parts, we have

$\displaystyle\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+1}e^{iKr}\,\mathrm{d}r$
$\displaystyle=\frac{1}{iK}\underbrace{\left[-\omega(r)r^{j+1}e^{iKr}\right]_{0}^{\frac{r_{c}}{2}}}_{=0}-\frac{1}{iK}\int_{0}^{\frac{r_{c}}{2}}(\omega(r)r^{j+1})^{\prime}e^{iKr}\,\mathrm{d}r$
$\displaystyle=-\frac{1}{iK}\int_{0}^{\frac{r_{c}}{2}}\omega^{\prime}(r)r^{j+1}e^{iKr}\,\mathrm{d}r-\frac{j+1}{iK}\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j}e^{iKr}\,\mathrm{d}r.$

The function $r\mapsto r^{j+1}\omega^{\prime}(r)$ belongs to
$C^{\infty}_{c}(0,\frac{r_{c}}{2})$ hence we have for any $n>j+2$

$\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+1}e^{iKr}\,\mathrm{d}r=-\frac{(j+1)}{iK}\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j}e^{iKr}\,\mathrm{d}r+o\left(\frac{1}{K^{n}}\right)$
(5.11)

By integrating by parts $j$ times and noticing that the functions $r\mapsto
r^{k}\omega^{\prime}(r)$, $k\in\mathbb{N}$ are in
$C^{\infty}_{c}(0,\frac{r_{c}}{2})$, we obtain

$\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+1}e^{iKr}\,\mathrm{d}r=(-1)^{j+1}\frac{(j+1)!}{(iK)^{j+2}}+o\left(\frac{1}{K^{n}}\right).$
(5.12)

Hence using (5.10), if $j$ is even, $a_{j+2,0}=o\left(\frac{1}{K^{n}}\right)$
for all positive integer $n$, otherwise
$a_{j+2,0}=(-1)^{j/2}\frac{(j+1)!}{K^{j+3}}$.

#### Iteration

Using the recurrence relation

$j_{\ell+1}(x)=-j_{\ell}^{\prime}(x)+\frac{\ell j_{\ell}(x)}{x},$ (5.13)

we have the following recurrence relation on $(a_{j,\ell})$:

$\displaystyle a_{j+2,\ell+1}$
$\displaystyle=\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+2}j_{\ell+1}(Kr)\,\mathrm{d}r$
(5.14)
$\displaystyle=\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+2}\left(-j_{\ell}^{\prime}(Kr)+\frac{\ell
j_{\ell}(Kr)}{Kr}\right)\,\mathrm{d}r$ (5.15)
$\displaystyle=\frac{\ell}{K}\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+1}j_{\ell}(Kr)\,\mathrm{d}r-\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+2}j_{\ell}^{\prime}(Kr)\,\mathrm{d}r.$
(5.16)

By integration by parts, using that $r\mapsto r^{j+2}\omega^{\prime}(r)\in
C^{\infty}_{c}(0,\frac{r_{c}}{2})$, we have

$\displaystyle\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+2}j_{\ell}^{\prime}(Kr)\,\mathrm{d}r$
$\displaystyle=-\frac{j+2}{K}\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+1}j_{\ell}(Kr)\,\mathrm{d}r+o\left(\frac{1}{K^{n}}\right).$

Thus, we have by iteration,

$\displaystyle a_{j+2,\ell+1}$
$\displaystyle=\frac{j+\ell+2}{K}\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+1}j_{\ell}(Kr)\,\mathrm{d}r+o\left(\frac{1}{K^{n}}\right)$
$\displaystyle=\frac{j+\ell+2}{K}a_{j+1,\ell}+o\left(\frac{1}{K^{n}}\right)$
$\displaystyle=\frac{\prod_{k=0}^{\ell+1}(j-\ell+2k)}{K^{\ell+1}}a_{j-\ell+1,0}+o\left(\frac{1}{K^{n}}\right).$

Hence if $j+\ell+1$ is odd, $a_{j+2,\ell+1}=o\left(\frac{1}{K^{n}}\right)$ for
all positive integer $n$, else if $j+\ell+1$ is even, for any $n\geq j+3$,

$a_{j+2,\ell+1}=(-1)^{(j-\ell+1)/2}(j-\ell)!\frac{\prod_{k=0}^{\ell+1}(j-\ell+2k)}{K^{j+3}}+o\left(\frac{1}{K^{n}}\right).$

∎

######  Lemma 5.5.

Let $N\in\mathbb{N}^{*}$. Let $c_{j}$ be the functions of the singular
expansion (1.19) at 0 of $\psi$. Let $\psi_{j\ell m}$ be the coefficients such
that

$c_{j}(\hat{\mathbf{r}})=\sum_{\ell=0}^{j}\sum_{|m|\leq\ell}\psi_{j\ell
m}Y_{\ell m}(\hat{\mathbf{r}}).$

Then, we have asymptotically as $K$ goes to $\infty$ and for any positive
integer $n$,

$\int_{[-\frac{1}{2},\frac{1}{2}]^{3}}\omega(r)\sum\limits_{j=0}^{N}c_{j}(\hat{\mathbf{r}})r^{j}e^{-i\mathbf{K}\cdot\mathbf{r}}\,\mathrm{d}\mathbf{r}=4\pi\sum\limits_{j=0}^{N}\sum_{\begin{subarray}{c}\ell=0\\\
j+\ell\ \textrm{odd}\end{subarray}}^{j}\sum_{|m|\leq\ell}i^{\ell}Y_{\ell
m}(-\hat{\mathbf{K}})\frac{\beta_{j,\ell}\psi_{j\ell
m}}{K^{j+3}}+o\left(\frac{1}{K^{n}}\right).$ (5.17)

###### Proof.

We have since $\mathrm{supp}(\omega)\subset(0,\frac{r_{c}}{2})$

$\displaystyle\int_{[-\frac{1}{2},\frac{1}{2}]^{3}}\omega(r)\sum\limits_{j=0}^{N}c_{j}(\hat{\mathbf{r}})r^{j}e^{-i\mathbf{K}\cdot\mathbf{r}}\,\mathrm{d}\mathbf{r}$
$\displaystyle=\int_{0}^{\frac{r_{c}}{2}}r^{2}\omega(r)\sum\limits_{j=0}^{N}r^{j}\sum_{\ell=0}^{j}\sum_{|m|\leq\ell}\psi_{j\ell
m}\int_{S(0,1)}Y_{\ell
m}(\hat{\mathbf{r}})e^{-i\mathbf{K}\cdot\mathbf{r}}\,\mathrm{d}\hat{\mathbf{r}}\,\mathrm{d}r.$
(5.18)

Using the scattering expansion (5.6) and applying Lemma 5.4, we get

$\displaystyle\int_{[-\frac{1}{2},\frac{1}{2}]^{3}}\omega(r)\sum\limits_{j=0}^{N}c_{j}(\hat{\mathbf{r}})r^{j}e^{-i\mathbf{K}\cdot\mathbf{r}}\,\mathrm{d}\mathbf{r}$
$\displaystyle=4\pi\sum\limits_{j=0}^{N}\sum_{\ell=0}^{j}\sum_{|m|\leq\ell}i^{\ell}Y_{\ell
m}^{*}(-\hat{\mathbf{K}})\int_{0}^{\frac{r_{c}}{2}}\omega(r)r^{j+2}\psi_{j\ell
m}j_{\ell}(Kr)\,\mathrm{d}r$ (5.19)
$\displaystyle=4\pi\sum\limits_{j=0}^{N}\sum_{\begin{subarray}{c}\ell=0\\\
j+\ell\ \textrm{odd}\end{subarray}}^{j}\sum_{|m|\leq\ell}i^{\ell}Y_{\ell
m}(-\hat{\mathbf{K}})\frac{\beta_{j,\ell}\psi_{j\ell
m}}{K^{j+3}}+o\left(\frac{1}{K^{n}}\right).$ (5.20)

∎

According to Lemma 5.5, the slowest decaying term is the term associated to
$j=1$ and $\ell=0$. Proposition 5.1 is then a direct consequence of Lemmas 5.5
and 5.6 stated below.

######  Lemma 5.6.

Let $\psi$ be an eigenfunction of $H\psi=E\psi$, with H defined in (1.3). Let
$n$ be the number of PAW functions associated to the angular momentum
$\ell=0,m=0$.

Then there exist coefficients $(\alpha_{k})_{1\leq k\leq n}$ and a positive
constant $C$ independent of $r_{c}$ such that

$\|\psi_{00}-\alpha^{T}{\mathcal{R}}\|_{L^{2}(0,{r_{c}})}\leq
C{r_{c}}^{1/2-\varepsilon+\min(2n,5)},$

where $\psi_{00}$ denotes the averaged function $\psi$ according to (5.5).
Moreover for these coefficients, we have

$\psi_{00}(0)-\alpha^{T}{\mathcal{R}}(0)=0.$

We first prove that Proposition 5.1 follows from Lemmas 5.5 and 5.6.

###### Proof of Proposition 5.1.

By definition of the PAW atomic functions $\phi_{k}=R_{k}Y_{00}$, these
functions have a cusp at 0 and satisfy $R_{k}^{\prime}(0)=-ZR_{k}(0)$. By
Theorem 1.5, we have $\psi_{100}=-Z\psi(0)$. Using that the PAW pseudo wave
functions are smooth in a neighbourdhood of 0 and Lemma 5.5 for $N=1$, we
obtain

$\left|\int_{\Gamma}\omega(r)\left(\sum\limits_{j=0}^{1}c^{I}_{j}(\hat{\mathbf{r}})r^{j}-\sum\limits_{k=1}^{n}\langle\widetilde{p}_{k}^{I}\,,\widetilde{\psi}\rangle(\phi^{I}_{k}(\mathbf{r})-\widetilde{\phi}^{I}_{k}(\mathbf{r}))\right)e^{-i\mathbf{K}\cdot\mathbf{r}}\,\mathrm{d}\mathbf{r}\right|\leq\frac{C}{K^{4}}|\psi(0)-\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\mathcal{R}(0)|,$

where $C>0$ is a constant independent of $r_{c}$ and $K$. We will now show
that
$|\psi(0)-\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\mathcal{R}(0)|\leq
C{r_{c}}^{\min(2n,5)-\varepsilon}$, for some constant $C>0$ independent of
$r_{c}$. By duality of $\widetilde{\phi}_{j}$ and $\widetilde{p}_{k}$, for any
$\alpha\in\mathbb{R}^{n}$ we have

$\psi(0)-\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\mathcal{R}(0)=\psi_{00}(0)-\alpha^{T}{\mathcal{R}}(0)-\langle\widetilde{p}\,,\widetilde{\psi}-\alpha^{T}\widetilde{\Phi}\rangle^{T}{\mathcal{R}}(0).$
(5.21)

First we rewrite
$\langle\widetilde{p}\,,\widetilde{\psi}-\alpha^{T}\widetilde{\Phi}\rangle$ in
a more convenient way. In a neighbordhood of 0, we have

$\psi-\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\Phi=\widetilde{\psi}-\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\widetilde{\Phi}.$

By multiplying by $\widetilde{p}_{k}$, $k=1,..,n$ and integrating over the
ball $B(0,{r_{c}})$, we obtain

$\langle\widetilde{p}_{k}\,,\psi\rangle-\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\langle\widetilde{p}_{k}\,,\Phi\rangle=0,$

so

$\langle\widetilde{p}\,,\widetilde{\psi}\rangle=A^{-1}\langle\widetilde{p}\,,\psi\rangle,$

where $A=(\langle\widetilde{p}_{j}\,,\phi_{k}\rangle)_{1\leq j,k\leq
n}=(\langle p_{j}\,,R_{k}\rangle_{[0,r_{c}]})_{1\leq j,k\leq n}$ which is
invertible by Assumption 1. By definition of the PAW functions, the
3D-integrals can be reduced to integrals on an interval

$\langle\widetilde{p}\,,\widetilde{\psi}\rangle=A^{-1}\langle
p\,,\psi_{00}\rangle_{[0,{r_{c}}]},$

where

$\langle f\,,g\rangle_{[0,r_{c}]}=\int_{0}^{r_{c}}f(r)g(r)r^{2}\,\mathrm{d}r.$

By duality of the PAW functions,
$\langle\widetilde{p}_{j}\,,\widetilde{\phi}_{k}\rangle=\delta_{jk}$, we have
$\langle\widetilde{p}\,,\alpha^{T}\widetilde{\Phi}\rangle=A^{-1}\langle
p\,,\alpha^{T}{\mathcal{R}}\rangle_{[0,r_{c}]}$. Hence,

$\langle\widetilde{p}\,,\widetilde{\psi}-\alpha^{T}\widetilde{\Phi}\rangle=A^{-1}\langle
p\,,\psi_{00}-\alpha^{T}{\mathcal{R}}\rangle_{[0,{r_{c}}]}.$

By Lemma A.7, there exists a constant $C$ independent of ${r_{c}}$ such that
for any ${r_{c}}>0$,

$\left|\langle
p\,,\psi_{00}-\alpha^{T}{\mathcal{R}}\rangle_{[0,r_{c}]}^{T}A^{-T}{\mathcal{R}}^{\prime}(0)\right|\leq\frac{C}{{r_{c}}^{3/2}}\|\psi_{00}-\alpha^{T}{\mathcal{R}}\|_{L^{2}(B(0,{r_{c}}))}\leq\frac{C}{{r_{c}}^{1/2}}\|\psi_{00}-\alpha^{T}{\mathcal{R}}\|_{L^{2}(0,{r_{c}})}.$
(5.22)

From Lemma 5.6, we know that there exists $\alpha\in\mathbb{R}^{n}$ such that

$\|\psi_{00}-\alpha^{T}{\mathcal{R}}\|_{L^{2}(0,{r_{c}})}\leq
C{r_{c}}^{\frac{1}{2}-\varepsilon+\min(2n,5)}\quad\text{ and
}\quad\alpha^{T}{\mathcal{R}}(0)=\psi(0).$

Inserting this equation into (5.21) finishes the proof.

∎

To show Proposition 5.1, it suffices to prove Lemma 5.6. In order to show this
lemma, we start with a few intermediary results and introduce some notation.
Let $\uppsi_{k}\in\mathbb{R}$ and $\zeta_{k}\in\mathbb{R}^{n}$ be,
respectively, the coefficients of the singular expansion of $\psi_{00}$ and
${\mathcal{R}}$:

$\displaystyle\psi_{00}(r)$
$\displaystyle=\sum\limits_{j=0}^{N}\uppsi_{j}r^{j}+\eta_{N+1}(r),\quad\eta_{N+1}\in{\mathcal{K}}^{\infty,\frac{5}{2}+N-\varepsilon}(\Gamma)$
(5.23) $\displaystyle{\mathcal{R}}(r)$
$\displaystyle=\sum\limits_{j=0}^{N}\zeta_{j}r^{j}+\xi_{N+1}(r),\quad\xi_{N+1}\in{\mathcal{K}}^{\infty,\frac{5}{2}+N-\varepsilon}(\Gamma).$
(5.24)

The potential $V^{I}$ defined in Equation (1.20) is smooth, hence we can find
$(v_{2k})_{0\leq k\leq N}\in\mathbb{R}^{N}$ such that

$\displaystyle
V_{00}^{I}(r)=\sum\limits_{k=0}^{N}v_{2k}r^{2k}+{\mathcal{O}}(r^{2N+2}).$
(5.25)

The atomic potential $W_{\mathrm{at}}$ in Equation 1.6 is also smooth so there
is $(w_{2k})_{0\leq k\leq N}\in\mathbb{R}^{N}$ such that

$(W_{\mathrm{at}})_{00}(r)=\sum\limits_{k=0}^{N}w_{2k}r^{2k}+{\mathcal{O}}(r^{2N+2}).$
(5.26)

Finally, we denote ${\mathcal{E}}$ the diagonal matrix with entries
$(\epsilon_{1},\dots,\epsilon_{n})$.

######  Lemma 5.7.

Let ${\mathcal{R}}=(R_{1},\dots,R_{n})^{T}$ where $R_{k}$ is defined in (1.7).

There exists $(\mu_{j}^{(k)})_{0\leq j\leq k\leq n-1}$ and
$(\nu_{j}^{(k)})_{0\leq j\leq k\leq n-1}$ such that

$\displaystyle\zeta_{2k}$
$\displaystyle=\sum\limits_{j=0}^{k}\mu_{j}^{(k)}{\mathcal{E}}^{j}\zeta_{0}$
(5.27) $\displaystyle\zeta_{2k+1}$
$\displaystyle=\sum\limits_{j=0}^{k}\nu_{j}^{(k)}{\mathcal{E}}^{j}\zeta_{0},$
(5.28)

with $\mu_{k}^{(k)}\not=0$ for any $0\leq k\leq n-1$.

Moreover, the vectors $(\zeta_{k})_{1\leq k\leq 5}$ satisfies

$\displaystyle\zeta_{1}$ $\displaystyle=-Z\zeta_{0}$ $\displaystyle\zeta_{2}$
$\displaystyle=-\frac{1}{3}{\mathcal{E}}\zeta_{0}+\frac{1}{3}(Z^{2}+w_{0})\zeta_{0}$
$\displaystyle\zeta_{3}$
$\displaystyle=\frac{2}{9}Z{\mathcal{E}}\zeta_{0}-\left(\frac{Z^{3}}{18}+\frac{2}{9}Zw_{0}\right)\zeta_{0}$
$\displaystyle\zeta_{4}$
$\displaystyle=\frac{1}{30}{\mathcal{E}}^{2}\zeta_{0}-\left(\frac{Z^{2}}{18}+\frac{w_{0}}{15}\right){\mathcal{E}}\zeta_{0}+\left(\frac{Z^{4}}{180}+\frac{Z^{2}w_{0}}{18}+\frac{w_{0}^{2}}{30}+\frac{w_{2}}{10}\right)\zeta_{0}$
$\displaystyle\zeta_{5}$
$\displaystyle=-\frac{23}{1350}Z{\mathcal{E}}^{2}\zeta_{0}+\left(\frac{Z^{3}}{135}+\frac{23}{675}w_{0}Z\right){\mathcal{E}}\zeta_{0}-\left(\frac{Z^{5}}{2700}+\frac{Z^{3}w_{0}}{135}+\frac{23}{1350}Zw_{0}^{2}+\frac{11}{150}Zw_{2}\right)\zeta_{0}.$

We give an explicit formula of $\zeta_{k}$ for $k\leq 5$ because we will show
that we cannot systematically improve the bound in Lemma 5.6 for $n\geq 2$.

###### Proof.

This lemma is proved by iteration using Equation (1.22) applied to each
$R_{k}$, $k=1,\dots,n$

$\forall N\geq 1,\
\frac{(N+1)(N+2)}{2}\zeta_{N+1}=-Z\zeta_{N}+(w*\zeta)_{N-1}-{\mathcal{E}}\zeta_{N-1}.$

Here $(w_{2k})_{k\geq 0}$ and $(\zeta_{2k})_{k\geq 0}$ are defined
respectively by (5.24) and (5.26). Since $W_{\mathrm{at}}$ is smooth, we have

$(w*\zeta)_{2k}=\sum\limits_{j=0}^{k}w_{2k-2j}\zeta_{2j},\quad(w*\zeta)_{2k+1}=\sum\limits_{j=0}^{k}w_{2k-2j}\zeta_{2j+1}.$

By iteration, one notices that $\mu_{k}^{(k)}\not=0$ for all $k\in\mathbb{N}$.
∎

###### Proof of Lemma 5.6.

To minimize $\|\psi_{00}-\alpha^{T}{\mathcal{R}}\|_{L^{2}(0,r_{c})}$ with
respect to $\alpha$ for $r_{c}$ small, we need to determine how many
successive terms in the singular expansion of $\psi_{00}$ can be canceled with
$n$ functions $R_{k},k=1,\dots,n$, _i.e._ we need to determine for which
$N_{\mathrm{max}}\leq 2n$ we have

$\forall\,0\leq k\leq N_{\mathrm{max}},\ \alpha^{T}\zeta_{k}=\uppsi_{k}.$
(5.29)

Expressing $\zeta_{k}$ (resp. $\uppsi_{k}$) as a linear combination of
$(\zeta_{0},\dots,{\mathcal{E}}^{\lceil\frac{N_{\mathrm{max}}}{2}\rceil-1}\zeta_{0})$
(resp.
$(\uppsi_{0},\dots,E^{\lceil\frac{N_{\mathrm{max}}}{2}\rceil-1}\uppsi_{0})$)
using Lemma 5.7, the linear system (5.29) can be reformulated as

$\begin{pmatrix}M_{1}\\\ M_{2}\end{pmatrix}\begin{pmatrix}\zeta_{0}^{T}\\\
\zeta_{0}^{T}{\mathcal{E}}\\\ \vdots\\\ \zeta_{0}^{T}{\mathcal{E}}^{\lceil
N_{\mathrm{max}}/2\rceil-1}\end{pmatrix}\alpha=\begin{pmatrix}N_{1}\\\
N_{2}\end{pmatrix}\begin{pmatrix}\uppsi_{0}\\\ E\uppsi_{0}\\\ \vdots\\\
E^{\lceil N_{\mathrm{max}}/2\rceil-1}\uppsi_{0}\end{pmatrix},$ (5.30)

where $M_{1}=(\mu_{j}^{(k)})_{0\leq j,k\leq\lceil
N_{\mathrm{max}}/2\rceil-1},$ and $M_{2}=(\nu_{j}^{(k)})_{0\leq j,k\leq\lceil
N_{\mathrm{max}}/2\rceil-1},$ with $\mu_{j}^{(k)}$ and $\nu_{j}^{(k)}$ given
by Lemma 5.7. $N_{1}$ and $N_{2}$ the same matrices as $M_{1}$ and $M_{2}$ but
where the coefficients $\mu_{j}^{(k)}$ and $\nu_{j}^{(k)}$ are generated using
$(v_{2k})$ instead of $(w_{2k})$.

We will show that if $n\leq 2$, then $N_{\mathrm{max}}=2n$, otherwise,
$N_{\mathrm{max}}=5$. Equation (5.30) is equivalent to

$\begin{cases}M_{1}(\zeta_{0}^{T}{\mathcal{E}}^{k})_{0\leq k\leq\lceil
N_{\mathrm{max}}/2\rceil-1}\alpha=N_{1}(E^{k}\uppsi_{0})_{0\leq k\leq\lceil
N_{\mathrm{max}}/2\rceil-1}\\\ M_{2}(\zeta_{0}^{T}{\mathcal{E}}^{k})_{0\leq
k\leq\lceil N_{\mathrm{max}}/2\rceil-1}\alpha=N_{2}(E^{k}\uppsi_{0})_{0\leq
k\leq\lceil N_{\mathrm{max}}/2\rceil-1},\end{cases}$

hence Equation (5.30) has a solution if and only if

$M_{2}M_{1}^{-1}N_{1}(E^{k}\uppsi_{0})_{0\leq k\leq\lceil
N_{\mathrm{max}}/2\rceil-1}=N_{2}(E^{k}\uppsi_{0})_{0\leq k\leq\lceil
N_{\mathrm{max}}/2\rceil-1}.$

Since this holds for any value $E$, a necessary and sufficient condition to
solve (5.30) is $M_{2}M_{1}^{-1}N_{1}=N_{2}$.

For $n=1$, $M_{1}=N_{1}=\begin{pmatrix}1\end{pmatrix}$ and
$M_{2}=N_{2}=\begin{pmatrix}-Z\end{pmatrix}$ hence (5.30) is solvable when
$N_{\mathrm{max}}=2$. The remainder $\psi_{00}-\alpha{\mathcal{R}}$ belongs to
${\mathcal{K}}^{\infty,\frac{7}{2}-\varepsilon}(\Gamma)$, hence by Lemma A.1,
there exists a constant $C$ independent of $r_{c}$ such that

$\|\psi_{00}-\alpha{\mathcal{R}}\|_{L^{2}(0,r_{c})}\leq
C{r_{c}}^{5/2-\varepsilon}.$

For $n=2$, we can check using Lemma 5.7 that

$M_{2}M_{1}^{-1}=\begin{pmatrix}-Z&0\\\
\frac{Z^{3}}{6}&-\frac{2Z}{3}\end{pmatrix}.$

Thus $M_{2}M_{1}^{-1}N_{1}=N_{2}$ and (5.30) has a solution
$(\alpha_{1},\alpha_{2})$ such that $\psi_{00}(r)-\alpha^{T}{\mathcal{R}}(r)$
belongs to ${\mathcal{K}}^{\infty,\frac{11}{2}-\varepsilon}(\Gamma)$, hence by
Lemma A.1, $\|\psi_{00}(r)-\alpha^{T}{\mathcal{R}}(r)\|_{L^{2}(0,r_{c})}\leq
C{r_{c}}^{9/2-\varepsilon}$.

For $n\geq 3$, the dependence on the atomic potential $W_{\mathrm{at}}$ does
not vanish in $M_{2}M^{-1}_{1}$. For example, the $(3,1)$ coordinate of the
matrix $M_{2}M^{-1}_{1}$ has a term equal to $\frac{391}{6075}Z^{3}w_{0}$
which is unlikely to be compensated in general. For $n\geq 3$, we thus have:

$\psi_{00}(r)-\alpha^{T}{\mathcal{R}}(r)={\mathcal{O}}(r^{5}),$

hence $\|\psi_{00}(r)-\alpha^{T}{\mathcal{R}}(r)\|_{L^{2}(0,r_{c})}\leq
C{r_{c}}^{11/2}$ for a constant $C$ independent of $r_{c}$. ∎

###  5.2 Proof of Proposition 5.2

The derivative jump of the $d$-th derivative jump at $r=r_{c}$ of
$\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\widetilde{{\mathcal{R}}}^{(d)}$
needs to be estimated.

######  Lemma 5.8.

There exists a positive constant $C$ independent of $r_{c}$ such that

$\forall\,0<{r_{c}}<r_{\mathrm{min}},\
\left|\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}[\widetilde{{\mathcal{R}}}^{(d)}]_{{r_{c}}}\right|\leq\frac{C}{{r_{c}}^{d-1}},$

where $[\widetilde{{\mathcal{R}}}^{(d)}]_{{r_{c}}}=\lim\limits_{h\to
0,h>0}\widetilde{{\mathcal{R}}}^{(d)}(r_{c}+h)-\widetilde{{\mathcal{R}}}^{(d)}(r_{c}-h)$.
For any $k\in\mathbb{N}$, we have

$\left\|\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\left({\mathcal{R}}^{(k)}-\widetilde{{\mathcal{R}}}^{(k)}\right)\right\|_{L^{\infty}(0,r_{c})}\leq\frac{C}{{r_{c}}^{k}}.$

The proof of this lemma is given in the appendix.

###### Proof of Proposition 5.2.

Using the relation (5.6), we have

$\int_{[-\frac{1}{2},\frac{1}{2}]^{3}}(1-\omega(r))\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}(\Phi-\widetilde{\Phi})e^{-i\mathbf{K}\cdot\mathbf{r}}\,\mathrm{d}\mathbf{r}=4\pi\int_{0}^{{r_{c}}}(1-\omega(r))\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}({\mathcal{R}}(r)-\widetilde{{\mathcal{R}}}(r))j_{0}(Kr)r^{2}\,\mathrm{d}r.$
(5.31)

Since $\omega$ is equal to 1 in a neighbourhood of 0, we can restrict the
integral in the equation above to the interval $({r_{c}}-\eta,{r_{c}})$ for
some $\eta>0$. Recall that

$j_{0}(x)=\frac{\sin(x)}{x},$

thus

$\int_{[-\frac{1}{2},\frac{1}{2}]^{3}}(1-\omega(r))\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}(\Phi-\widetilde{\Phi})e^{-i\mathbf{K}\cdot\mathbf{r}}\,\mathrm{d}\mathbf{r}=\frac{4\pi}{K}\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\int_{{r_{c}}-\eta}^{{r_{c}}}(1-\omega(r))({\mathcal{R}}(r)-\widetilde{{\mathcal{R}}}(r))\sin(Kr)r\,\mathrm{d}r.$
(5.32)

We denote by $f$ the function $r\mapsto
r(1-\omega(r))({\mathcal{R}}(r)-\widetilde{{\mathcal{R}}}(r))$ and use

$\int_{{r_{c}}-\eta}^{{r_{c}}}f(r)\sin(Kr)\,\mathrm{d}r=\mathrm{Im}\left(\int_{{r_{c}}-\eta}^{{r_{c}}}f(r)e^{iKr}\,\mathrm{d}r\right).$

By definition of the cut-off function, for any $k\in\mathbb{N}$, we have

$f^{(k)}({r_{c}}-\eta)=0,$ (5.33)

and for $k\in\mathbb{N}^{*}$, $(1-\omega)^{(k)}({r_{c}})=0$. Thus by
integration by parts,

$\int_{{r_{c}}-\eta}^{{r_{c}}}f(r)e^{iKr}\,\mathrm{d}r=\left[f(r)\frac{e^{iKr}}{iK}\right]_{{r_{c}}-\eta}^{r_{c}}-\frac{1}{iK}\int_{{r_{c}}-\eta}^{{r_{c}}}f^{\prime}(r)e^{iKr}\,\mathrm{d}r=\frac{i}{K}\int_{{r_{c}}-\eta}^{{r_{c}}}f^{\prime}(r)e^{iKr}\,\mathrm{d}r.$

As ${\mathcal{R}}-\widetilde{{\mathcal{R}}}$ is $C^{d-1}$ but not $C^{d}$ at
${r_{c}}$, by integrating by parts $d$ times, we have:

$\int_{{r_{c}}-\eta}^{{r_{c}}}f(r)e^{iKr}\,\mathrm{d}r=\frac{i^{d+1}{r_{c}}}{K^{d+1}}[\widetilde{{\mathcal{R}}}^{(d)}]_{{r_{c}}}e^{iK{r_{c}}}-\frac{i^{d+1}}{K^{d+1}}\int_{{r_{c}}-\eta}^{{r_{c}}}f^{(d+1)}(r)e^{iKr}\,\mathrm{d}r.$
(5.34)

Thus inserting the last equation in (5.32), we obtain

$\int_{[-\frac{1}{2},\frac{1}{2}]^{3}}(1-\omega(r))\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}(\Phi-\widetilde{\Phi})e^{-i\mathbf{K}\cdot\mathbf{r}}\,\mathrm{d}\mathbf{r}=\mathrm{Im}(i^{d+1}e^{iK{r_{c}}})\frac{4\pi{r_{c}}}{K^{d+2}}\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}[\widetilde{{\mathcal{R}}}^{(d)}]_{{r_{c}}}\\\
-\mathrm{Im}\left(\frac{4\pi
i^{d+1}}{K^{d+2}}\int_{{r_{c}}-\eta}^{{r_{c}}}\langle\widetilde{p}\,,f\rangle^{T}f^{(d+1)}(r)e^{iKr}\,\mathrm{d}r\right).$

According to Lemma 5.8,

$\left|\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}[\widetilde{{\mathcal{R}}}^{(d)}]_{{r_{c}}}\right|\leq\frac{C}{{r_{c}}^{d-1}}.$

Furthermore, we have

$\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}f^{(d+1)}(r)=r\left((1-\omega)\langle\widetilde{p}\,,\widetilde{\psi}\rangle({\mathcal{R}}-\widetilde{{\mathcal{R}}})\right)^{(d+1)}+\left((1-\omega)\langle\widetilde{p}\,,\widetilde{\psi}\rangle({\mathcal{R}}-\widetilde{{\mathcal{R}}})\right)^{(d)}.$

By assumption on $\omega$, we know that
$\|\omega^{(k)}\|_{L^{\infty}(0,r_{c})}\leq C{r_{c}}^{-k}$, hence

$\displaystyle\left\|\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}f^{(d+1)}\right\|_{L^{\infty}(0,r_{c})}$
$\displaystyle\leq
r_{c}\sum\limits_{k=0}^{d+1}\binom{d+1}{k}\left\|(1-\omega)^{(k)}({\mathcal{R}}-\widetilde{{\mathcal{R}}})^{(d+1-k)}\right\|_{L^{\infty}(0,r_{c})}$
$\displaystyle\qquad\qquad+\sum\limits_{k=0}^{d}\binom{d}{k}\left\|(1-\omega)^{(k)}({\mathcal{R}}-\widetilde{{\mathcal{R}}})^{(d-k)}\right\|_{L^{\infty}(0,r_{c})}$
$\displaystyle\leq\frac{C}{{r_{c}}^{d}},$

where we used Lemma 5.8. Thus, we obtain

$\displaystyle\left|\int_{{r_{c}}-\eta}^{{r_{c}}}\langle\widetilde{p}\,,f\rangle^{T}f^{(d+1)}(r)e^{iKr}\,\mathrm{d}r\right|$
$\displaystyle\leq
r_{c}\left\|\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}f^{(d+1)}\right\|_{L^{\infty}(0,r_{c})}$
$\displaystyle\leq\frac{C}{{r_{c}}^{d-1}},$

which finishes the proof of this proposition. ∎

###  5.3 Convergence theorem

To prove the estimate on the eigenvalues, we will use the following classical
result ([33], p. 68).

######  Proposition 5.9.

Let $H$ be a self-adjoint coercive $H^{1}$-bounded operator,
$E_{1}\leq\dots\leq E_{n}$ be the lowest eigenvalues of $H$ and
$\psi_{1},\dots,\psi_{n}$ be $L^{2}$-normalized associated eigenfunctions. Let
$E^{(M)}_{1}\leq\dots\leq E^{(M)}_{n}$ be the lowest eigenvalues of the
Rayleigh quotient of $H$ restricted to the subspace $V_{M}$ of dimension $M$.

Let $w_{k}\in V_{M}$ for $1\leq k\leq n$ be such that

$\sum\limits_{k=1}^{n}\|w_{k}-\psi_{k}\|^{2}_{H^{1}}<1$

Then there exists a positive constant $C$ which depends on the $H^{1}$ norm of
$H$ and the coercivity constant such that for all $1\leq k\leq n$

$\left|E^{(M)}_{k}-E_{k}\right|\leq
C\sum\limits_{k=1}^{n}\|w_{k}-\psi_{k}\|^{2}_{H^{1}}$

We would like to apply this result to
$\psi_{K}=(\mathrm{Id}+T)\widetilde{\psi}_{K}$ where $\widetilde{\psi}_{K}$ is
the truncation of the plane-wave expansion of $f$ to the wave number ${K}$. In
order to do this, we need to show that the $H^{1}_{\mathrm{per}}$-norm of
$(\mathrm{Id}+T)\widetilde{\psi}$ is bounded by a the
$H^{1}_{\mathrm{per}}$-norm of $\widetilde{\psi}$ independently of the cut-off
radius $r_{c}$.

######  Lemma 5.10.

There exists a positive constant $C$ independent of $r_{c}$ such that for any
function $f\in H^{1}_{\mathrm{per}}([-\frac{1}{2},\frac{1}{2}]^{3})$

$\left\|(\mathrm{Id}+T)f\right\|_{H^{1}_{\mathrm{per}}}\leq
C\|f\|_{H^{1}_{\mathrm{per}}}.$

###### Proof.

By definition, we have
$(\mathrm{Id}+T)f=f+\langle\widetilde{p}\,,f\rangle^{T}(\Phi-\widetilde{\Phi})$.
By Lemma A.6, we have

$\langle\widetilde{p}\,,f\rangle^{T}({\mathcal{R}}(r)-\widetilde{{\mathcal{R}}}(r))=C_{r_{c}}^{T}(C_{r_{c}}G(P)C_{r_{c}}^{T})^{-1}C_{r_{c}}\int_{0}^{1}\chi(t)P(t)f_{00}(r_{c}t)t^{2}\,\mathrm{d}t\cdot\left(\begin{pmatrix}C_{1}^{-1}\\\
0\end{pmatrix}{\mathcal{R}}(r)-P(\tfrac{r}{r_{c}})\right),$

with $\left\|\begin{pmatrix}C_{1}^{-1}\\\
0\end{pmatrix}{\mathcal{R}}(r)\right\|_{L^{\infty}(0,r_{c})}$ uniformly
bounded with respect to $r_{c}$. Thus it suffices to bound
$\int_{0}^{1}t^{2}\chi(r_{c}t)f_{00}(r_{c}t)P(t)\,\mathrm{d}t$:

$\displaystyle\left|\int_{0}^{1}t^{2}\chi(r_{c}t)f_{00}(r_{c}t)P(t)^{2}\,\mathrm{d}t\right|$
$\displaystyle\leq\left(\int_{0}^{1}t^{2}\chi(r_{c}t)^{2}P(t)\,\mathrm{d}t\right)^{1/2}\left(\int_{0}^{1}t^{2}f_{00}(r_{c}t)^{2}\,\mathrm{d}t\right)^{1/2}$
$\displaystyle\leq\frac{C}{{r_{c}}^{3/2}}\left(\int_{0}^{r_{c}}r^{2}f_{00}(r)^{2}\,\mathrm{d}r\right)^{1/2}$
$\displaystyle\leq\frac{C}{{r_{c}}^{3/2}}\left(\int_{0}^{r_{c}}r^{2}f_{00}(r)^{6}\,\mathrm{d}r\right)^{1/6}\left(\int_{0}^{r_{c}}r^{2}\,\mathrm{d}r\right)^{1/3}$
$\displaystyle\leq\frac{C}{{r_{c}}^{1/2}}\|f\|_{L^{6}(B_{r_{c}})}$
$\displaystyle\leq\frac{C}{{r_{c}}^{1/2}}\|f\|_{H^{1}(B_{r_{c}})}.$

We obtain

$\left|\langle
p\,,f_{00}\rangle^{T}_{[0,r_{c}]}({\mathcal{R}}(r)-\widetilde{{\mathcal{R}}}(r))\right|\leq\frac{C}{{r_{c}}^{1/2}}\|f\|_{H^{1}(B_{r_{c}})}.$

Since $\left\|\begin{pmatrix}C_{1}^{-1}\\\
0\end{pmatrix}{\mathcal{R}}^{\prime}(r)\right\|_{L^{\infty}(0,r_{c})}={\mathcal{O}}\left(\frac{1}{r_{c}}\right)$,
we can prove similarly

$\left|\langle
p\,,f_{00}\rangle^{T}_{[0,r_{c}]}({\mathcal{R}}^{\prime}(r)-\widetilde{{\mathcal{R}}}^{\prime}(r))\right|\leq\frac{C}{{r_{c}}^{3/2}}\|f\|_{H^{1}(B_{r_{c}})}.$

Thus we get

$\left\|(\mathrm{Id}+T)f\right\|_{H^{1}_{\mathrm{per}}}\leq
C\|f\|_{H^{1}_{\mathrm{per}}}.$

∎

We have all the elements to prove the main convergence theorem.

###### Proof of Theorem 2.1.

Using Equation (5.1), Propositions 5.1, 5.2 and 5.3, we have

$\displaystyle\|\widetilde{\psi}_{M}-\widetilde{\psi}\|_{H^{1}_{\mathrm{per}}}$
$\displaystyle\leq\left\|\widetilde{\psi}_{M}-\widetilde{\psi}-\eta_{M}+\eta\right\|_{H^{1}_{\mathrm{per}}}+\left\|\eta_{M}-\eta\right\|_{H^{1}_{\mathrm{per}}}$
$\displaystyle\leq\left(\sum\limits_{K\geq
M,\mathbf{K}\in(2\pi\mathbb{Z})^{3}}(1+K^{2})\left(\frac{{r_{c}}^{\min(n,3)+1-\varepsilon}}{K^{4}}+\frac{1}{{r_{c}}^{d-1}K^{d+2}}\right)^{2}\right)^{\frac{1}{2}}+o\left(\frac{1}{M^{5/2-\varepsilon}}\right).$

By Proposition 5.9, we obtain

$\displaystyle|E_{M}-E|$ $\displaystyle\leq
C\|(\mathrm{Id}+T)(\widetilde{\psi}_{M}-\widetilde{\psi})\|_{H^{1}_{\mathrm{per}}}^{2}$
$\displaystyle\leq
C\|\widetilde{\psi}_{M}-\widetilde{\psi}\|_{H^{1}_{\mathrm{per}}}^{2}$
$\displaystyle\leq C\left(\left(\sum\limits_{K\geq
M,\mathbf{K}\in(2\pi\mathbb{Z})^{3}}(1+K^{2})\left(\frac{{{r_{c}}}^{\min(2n,5)-\varepsilon}}{K^{4}}+\frac{1}{{{r_{c}}}^{d-1}K^{d+2}}\right)^{2}\right)^{1/2}+o\left(\frac{1}{M^{5/2-\varepsilon}}\right)\right)^{2}$
$\displaystyle\leq
C\left(\frac{{{r_{c}}}^{2\min(2n,5)-2\varepsilon}}{M^{3}}+\frac{{{r_{c}}}^{\min(2n,5)-\varepsilon}}{M^{4-\varepsilon}}+\frac{1}{{{r_{c}}}^{2d-2}M^{2d-1}}+o\left(\frac{1}{M^{5-\varepsilon}}\right)\right).$

∎

##  6 Acknowledgements

The author is grateful to Xavier Blanc, Antoine Levitt, Eric Cancès and Marc
Torrent for helpful discussion and fruitful comments.

##  Appendix A Appendix

We have gathered in this section proofs of some technical lemmas, most of
which are simple transpositions of lemmas that can be found in [11, 3].

###  A.1 Results related to the weighted Sobolev space
${\mathcal{K}}^{\infty,a}(\Gamma)$

######  Lemma A.1.

Let $f\in{\mathcal{K}}^{\infty,a}(\Gamma)$ and $0<R<1$. Let
$\ell\in\mathbb{N}$ and $m\in\mathbb{N}$ such that $|m|\leq\ell$. Then there
exists a constant $C$ independent of $f$ such that

$\int_{0}^{R}|f_{\ell m}(r)|r^{2}\,\mathrm{d}r\leq
CR^{a+\frac{3}{2}}\|f\|_{{\mathcal{K}}^{\infty,a}},$

and for $a\geq 1$

$\int_{0}^{R}|f_{\ell m}(r)|^{2}\,\mathrm{d}r\leq
CR^{2a-2}\|f\|_{{\mathcal{K}}^{\infty,a}}.$

###### Proof.

Since $Y_{\ell m}\in L^{\infty}(S(0,1))$, we have

$\displaystyle\int_{0}^{R}|f_{\ell m}(r)|r^{2}\,\mathrm{d}r$
$\displaystyle\leq C\int_{B(0,R)}|f(\mathbf{r})|\,\mathrm{d}\mathbf{r}$
$\displaystyle\leq
C\int_{B(0,R)}r^{a}r^{-a}|f(\mathbf{r})|\,\mathrm{d}\mathbf{r}$
$\displaystyle\leq
C\left(\int_{0}^{R}r^{2a+2}\,\mathrm{d}r\right)^{1/2}\left(\int_{B(0,R)}r^{-2a}|f(\mathbf{r})|^{2}\,\mathrm{d}\mathbf{r}\right)^{1/2}$
$\displaystyle\leq CR^{a+\frac{3}{2}}\|f\|_{{\mathcal{K}}^{\infty,a}},$

where in the fourth inequality we used the definition of the weighted Sobolev
space ${\mathcal{K}}^{\infty,a}$. The second identity is proved the same way.
∎

######  Lemma A.2.

Let $N\in\mathbb{N}^{*}$ and $\eta$ a radial function such that
$\eta\in{\mathcal{K}}^{\infty,5/2+N-\varepsilon}(\Gamma)$, $\varepsilon>0$.
Then for $R$ sufficiently small, we have

$\|\eta\|_{L^{\infty}(0,R)}\leq\|f\|_{{\mathcal{K}}^{\infty,a}}R^{N+\frac{1}{2}-\varepsilon}.$

###### Proof.

By definition of the weighted Sobolev space, we have for $R$ sufficiently
small,

$\int_{B(0,R)}|\eta(r)|^{2}r^{-5-2N+2\varepsilon}\,\mathrm{d}\mathbf{r}<\infty,$

hence

$\int_{0}^{R}|\eta(r)|^{2}r^{-3-2N+2\varepsilon}\,\mathrm{d}r<\infty,$

Similarly we have

$\int_{0}^{R}|\eta^{\prime}(r)|^{2}r^{-1-2N+2\varepsilon}\,\mathrm{d}r<\infty.$

Therefore,

$\begin{cases}\int_{0}^{R}|\eta(r)|^{2}\,\mathrm{d}r\leq
R^{2N+3-2\varepsilon}\int_{0}^{R}|\eta(r)|^{2}r^{-1-2N+2\varepsilon}\,\mathrm{d}r,\\\
\int_{0}^{R}|\eta^{\prime}(r)|^{2}\,\mathrm{d}r\leq
R^{2N+1-2\varepsilon}\int_{0}^{R}|\eta^{\prime}(r)|^{2}r^{-1-2N+2\varepsilon}\,\mathrm{d}r.\end{cases}$

By the Sobolev embedding theorem, we have the result. ∎

######  Remark A.3.

Lemma A.2 implies that the remainder $\eta_{N}$ of the singularity expansion
(1.19) of radial functions are bounded :
$\|\eta_{N}\|_{L^{\infty}(0,r_{c})}\leq C{r_{c}}^{N+1}$, where the constant is
independent of $r_{c}$.

###  A.2 Validity of the assumptions for the hydrogenoid atom

We show in this subsection that the Assumptions 1 and 3 hold in the particular
case of the hydrogenoid atom, _i.e._ where in (1.6) $W_{\mathrm{at}}=0$. The
eigenfunctions of the hydrogenoid atom can be written

$\varphi_{n\ell m}(\mathbf{r})=R_{n\ell}(r)Y_{\ell m}(\hat{\mathbf{r}}),\quad
n\geq\ell+1,$

with

$R_{n\ell}(r)=\sqrt{{\left(\frac{2Z}{n}\right)}^{3}\frac{(n-\ell-1)!}{2n(n+\ell)!}}e^{-Zr/{n}}\left(\frac{2Zr}{n}\right)^{\ell}L_{n-\ell-1}^{(2\ell+1)}\left(\frac{2Zr}{n}\right),$
(A.1)

where $L_{n-\ell-1}^{(2\ell+1)}$ denotes the generalized Laguerre polynomials.
The constant term of $L_{n-\ell-1}^{(2\ell+1)}$ is equal to
$\binom{n+\ell}{n-\ell-1}>0$, hence Assumption 3 holds for the hydrogenoid
atom wave function.

######  Lemma A.4.

Let $R_{n\ell}(r)=L_{n-1}(\frac{2Z}{n})e^{-\frac{Zr}{n}}$, where
$\mathrm{deg}(L_{n-1})=n-1$. Let $r>0$ and
${\mathcal{R}}_{k}=(R_{1\ell}^{(k)}(r),\dots,R_{n\ell}^{(k)}(r))^{T}$ for
$0\leq k\leq n-1$. Then the matrix
$({\mathcal{R}}_{0},\dots,{\mathcal{R}}_{n-1})$ is invertible.

###### Proof.

We have

$R_{n\ell}^{(k)}(r)=\left(\frac{Z}{n}\right)^{k}e^{-\frac{Zr}{n}}\sum\limits_{j=0}^{k}2^{j}\tbinom{k}{j}L_{n-1}^{(j)}\left(\frac{2Zr}{n}\right).$

Let $\mathcal{P}=\left(L_{j}^{(k)}(\frac{2Zr}{j+1})\right)_{0\leq j,k\leq
n-1}$, $\mathcal{M}=(2^{j}\binom{k}{j})_{0\leq j,k\leq n-1}$ and
$\mathcal{Z}=\left((\frac{Z}{k+1})^{j}\right)_{0\leq j,k\leq n-1}$. One can
check that

$({\mathcal{R}}_{0},\dots,{\mathcal{R}}_{n-1})=\mathcal{Z}\
\mathrm{diag}(e^{-Zr},\dots,e^{-\frac{Zr}{n}})\mathcal{P}\mathcal{M}.$

$\mathcal{P}$ and $\mathcal{M}$ are both triangular with no null entry on the
diagonal. $\mathcal{Z}$ is a Vandermonde matrix, hence
$({\mathcal{R}}_{0},\dots,{\mathcal{R}}_{n-1})$ is invertible. ∎

###  A.3 Lemmas related to PAW functions

Let $P_{k}$, $k\in\mathbb{N}$ the polynomials defined by

$P_{k}(t)=\frac{1}{2^{k}k!}(t^{2}-1)^{k}.$ (A.2)

By definition, these polynomials form a basis of even polynomials and satisfy

$\begin{cases}P_{k}^{(j)}(1)=0,\quad 0\leq j\leq k-1,\\\
P_{k}^{(k)}(1)=1.\end{cases}$

Let $P$ be the vector $(P_{0},\dots,P_{d-1})^{T}$ and
$C_{r_{c}}\in\mathbb{R}^{n\times d}$ be the matrix such that

$\widetilde{{\mathcal{R}}}(t)=C_{r_{c}}P(\tfrac{t}{r_{c}}).$ (A.3)

The following lemma summarizes the main properties of the matrix $C_{r_{c}}$.

######  Lemma A.5.

Let $C_{1}\in\mathbb{R}^{n\times n}$ and $C_{2}\in R^{n\times(n-d)}$ be the
matrices such that

$\widetilde{\mathcal{R}}(r)=\Big{(}C_{1}\ \Big{|}\
C_{2}\Big{)}P(\tfrac{r}{r_{c}}).$

Moreover, $C_{1}$ is invertible, the norm of $C_{1}^{-1}C_{2}$ is uniformly
bounded with respect to $r_{c}$ and

$C_{2}^{T}C_{1}^{-T}e_{0}={\mathcal{O}}(r_{c}).$

###### Proof.

Let $C_{r_{c}}$ be the matrix $\Big{(}C_{1}\ \Big{|}\ C_{2}\Big{)}$. Let
$c_{j}$ be the columns of $C_{r_{c}}$. By continuity of
$\widetilde{{\mathcal{R}}}$ and of its derivatives at $r_{c}$, and by our
choice of the polynomials $P_{k}$, the columns of $C_{r_{c}}$ satisfy

$\forall\,0\leq j\leq d-1,\
c_{j}={r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c})-\sum\limits_{k=0}^{j-1}P_{k}^{(j)}c_{k}.$
(A.4)

Hence $c_{k}$ is a linear combination of the vectors
${r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c})$ for $j\leq k$ with coefficients that
are _independent_ of $r_{c}$. Moreover we can deduce that the transformation
of $(c_{j})_{0\leq j\leq n-1}$ to
$({r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c}))_{0\leq j\leq n-1}$ is invertible. If
$r_{c}$ is sufficiently small, by Assumption 1, the family
$({r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c}))_{0\leq j\leq n-1}$ is linearly
independent, thus we can define $(g_{j})_{0\leq j\leq n-1}$ to be the dual
family to $(c_{j})_{0\leq j\leq n-1}$ (_i.e._ $c_{j}^{T}g_{k}=\delta_{jk}$)
and we have $\|g_{j}\|={\mathcal{O}}\left(\tfrac{1}{{r_{c}}^{n-1}}\right)$.
Hence, using the recurrence (A.4), we can show that the norm of
$C^{-1}_{1}C_{2}$ is uniformly bounded with respect to $r_{c}$.

To prove $C_{2}^{T}C_{1}^{-T}e_{0}={\mathcal{O}}(r_{c})$, first notice that
$C_{1}^{-T}e_{0}=g_{1}$. Since $P_{0}$ is a constant polynomial, for $j\geq
1$, we have

$c_{j}={r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c})-\sum\limits_{k=1}^{j-1}P_{k}^{(j)}c_{k}.$

Thus for $j=n$, we have

$c_{n}^{T}g_{0}={r_{c}}^{n}{\mathcal{R}}^{(n)}(r_{c})^{T}g_{0}-\sum\limits_{k=1}^{j-1}P_{k}^{(j)}c_{k}^{T}g_{0}={\mathcal{O}}(r_{c}),$

and by iteration, we can check that

$\forall\,n\leq j\leq d-1,\ c_{j}^{T}g_{0}={\mathcal{O}}(r_{c}).$

∎

######  Lemma A.6.

We have

$\langle\widetilde{p}\,,f\rangle^{T}({\mathcal{R}}(r)-\widetilde{{\mathcal{R}}}(r))=C_{r_{c}}^{T}(C_{r_{c}}G(P)C_{r_{c}}^{T})^{-1}C_{r_{c}}\int_{0}^{1}\chi(t)P(t)f_{00}(r_{c}t)t^{2}\,\mathrm{d}t\cdot\left(\begin{pmatrix}C_{1}^{-1}\\\
0\end{pmatrix}{\mathcal{R}}(r)-P(\tfrac{r}{r_{c}})\right),$

with $P$ being the vector of the polynomials $P_{k}$ defined in (A.2),
$C_{r_{c}}$ the matrix of coefficients of $\widetilde{{\mathcal{R}}}$ in the
basis $(P_{k})$ given in (A.3) and $G(P)$ the matrix

$G(P)=\int_{0}^{1}\chi(t)P(t)P(t)^{T}t^{2}\,\mathrm{d}t.$

The norm of the matrix
$C_{r_{c}}^{T}(C_{r_{c}}G(P)C_{r_{c}}^{T})^{-1}C_{r_{c}}$ is uniformly bounded
as $r_{c}$ goes to 0 and we have

$\left\|\begin{pmatrix}C_{1}^{-1}\\\
0\end{pmatrix}{\mathcal{R}}(r)\right\|_{L^{\infty}(0,r_{c})}\leq C\quad\text{
and }\quad\left\|\begin{pmatrix}C_{1}^{-1}\\\
0\end{pmatrix}{\mathcal{R}}^{\prime}(r)\right\|_{L^{\infty}(0,r_{c})}\leq\frac{C}{r_{c}}.$

###### Proof.

Since $G(P)$ is positive-definite, $G(P)^{1/2}$ exists. Note that the matrix
$G(P)^{1/2}C_{r_{c}}^{T}(C_{r_{c}}G(P)C_{r_{c}}^{T})^{-1}C_{r_{c}}G(P)^{1/2}$
is symmetric and is a projector, hence its norm is independent of $r_{c}$.

Writing down the Taylor expansion of ${\mathcal{R}}$ at $r_{c}$, we obtain

$\displaystyle{\mathcal{R}}(r)$
$\displaystyle=\sum\limits_{k=0}^{n-1}\frac{(r-r_{c})^{k}}{k!}{\mathcal{R}}^{(k)}(r_{c})+{\mathcal{O}}((r-r_{c})^{n})$
$\displaystyle=\sum\limits_{k=0}^{n-1}\frac{1}{k!}\left(\frac{r}{r_{c}}-1\right)^{k}{r_{c}}^{k}{\mathcal{R}}^{(k)}(r_{c})+{\mathcal{O}}((r-r_{c})^{n}).$

By Lemma A.5, $\|C_{1}^{-1}{r_{c}}^{k}{\mathcal{R}}^{(k)}(r_{c})\|$ is
uniformly bounded as $r_{c}$ goes to 0, which concludes the proof of the
lemma. ∎

######  Lemma A.7.

Let $f\in L^{2}(\Gamma)$. Let $\ell\geq 0$, $|m|\leq\ell$ be integers. Let $n$
be the number of PAW functions associated to the angular momentum $\ell,m$ for
a cut-off radius $r_{c}$. There exists a constant independent of $r_{c}$ and
$f$ such that

$|\langle\widetilde{p}\,,f\rangle^{T}{\mathcal{R}}(0)|\leq\frac{C}{{r_{c}}^{3/2}}\|f\|_{L^{2}(B_{r_{c}})}.$

###### Proof.

The proof of this lemma is similar to the proof of Proposition 3.2 in [3]. We
give the proof in case $\ell=m=0$. First, it is possible to show that

$\langle\widetilde{p}\,,f\rangle=\int_{0}^{1}\chi(t)(C_{r_{c}}^{(Q)}G_{r_{c}})^{-1}C_{r_{c}}^{(Q)}Q(t)f(r_{c}t)t^{2}\,\mathrm{d}t,$
(A.5)

where $Q(t)=(Q_{0}(t),\dots,Q_{d-1}(t))^{T}$ is a vector of even polynomials
which forms a basis of even polynomials of degree at most $2d-2$,

$\widetilde{{\mathcal{R}}}(x)=C_{r_{c}}^{(Q)}Q(\tfrac{x}{r_{c}}),$

with $C_{r_{c}}^{(Q)}\in\mathbb{R}^{n\times d}$ and

$G_{r_{c}}=\int_{0}^{1}\chi(t)Q(t)\mathcal{R}(r_{c}t)^{T}t^{2}\,\mathrm{d}t\in\mathbb{R}^{d\times
n}.$

By Lemma 5.7, we have that $(\zeta_{2k})_{0\leq k\leq n-1}$ and
$(\zeta_{2k+1})_{0\leq k\leq n-1}$ defined by the singular expansion of
${\mathcal{R}}$ :

$\mathcal{R}(t)=\sum\limits_{k=0}^{n-1}\zeta_{2k}t^{2k}+\zeta_{2k+1}t^{2k+1}+\eta_{2n}(t),$
(A.6)

satisfy

$\zeta_{2k}=\sum\limits_{j=0}^{k}\mu_{j}^{(k)}{\mathcal{E}}^{j}\zeta_{0}\quad\text{
and
}\quad\zeta_{2k+1}=\sum\limits_{j=0}^{k}\nu_{j}^{(k)}{\mathcal{E}}^{j}\zeta_{0},$

where $\mu_{k}^{(k)}\not=0$ and ${\mathcal{E}}$ is the diagonal matrix of the
eigenvalues $(\epsilon_{1},\dots,\epsilon_{n})$. By Assumption 3, $\zeta_{0}$
has no null entry. The eigenvalues of the atomic operator (1.6) for a fixed
$\ell,m$ are simple, hence $({\mathcal{E}}^{j}\zeta_{0})_{0\leq j\leq n-1}$ is
a linearly independent family. Hence, $(\zeta_{2k})_{0\leq k\leq n-1}$ is a
basis of $\mathbb{R}^{n}$. Let $(h_{k})$ be the dual basis to
$({r_{c}}^{2j}{\mathcal{E}}^{j}\zeta_{0})_{0\leq j\leq n-1}$, _i.e._
$h_{k}^{T}{r_{c}}^{2j}{\mathcal{E}}^{j}\zeta_{0}=\delta_{kj}$.

Injecting (A.6) in the definition of $G_{r_{c}}$, we obtain

$\displaystyle G_{r_{c}}^{T}$
$\displaystyle=\int_{0}^{1}\chi(t){\mathcal{R}}(r_{c}t)Q(t)^{T}t^{2}\,\mathrm{d}t$
(A.7)
$\displaystyle=\int_{0}^{1}\chi(t)\left(\sum\limits_{k=0}^{n-1}{r_{c}}^{2k}{\mathcal{E}}^{k}\zeta_{0}\sum\limits_{j=k}^{n-1}\mu_{k}^{(j)}{r_{c}}^{2j-2k}t^{2j}+r_{c}\sum\limits_{k=0}^{n-1}{r_{c}}^{2k}{\mathcal{E}}^{k}\zeta_{0}\sum\limits_{j=k}^{n-1}\nu_{k}^{(j)}{r_{c}}^{2j-2k}t^{2j+1}+\eta_{2n}(r_{c}t)\right)Q(t)^{T}t^{2}\,\mathrm{d}t.$
(A.8)

By Lemma A.1, we have

$\left|\int_{0}^{1}\chi(t)\eta_{2n}(r_{c}t)Q(t)^{T}t^{2}\,\mathrm{d}t\right|=\frac{1}{r_{c}^{3}}\left|\int_{0}^{r_{c}}\chi(\tfrac{t}{r_{c}})\eta_{2n}(t)Q(\tfrac{t}{r_{c}})^{T}t^{2}\,\mathrm{d}t\right|\leq
C{r_{c}}^{2n}.$ (A.9)

Let $(Q_{k})_{0\leq k\leq d-1}$ be the even polynomials such that

$\int_{0}^{1}\chi(t)t^{2j}Q_{k}(t)t^{2}\,\mathrm{d}t=\delta_{jk}.$

Such polynomials exist since the Gram matrix
$(\int_{0}^{1}\chi(t)t^{2j+2k+2}\,\mathrm{d}t)_{0\leq j,k\leq n-1}$ is
invertible. Let

$X_{j}=\int_{0}^{1}\chi(t)t^{2j+1}Q(t)t^{2}\,\mathrm{d}t,$

and

$\mathcal{H}=\begin{pmatrix}h_{0}^{T}\\\ \vdots\\\
h_{n-1}^{T}\end{pmatrix}\in\mathbb{R}^{n\times n}.$ (A.10)

Then denoting by $e_{k}$ the $k$-th canonical vector, we have

$\displaystyle\mathcal{H}G_{r_{c}}^{T}$
$\displaystyle=\sum\limits_{k=0}^{n-1}\mu_{k}^{(k)}e_{k}e_{k}^{T}+\sum\limits_{k=0}^{n-1}\sum\limits_{j=k+1}^{n-1}\mu_{k}^{(j)}{r_{c}}^{2j-2k}e_{k}e_{j}^{T}+r_{c}\sum\limits_{k=0}^{n-1}\nu_{k}^{(j)}{r_{c}}^{2j-2k}e_{k}X_{j}^{T}+{\mathcal{O}}({r_{c}}^{2})$
$\displaystyle=\sum\limits_{k=0}^{n-1}\mu_{k}^{(k)}e_{k}e_{k}^{T}+{\mathcal{O}}({r_{c}}).$

Let

$\mathcal{A}=\sum\limits_{k=0}^{n-1}\mu_{k}^{(k)}e_{k}e_{k}^{T}\in\mathbb{R}^{n\times
n},$ (A.11)

and $\Pi$ the transition matrix such that

$C_{r_{c}}^{(Q)}=C_{r_{c}}\Pi,$

where $C_{r_{c}}$ is defined in (A.3). Hence we have

$\displaystyle(C_{r_{c}}^{(Q)})^{T}(C_{r_{c}}^{(Q)}G_{r_{c}})^{-T}{\mathcal{R}}(0)$
$\displaystyle=(C_{r_{c}}^{(Q)})^{T}(C_{r_{c}}^{(Q)}G_{r_{c}})^{-T}\mathcal{H}^{-1}e_{0}$
$\displaystyle=(C_{r_{c}}^{(Q)})^{T}(C_{r_{c}}^{(Q)}G_{r_{c}}\mathcal{H}^{T})^{-T}e_{0}$
$\displaystyle=\Pi^{T}\begin{pmatrix}I_{n}\\\
M^{T}+{\mathcal{O}}(r_{c})\end{pmatrix}\left(C^{-1}C_{r_{c}}\Pi\left(\begin{pmatrix}\mathcal{A}\\\
0\end{pmatrix}+{\mathcal{O}}(r_{c})\right)\right)^{-T}e_{0}$
$\displaystyle=\Pi^{T}\begin{pmatrix}I_{n}\\\
M^{T}+{\mathcal{O}}(r_{c})\end{pmatrix}\left(\Big{(}\ I_{n}\Big{|}\
M+{\mathcal{O}}(r_{c})\Big{)}\Pi\left(\begin{pmatrix}\mathcal{A}\\\
0\end{pmatrix}+{\mathcal{O}}(r_{c})\right)\right)^{-T}e_{0},$

where we used Lemma A.5 in the third and fourth inequality. Decomposing $\Pi$
into four blocks

$\Pi=\begin{pmatrix}\Pi_{1}&\Pi_{2}\\\
\Pi_{3}&\Pi_{4}\end{pmatrix},\quad\text{with }\Pi_{1}\in\mathbb{R}^{n\times
n},$

we obtain

$\displaystyle(C_{r_{c}}^{(Q)})^{T}(C_{r_{c}}^{(Q)}G_{r_{c}})^{-T}{\mathcal{R}}(0)$
$\displaystyle=\left(\begin{pmatrix}\Pi_{1}^{T}+\Pi_{3}^{T}M^{T}\\\
\Pi_{2}^{T}+\Pi_{4}^{T}M^{T}\end{pmatrix}+{\mathcal{O}}(r_{c})\right)\left(\Pi_{1}\mathcal{A}+M\Pi_{3}\mathcal{A}+{\mathcal{O}}(r_{c})\right)e_{0}$
$\displaystyle=\begin{pmatrix}\mathcal{A}^{-1}\\\
(\Pi_{2}^{T}+\Pi_{4}^{T}M^{T})(\Pi_{1}+M\Pi_{3})^{-1}\end{pmatrix}e_{0}+{\mathcal{O}}(r_{c}).$

Hence
$\|(C_{r_{c}}^{(Q)})^{T}(C_{r_{c}}^{(Q)}G_{r_{c}})^{-T}{\mathcal{R}}(0)\|$ is
uniformly bounded as $r_{c}$ goes to 0. Thus, there exists a constant $C$
independent of $r_{c}$ and $f$ such that :

$\displaystyle\left|\langle\widetilde{p}\,,f\rangle^{T}{\mathcal{R}}(0)\right|$
$\displaystyle\leq C\int_{0}^{1}|f(r_{c}t)|t^{2}\,\mathrm{d}t$
$\displaystyle\leq\frac{C}{{r_{c}}^{3/2}}\|f\|_{L^{2}_{\mathrm{per}}}.$

∎

We can now prove Lemma 5.8.

###### Proof of Lemma 5.8.

We start with the proof of the estimate of
$\left[\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\widetilde{{\mathcal{R}}}^{(d)}\right]_{r_{c}}$.
We have using (A.5) and (A.3)

$\displaystyle\left[\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\widetilde{{\mathcal{R}}}^{(d)}\right]_{r_{c}}$
$\displaystyle=\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}({\mathcal{R}}^{(d)}(r_{c})-\widetilde{{\mathcal{R}}}^{(d)}(r_{c}))$
(A.12)
$\displaystyle=\frac{1}{{r_{c}}^{d}}\int_{0}^{1}\chi(t)\psi_{00}(r_{c}t)Q(t)t^{2}\,\mathrm{d}t\cdot(C^{(Q)}_{r_{c}})^{T}(C^{(Q)}_{r_{c}}G_{r_{c}})^{-T}\left({r_{c}}^{d}{\mathcal{R}}^{(d)}(r_{c})-C_{r_{c}}P^{(d)}(1)\right).$
(A.13)

First, we prove that

$(C^{(Q)}_{r_{c}})^{T}(C^{(Q)}_{r_{c}}G_{r_{c}})^{-T}\left({r_{c}}^{d}{\mathcal{R}}^{(d)}(r_{c})-C_{r_{c}}P^{(d)}(1)\right)=\begin{pmatrix}0\\\
*\end{pmatrix}+{\mathcal{O}}(r_{c}),$ (A.14)

then

$\int_{0}^{1}\chi(t)\psi_{00}(r_{c}t)Q(t)t^{2}\,\mathrm{d}t=\psi(0)e_{0}+{\mathcal{O}}(r_{c}).$
(A.15)

If both statements are true, then we deduce that there exists a constant $C$
independent of $r_{c}$ such that

$\left|\left[\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}\widetilde{{\mathcal{R}}}^{(d)}\right]_{r_{c}}\right|\leq\frac{C}{{r_{c}}^{d-1}}.$

#### Step 1 (proof of (A.14))

By (5.24), we have for $0\leq j\leq 2n-1$ even

$\displaystyle{r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c})$
$\displaystyle=\sum\limits_{k=j/2}^{n-1}\zeta_{2k}\frac{(2k)!}{(2k-j)!}{r_{c}}^{2k}+\zeta_{2k+1}\frac{(2k+1)!}{(2k+1-j)!}{r_{c}}^{2k+1}+{\mathcal{O}}({r_{c}}^{2n})$
$\displaystyle=\sum\limits_{k=j/2}^{n-1}\sum\limits_{\ell=0}^{k}\mu_{\ell}^{(k)}{r_{c}}^{2\ell}{\mathcal{E}}^{\ell}\zeta_{0}\frac{(2k)!}{(2k-j)!}{r_{c}}^{2k-2\ell}+r_{c}\sum\limits_{k=j/2}^{n-1}\sum\limits_{\ell=0}^{k}\nu_{\ell}^{(k)}{r_{c}}^{2\ell}{\mathcal{E}}^{\ell}\zeta_{0}\frac{(2k+1)!}{(2k+1-j)!}{r_{c}}^{2k-2\ell}+{\mathcal{O}}({r_{c}}^{2n}),$

where we applied Lemma A.2 to estimate the remainder of the singular
expansion. By noticing that
$\mathcal{H}{r_{c}}^{2j}{\mathcal{E}}^{j}\zeta_{0}=e_{j}$ with $\mathcal{H}$
defined in (A.10), we have

$(C^{(Q)}_{r_{c}})^{T}(C^{(Q)}_{r_{c}}G_{r_{c}})^{-T}{r_{c}}^{j}{\mathcal{E}}^{j}\zeta_{0}=\begin{pmatrix}\mathcal{A}^{-1}\\\
*\end{pmatrix}e_{j}+{\mathcal{O}}(r_{c}).$

Using
$\|(C^{(Q)}_{r_{c}})(C^{(Q)}_{r_{c}}G_{r_{c}})^{-1}\|={\mathcal{O}}\left(\frac{1}{{r_{c}}^{2n-2}}\right)$,
we thus get

$\displaystyle(C^{(Q)}_{r_{c}})^{T}(C^{(Q)}_{r_{c}}G_{r_{c}})^{-T}{r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c})$
$\displaystyle=\begin{pmatrix}I_{n}\\\
*\end{pmatrix}\sum\limits_{k=j/2}^{n-1}e_{k}+{\mathcal{O}}(r_{c}).$ (A.16)

For $0\leq j\leq 2n-1$ odd, we have

${r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c})={r_{c}}^{j}\zeta_{j}+\sum\limits_{k=\frac{j+1}{2}}^{n-1}\zeta_{2k}\frac{(2k)!}{(2k-j)!}{r_{c}}^{2k}+\zeta_{2k+1}\frac{(2k+1)!}{(2k+1-j)!}{r_{c}}^{2k+1}+{\mathcal{O}}({r_{c}}^{2n}),$

similarly to the even case, we thus obtain

$\displaystyle(C^{(Q)}_{r_{c}})^{T}(C^{(Q)}_{r_{c}}G_{r_{c}})^{-T}{r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c})$
$\displaystyle=\begin{pmatrix}\mathcal{A}^{-1}\\\
*\end{pmatrix}\sum\limits_{k=\frac{j+1}{2}}^{n-1}\nu_{k}^{(k)}e_{k}+{\mathcal{O}}(r_{c}).$
(A.17)

For $j\geq 2n$, using
$\|(C^{(Q)}_{r_{c}})(C^{(Q)}_{r_{c}}G_{r_{c}})^{-1}\|={\mathcal{O}}\left(\frac{1}{{r_{c}}^{2n-2}}\right)$,
then

$(C^{(Q)}_{r_{c}})^{T}(C^{(Q)}_{r_{c}}G_{r_{c}})^{-T}{r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c})={\mathcal{O}}(r_{c}).$
(A.18)

From (A.16) (when $d\leq 2n-1$ and $d$ is even), (A.17) (when $d\leq 2n-1$ and
$d$ is odd) or (A.18) (when $d\geq 2n$), we have

$(C^{(Q)}_{r_{c}})^{T}(C^{(Q)}_{r_{c}}G_{r_{c}})^{-T}{r_{c}}^{d}{\mathcal{R}}^{(d)}(r_{c})=\begin{pmatrix}0\\\
*\end{pmatrix}+{\mathcal{O}}(r_{c}).$

It remains to prove the same statement for the other part. By definition of
the polynomials $P_{k}$ (A.2), we have

$P^{(d)}(1)=(\underbrace{0,\dots,0}_{\lfloor\frac{d}{2}\rfloor},*,\dots,*)^{T},$

so $C_{r_{c}}P^{(d)}(1)$ is a linear combination of the last
$\lceil\frac{d}{2}\rceil$ columns of $C_{r_{c}}$. However, by Lemma A.5, we
know that except the first column of $C_{r_{c}}$, the columns of $C_{r_{c}}$
do not depend on ${\mathcal{R}}(r_{c})$ and by (A.16) and (A.17), for $j\geq
1$,

$e_{0}^{T}(C^{(Q)}_{r_{c}})^{T}(C^{(Q)}_{r_{c}}G_{r_{c}})^{-T}{r_{c}}^{j}{\mathcal{R}}^{(j)}(r_{c})={\mathcal{O}}(r_{c}),$

which finishes the proof of the first step.

#### Step 2 (proof of (A.15))

Since $\psi\in H^{2}_{\mathrm{per}}(\Gamma)$, by Sobolev embedding theorem,
$\psi$ is continuous, hence $\psi(0)$ is finite. Thus

$\displaystyle\int_{0}^{1}\chi(t)\psi_{00}(r_{c}t)Q(t)t^{2}\,\mathrm{d}t$
$\displaystyle=\int_{0}^{1}\chi(t)\left(\psi(0)+\int_{0}^{r_{c}t}\psi_{00}^{\prime}(u)\,\mathrm{d}u\right)Q(t)t^{2}\,\mathrm{d}t$
$\displaystyle=\psi(0)e_{0}+\int_{0}^{1}\chi(t)\int_{0}^{r_{c}t}\psi_{00}^{\prime}(u)\,\mathrm{d}u\,Q(t)t^{2}\,\mathrm{d}t,$

by definition of the polynomials $Q_{k}$.

We have ($C$ denotes a constant independent of $r_{c}$)

$\displaystyle\left|\int_{0}^{1}\chi(t)\int_{0}^{r_{c}t}\psi_{00}^{\prime}(u)\,\mathrm{d}u\
Q(t)t^{2}\,\mathrm{d}t\right|$
$\displaystyle\leq\left(\int_{0}^{1}\chi(t)^{2}Q(t)^{2}t^{4}\,\mathrm{d}t\int_{0}^{1}\left(\int_{0}^{r_{c}t}\psi_{00}^{\prime}(u)\,\mathrm{d}u\right)^{2}\,\mathrm{d}t\right)^{1/2}$
$\displaystyle\leq
Cr_{c}\left(\int_{0}^{1}\left(\frac{1}{u}\psi_{00}^{\prime}\right)^{2}u^{2}\,\mathrm{d}u\right)^{1/2}.$

Using Hardy inequality, we get

$\displaystyle\left|\int_{0}^{1}\chi(t)\int_{0}^{r_{c}t}\psi_{00}^{\prime}(u)\,\mathrm{d}u\
Q(t)t^{2}\,\mathrm{d}t\right|$ $\displaystyle\leq
Cr_{c}\|\psi\|_{H^{2}_{\mathrm{per}}},$

which ends the proof of (A.15).

The proof of the bound on
$\langle\widetilde{p}\,,\widetilde{\psi}\rangle^{T}({\mathcal{R}}^{(k)}-\widetilde{{\mathcal{R}}}^{(k)})$
is a direct extension of the proof of (A.14). ∎

## References

  * [1] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing, SIAM review, 59 (2017), pp. 65–98.
  * [2] X. Blanc, E. Cancès, and M.-S. Dupuy, Variational projector augmented-wave method, Comptes Rendus Mathematique, 355 (2017), pp. 665 – 670\.
  * [3] X. Blanc, E. Cancès, and M.-S. Dupuy, Variational projector augmented-wave method: theoretical analysis and preliminary numerical results, Numer. Math., 144 (2020), pp. 271–321.
  * [4] P. E. Blochl, Projector augmented-wave method, Phys. Rev. B, 50 (1994), pp. 17953–17979.
  * [5] E. Cancès and N. Mourad, A mathematical perspective on density functional perturbation theory, Nonlinearity, 27 (2014), pp. 1999–2033.
  * [6] E. Cances and N. Mourad, Existence of a type of optimal norm-conserving pseudopotentials for Kohn–Sham models, Communications in Mathematical Sciences, 14 (2016), pp. 1315–1352.
  * [7] I. Catto, C. Le Bris, and P.-L. Lions, On the thermodynamic limit for Hartree-Fock type models, vol. 18, 2001, pp. 687–760.
  * [8] , On some periodic Hartree-type models for crystals, vol. 19, 2002, pp. 143–190.
  * [9] H. Chen and R. Schneider, Numerical analysis of augmented plane wave methods for full-potential electronic structure calculations, ESAIM Math. Model. Numer. Anal., 49 (2015), pp. 755–785.
  * [10] M.-S. Dupuy, Analysis of the projector augmented-wave method for electronic structure calculations in periodic settings, theses, Université Sorbonne Paris Cité, Sept. 2018.
  * [11] , Projector augmented-wave method: an analysis in a one-dimensional setting, ESAIM Math. Model. Numer. Anal., 54 (2020), pp. 25–58.
  * [12] M. S. P. Eastham, The spectral theory of periodic differential equations, Texts in Mathematics (Edinburgh), Scottish Academic Press, Edinburgh; Hafner Press, New York, 1973.
  * [13] Y. V. Egorov and B.-W. Schulze, Pseudo-differential operators, singularities, applications, vol. 93 of Operator Theory: Advances and Applications, Birkhäuser Verlag, Basel, 1997.
  * [14] H.-J. Flad, R. Schneider, and B.-W. Schulze, Asymptotic regularity of solutions to Hartree-Fock equations with Coulomb potential, Math. Methods Appl. Sci., 31 (2008), pp. 2172–2201.
  * [15] D. Gontier, Mathematical contributions to the calculations of electronic structures, PhD thesis, Université Paris-Est, 2015.
  * [16] E. Hunsicker, V. Nistor, and J. O. Sofo, Analysis of periodic Schrödinger operators: regularity and approximation of eigenfunctions, J. Math. Phys., 49 (2008), pp. 083501, 21.
  * [17] F. Jollet, M. Torrent, and N. Holzwarth, Generation of Projector Augmented-Wave atomic data: A 71 element validated table in the XML format, Computer Physics Communications, 185 (2014), pp. 1246–1254.
  * [18] T. Kato, On the eigenfunctions of many-particle systems in quantum mechanics, Comm. Pure Appl. Math., 10 (1957), pp. 151–177.
  * [19] D. Koelling and G. Arbman, Use of energy derivative of the radial solution in an augmented plane wave method: application to copper, Journal of Physics F: Metal Physics, 5 (1975), p. 2041.
  * [20] V. A. Kozlov, V. Mazia, and J. Rossmann, Elliptic boundary value problems in domains with point singularities, vol. 52, American Mathematical Soc., 1997.
  * [21] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B, 54 (1996), pp. 11169–11186.
  * [22] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B, 59 (1999), pp. 1758–1775.
  * [23] P. A. Kuchment, Floquet theory for partial differential equations, vol. 60, Birkhäuser, 2012.
  * [24] A. Levitt and M. Torrent, Parallel eigensolvers in plane-wave Density Functional Theory, Computer Physics Communications, 187 (2015), pp. 98–105.
  * [25] X. Li and H. Chen, A discontinuous Galerkin scheme for full-potential electronic structure calculations, J. Comput. Phys., 385 (2019), pp. 33–50.
  * [26] Y. Maday and C. Marcati, Analyticity and hp discontinuous Galerkin approximation of nonlinear Schrödinger eigenproblems, arXiv preprint arXiv:1912.07483, (2019).
  * [27] Y. Maday and C. Marcati, Regularity and $hp$ discontinuous Galerkin finite element approximation of linear elliptic eigenvalue problems with singular potentials, Math. Models Methods Appl. Sci., 29 (2019), pp. 1585–1617.
  * [28] R. B. Melrose, The Atiyah-Patodi-Singer index theorem, vol. 4 of Research Notes in Mathematics, A K Peters, Ltd., Wellesley, MA, 1993.
  * [29] M. Reed and B. Simon, Methods of modern mathematical physics. IV. Analysis of operators, Academic Press [Harcourt Brace Jovanovich, Publishers], New York-London, 1978.
  * [30] J. C. Slater, Wave functions in a periodic potential, Physical Review, 51 (1937), p. 846.
  * [31] J. P. Solovej, Proof of the ionization conjecture in a reduced Hartree-Fock model, Invent. Math., 104 (1991), pp. 291–311.
  * [32] M. Torrent, F. Jollet, F. Bottin, G. Zérah;, and X. Gonze, Implementation of the projector augmented-wave method in the ABINIT code: Application to the study of iron under pressure, Computational Materials Science, 42 (2008), pp. 337 – 351.
  * [33] H. F. Weinberger, Variational methods for eigenvalue approximation, SIAM, Philadelphia, Pa., 1974.