A numerical study of the extended Kohn-Sham ground states of atoms

TL;DR

This paper numerically investigates the extended Kohn-Sham model for atoms, computing energy levels and responses to electric fields, and compares results with theoretical assumptions and perturbation theory.

Contribution

It provides a detailed numerical analysis of the extended Kohn-Sham model for atoms, including discretization, algorithm development, and validation against theoretical predictions.

Findings

Computed energy levels for atoms in the first four rows of the periodic table.

Observed accidental degeneracies between atomic shells.

Validated density response to electric fields against perturbation theory.

Abstract

In this article, we consider the extended Kohn-Sham model for atoms subjected to cylindrically-symmetric external potentials. The variational approximation of the model and the construction of appropriate discretization spaces are detailed together with the algorithm to solve the discretized Kohn-Sham equations used in our code. Using this code, we compute the occupied and unoccupied energy levels of all the atoms of the first four rows of the periodic table for the reduced Hartree-Fock (rHF) and the extended Kohn-Sham X$\alpha$ models. These results allow us to test numerically the assumptions on the negative spectra of atomic rHF and Kohn-Sham Hamiltonians used in our previous theoretical works on density functional perturbation theory and pseudopotentials. Interestingly, we observe accidental degeneracies between s and d shells or between p and d shells at the Fermi level of some atoms. We also consider the case of an atom subjected to a uniform electric-field. For various magnitudes of the electric field, we compute the response of the density of the carbon atom confined in a large ball with Dirichlet boundary conditions, and we check that, in the limit of small electric fields, the results agree with the ones obtained with first-order density functional perturbation theory.

Equations

$$E_{z,N}(\gamma ):=\mathrm{Tr}\left(-\frac{1}{2}\Delta \gamma \right)-z{\int }_{{\mathbb{R}}^3}\frac{{\rho }_{\gamma }}{|\cdot |}+\frac{1}{2}D({\rho }_{\gamma },{\rho }_{\gamma })+E_{\mathrm{xc}}({\rho }_{\gamma })$$

$${\mathcal{K}}_N:=\left\{\gamma \in \mathcal{S}(L^2({\mathbb{R}}^3))|0\le \gamma \le 2,\mathrm{Tr}(\gamma )=N,\mathrm{Tr}(-\Delta \gamma )<\infty \right\},$$

$${\mathfrak{S}}_{1,1}:=\left\{T\in {\mathfrak{S}}_1||\nabla |T|\nabla |\in {\mathfrak{S}}_1\right\},$$

$$\Vert T{\Vert }_{{\mathfrak{S}}_{1,1}}:=\Vert T{\Vert }_{{\mathfrak{S}}_1}+\Vert |\nabla |T|\nabla |{\Vert }_{{\mathfrak{S}}_1}.$$

$${\rho }_{\gamma }\ge 0,{\int }_{{\mathbb{R}}^3}{\rho }_{\gamma }=N,{\int }_{{\mathbb{R}}^3}|\nabla \sqrt{{\rho }_{\gamma }}{|}^2\le \mathrm{Tr}(-\Delta \gamma )<\infty .$$

$$V^{\mathrm{H}}(\rho )=\rho \star |\cdot {|}^{-1}.$$

$$\mathcal{C}:=\left\{\rho \in {\mathcal{S}}^{\prime }({\mathbb{R}}^3)|\hat{\rho }\in L_{\mathrm{loc}}^1({\mathbb{R}}^3,\mathbb{C}),|\cdot {|}^{-1}\hat{\rho }\in L^2({\mathbb{R}}^3,\mathbb{C})\right\},({\rho }_1,{\rho }_2{)}_{\mathcal{C}}=4\pi {\int }_{{\mathbb{R}}^3}\frac{{\hat{\rho }}_1(\mathbf{k}{)}^{\ast }{\hat{\rho }}_2(\mathbf{k})}{|\mathbf{k}{|}^2}d\mathbf{k},$$

$${\mathcal{C}}^{\prime }:=\left\{v\in L^6({\mathbb{R}}^3)|\nabla v\in (L^2({\mathbb{R}}^3){)}^3\right\},(v_1,v_2{)}_{{\mathcal{C}}^{\prime }}=\frac{1}{4\pi }{\int }_{{\mathbb{R}}^3}\nabla v_1\nabla v_2=\frac{1}{4\pi }{\int }_{{\mathbb{R}}^3}|\mathbf{k}{|}^2{\hat{v}}_1(\mathbf{k}{)}^{\ast }{\hat{v}}_2(\mathbf{k})d\mathbf{k}.$$

$$E_{\mathrm{xc}}^{\mathrm{rHF}}=0,$$

$$E_{\mathrm{xc}}^{\mathrm{LDA}}(\rho )={\int }_{{\mathbb{R}}^3}{ϵ}_{\mathrm{xc}}(\rho (r))dr,$$

$${ϵ}_{\mathrm{xc}}(0)=0,{ϵ}_{\mathrm{xc}}^{\prime }(0)\le 0,$$

$$\exists 0<{\beta }_{-}\le {\beta }_{+}<\frac{2}{3}\text{mbox}s.t\underset{\rho \in {\mathbb{R}}_{+}}{\sup }\frac{|{ϵ}_{\mathrm{xc}}^{\prime }(\rho )|}{{\rho }^{{\beta }_{-}}+{\rho }^{{\beta }_{+}}}<\infty ,$$

$$\exists 1\le \alpha <\frac{3}{2}\text{mbox}s.t\underset{\rho \to 0_{+}}{\mathrm{limsup}}\frac{{ϵ}_{\mathrm{xc}}(\rho )}{{\rho }^{\alpha }}<0,$$

$$\exists \lambda >-1\text{mbox}s.t{ϵ}_{\mathrm{xc}}^{\prime \prime }(\rho )\underset{\rho \to 0_{+}}{\sim }c{\rho }^{\lambda }.$$

$${\gamma }_0:=\frac{1}{4}\left({\gamma }^{\uparrow \uparrow }+\overline{{\gamma }^{\uparrow \uparrow }}+{\gamma }^{\downarrow \downarrow }+\overline{{\gamma }^{\downarrow \downarrow }}\right),$$

$${\mathcal{I}}_{z,N}^{\mathrm{rHF}}:=\inf \left\{E_{z,N}^{\mathrm{rHF}}(\gamma ),\gamma \in {\mathcal{K}}_N\right\},$$

$$E_{z,N}^{\mathrm{rHF}}(\gamma ):=\mathrm{Tr}\left(-\frac{1}{2}\Delta \gamma \right)-z{\int }_{{\mathbb{R}}^3}\frac{{\rho }_{\gamma }}{|\cdot |}+\frac{1}{2}D({\rho }_{\gamma },{\rho }_{\gamma }),$$

$${\mathcal{I}}_{z,N}^{\mathrm{LDA}}:=\inf \left\{E_{z,N}^{\mathrm{LDA}}(\gamma ),\gamma \in {\mathcal{K}}_N\right\},$$

$$E_{z,N}^{\mathrm{LDA}}(\gamma ):=\mathrm{Tr}\left(-\frac{1}{2}\Delta \gamma \right)-z{\int }_{{\mathbb{R}}^3}\frac{{\rho }_{\gamma }}{|\cdot |}+\frac{1}{2}D({\rho }_{\gamma },{\rho }_{\gamma })+E_{\mathrm{xc}}^{\mathrm{LDA}}({\rho }_{\gamma }).$$

$$H_{z,N}^{0,\mathrm{rHF}}:=-\frac{1}{2}\Delta -\frac{z}{|\cdot |}+V^{\mathrm{H}}({\rho }_{z,N}^{0,\mathrm{rHF}}),$$

$${\gamma }_{z,N}^{0,\mathrm{rHF}}=2{\text{mathds}1}_{(-\infty ,{ϵ}_{z,N,\mathrm{F}}^{0,\mathrm{rHF}})}(H_{z,N}^{0,\mathrm{rHF}})+{\delta }_{z,N}^{0,\mathrm{rHF}},$$

$${\gamma }_{z,z}^{0,\mathrm{rHF}}=2{\text{mathds}1}_{(-\infty ,{ϵ}_{z,z,\mathrm{F}}^{0,\mathrm{rHF}})}(H_{z,z}^{0,\mathrm{rHF}})$$

$$+\sum _{m=-2}^2{\gamma }_m\left(|{\varphi }_s⟩⟨{\varphi }_{d,m}|+|{\varphi }_{d,m}⟩⟨{\varphi }_s|\right),$$

$${\varphi }_s(\mathbf{r})=f_{ns}(r),{\varphi }_{d,m}(\mathbf{r})=r^2{f}_{n^{\prime }d}(r){\tilde{Y}}_2^m(\theta ,\phi ).$$

$${\alpha }^2{f}_{ns}(r{)}^2+\frac{\sqrt{15}}{\pi }{f}_{ns}(r)f_{n^{\prime }d}(r)\left({\gamma }_{-2}xy+{\gamma }_{-1}yz+{\gamma }_0\frac{2z^2-x^2-y^2}{\sqrt{3}}+{\gamma }_{1}xz+{\gamma }_2\frac{x^2-y^2}{2}\right)+\frac{15}{4\pi }{f}_{n^{\prime }d}(r{)}^2$$

$$\displaystyle\times\bigg{(}\beta_{-2,-2}x^{2}y^{2}+\beta_{-1,-1}y^{2}z^{2}+\beta_{0,0}\frac{(2z^{2}-x^{2}-y^{2})^{2}}{\sqrt{3}}+\beta_{1,1}x^{2}z^{2}+\beta_{2,2}\frac{(x^{2}-y^{2})^{2}}{4}+2\beta_{-2,-1}xy^{2}z$$

$$+{\beta }_{-2,0}\frac{xy(2z^2-x^2-y^2)}{\sqrt{3}}+2{\beta }_{-2,1}{x}^{2}yz+{\beta }_{-2,2}xy(x^2-y^2)+{\beta }_{-1,0}\frac{yz(2z^2-x^2-y^2)}{12}+2{\beta }_{-1,1}xy{z}^2$$

$$\displaystyle\qquad+\beta_{-1,2}yz(x^{2}-y^{2})+\beta_{0,1}\frac{xz(2z^{2}-x^{2}-y^{2})}{\sqrt{3}}+\beta_{0,2}\frac{(x^{2}-y^{2})(2z^{2}-x^{2}-y^{2})}{2\sqrt{3}}+\beta_{1,2}xz(x^{2}-y^{2})\bigg{)}$$

$${\gamma }_{z,z}^{0,\mathrm{rHF}}=2{\text{mathds}1}_{(-\infty ,{ϵ}_{z,z,\mathrm{F}}^{0,\mathrm{rHF}})}(H_{z,z}^{0,\mathrm{rHF}})$$

$$+\sum _{m=-1}^1\sum _{m^{\prime }=-2}^2{\gamma }_{m,m^{\prime }}\left(|{\varphi }_{p,m}⟩⟨{\varphi }_{d,m^{\prime }}|+|{\varphi }_{d,m^{\prime }}⟩⟨{\varphi }_{p,m}|\right)$$

$${\varphi }_{p,m}(\mathbf{r})=r{f}_{np}(r){\tilde{Y}}_1^m(\theta ,\varphi ),{\varphi }_{d,m}(\mathbf{r})=r^2{f}_{n^{\prime }d}(r){\tilde{Y}}_2^m(\theta ,\phi ).$$

$$\displaystyle\frac{3}{4\pi}f_{np}(r)^{2}\big{(}\alpha_{-1,-1}y^{2}+\alpha_{0,0}z^{2}+\alpha_{1,1}x^{2}+2\alpha_{-1,0}yz+2\alpha_{-1,1}xy+2\alpha_{0,1}xz\big{)}+\frac{3\sqrt{5}}{2\pi}f_{np}(r)f_{n^{\prime}d}(r)$$

$$\displaystyle\quad\times\bigg{(}\gamma_{-1,-2}xy^{2}+\gamma_{-1,-1}y^{2}z+\gamma_{-1,0}\frac{y(2z^{2}-x^{2}-y^{2})}{2\sqrt{3}}+\gamma_{-1,1}xyz+\gamma_{-1,2}\frac{y(x^{2}-y^{2})}{2}$$

$${\gamma }_{0,-2}xyz+{\gamma }_{0,-1}y{z}^2+{\gamma }_{0,0}\frac{z(2z^2-x^2-y^2)}{2\sqrt{3}}+{\gamma }_{0,1}x{z}^2+{\gamma }_{0,2}\frac{z(x^2-y^2)}{2}$$

$$\displaystyle\qquad\quad\gamma_{1,-2}x^{2}y+\gamma_{1,-1}xyz+\gamma_{1,0}\frac{x(2z^{2}-x^{2}-y^{2})}{2\sqrt{3}}+\gamma_{1,1}x^{2}z+\gamma_{1,2}\frac{x(x^{2}-y^{2})}{2}\bigg{)}+\frac{15}{4\pi}f_{n^{\prime}d}(r)^{2}$$

$$\displaystyle\times\bigg{(}\beta_{-2,-2}x^{2}y^{2}+\beta_{-1,-1}y^{2}z^{2}+\beta_{0,0}\frac{(2z^{2}-x^{2}-y^{2})^{2}}{12}+\beta_{1,1}x^{2}z^{2}+\beta_{2,2}\frac{(x^{2}-y^{2})^{2}}{4}+2\beta_{-2,-1}xy^{2}z$$

Full text

1 Introduction

This article is concerned with the numerical computation of the extended Kohn-Sham ground states of atoms for the reduced Hartree-Fock (rHF, also called Hartree) and LDA (local density approximation) models [11, 12]. We consider the case of an isolated atom, as well as the case of an atom subjected to cylindrically symmetric external potential. We notably have in mind Stark potentials, that are potentials of the form $W(\mathbf{r})=-{\mathcal{E}}\cdot\mathbf{r}$ generated by a uniform electric field ${\mathcal{E}}\neq 0$.

We first propose a method to accurately solve the extended Kohn-Sham problem for cylindrically symmetric systems, using spherical coordinates and a separation of variables. This approach is based on the fact that, for such systems, the Kohn-Sham Hamiltonian commutes with $L_{\mathbf{z}}$, the $\mathbf{z}$-component of the angular momentum operator, $\mathbf{z}$ denoting the symmetry axis of the system. We obtain in this way a family of 2D elliptic eigenvalue problems in the $r$ and $\theta$ variables, indexed by the eigenvalue $m\in{\mathbb{Z}}$ of $L_{\mathbf{z}}$, all these problems being coupled together through the self-consistent density. To discretize the 2D eigenvalue problems, we use harmonic polynomials in $\theta$ (or in other words, spherical harmonics $Y_{l}^{0}$, which only depend on $\theta$) to discretize along the angular variable, and high-order finite element methods to discretize along the radial variable $r\in[0,L_{e}]$. We then apply this approach to study numerically two kinds of systems.

First, we provide accurate approximations of the extended Kohn-Sham ground states of all the atoms of the first four rows of the periodic table. These results allow us to test numerically the assumptions on the negative spectra of atomic rHF and Kohn-Sham LDA Hamiltonians that we used in previous theoretical works on density functional perturbation theory [5] and norm-conserving semilocal pseudopotentials [6]. We show in particular that for most atoms of the first four rows of the periodic table, the Fermi level is negative and is not an accidentally degenerate eigenvalue of the rHF Hamiltonian. We also observe that there seems to be no unoccupied orbitals with negative energies. On the other hand, for some chemical elements, the Fermi level seems to be an accidentally degenerate eigenvalue (for example the rHF 5s and 4d states of the palladium atom seem to be degenerate). For a few of them, this accidentally degenerate eigenvalue is so close to zero that our calculations do not allow us to know whether it is slightly negative or equal to zero. For instance, our simulations seem to show that the 5s and 3d states of the iron atom seem to be degenerate at the rHF level of theory, and the numerical value of their energy we obtain with our code is about $-10^{-5}$ Ha.

Second, we study an atom subjected to a uniform electric field (Stark effect). In this case, the system has no ground state (the Kohn-Sham energy functional is not bounded below), but density functional perturbation theory (see [5, 6] for a mathematical analysis) can be used to compute the polarization of the electronic cloud caused by the external electric field. The polarized electronic state is not a steady state, but a resonant state, and the smaller the electric field, the longer its life time. Another way to compute the polarization of the electronic cloud is to compute the ground state for a small enough electric field in a basis set consisting of functions decaying fast enough at infinity for the electrons to stay close to the nuclei. The Gaussian basis functions commonly used in quantum chemistry satisfy this decay property. However, it is not easy to obtain very accurate results with Gaussian basis sets, since they are not systematically improvable (over-completeness issues). Here we consider instead basis functions supported in a ball $B_{L_{e}}$, where $L_{e}$ is a numerical parameter chosen large enough to obtain accurate results and small enough to prevent electrons from escaping to infinity (for a given, small, value of the external electric field ${\mathcal{E}}$). We study the ground state energy and density as functions of the cut-off radius $L_{e}$, and observe that for a given, small enough, uniform electric field, there is a plateau $[L_{e,\rm min},L_{e,\rm max}]$ on which these quantities hardly vary. For $L_{e}<L_{e,\rm min}$, the simulated system is too much confined, which artificially increases its energy, while for $L_{e}>L_{e,\rm max}$, a noticeable amount of charge accumulates at the boundary of the simulation domain, in the direction of ${\mathcal{E}}$ (where the potential energy is very negative). On the other hand, for $L_{e,\rm min}\leq L_{e}\leq L_{e,\rm max}$, the simulation provides a fairly accurate approximation of the polarization energy and the polarized density.

The article is organized as follows. In Section 2, we recall the mathematical formulation of the extended Kohn-Sham model, and some theoretical results about the rHF and LDA ground states of isolated atoms and of atoms subjected to an external cylindrically symmetric potential. In Section 3, we describe the discretization method and the algorithms used in this work to compute the extended Kohn-Sham ground states of atoms subjected to cylindrically symmetric external potentials. Some numerical results are presented in Section 4.

2 Modeling

In this article, we consider a molecular system consisting of a single nucleus of atomic charge $z\in{\mathbb{N}}^{\ast}$ and of $N$ electrons. For $N=z$, this system is the neutral atom with nuclear charge $z$, which we call atom $z$ for convenience.

2.1 Kohn-Sham models for atoms

In the framework of the (extended) Kohn-Sham model [7], the ground state energy of a system with one nucleus with charge $z$ and $N$ electrons is obtained by minimizing an energy functional of the form

[TABLE]

over the set

[TABLE]

where ${\mathcal{S}}(L^{2}({\mathbb{R}}^{3}))$ is the space of the self-adjoint operators on $L^{2}({\mathbb{R}}^{3}):=L^{2}({\mathbb{R}}^{3},{\mathbb{R}})$ and ${\rm Tr\,}(-\Delta\gamma):={\rm Tr\,}(|\nabla|\gamma|\nabla|)$. Note that, ${\mathcal{K}}_{N}$ is a closed convex subset of the space ${\mathfrak{S}}_{1,1}$ defined by

[TABLE]

endowed with norm

[TABLE]

The function $-\frac{z}{|.|}$ is the attraction potential induced on the electrons by the nucleus, and $\rho_{\gamma}$ is the density associated with the one-body density matrix $\gamma$. For $\gamma\in{\mathcal{K}}_{N}$, we have

[TABLE]

The last result is the Hoffmann-Ostenhof inequality [10]. Therefore $\sqrt{\rho_{\gamma}}\in H^{1}({\mathbb{R}}^{3})$, and in particular, $\rho_{\gamma}\in L^{1}({\mathbb{R}}^{3})\cap L^{3}({\mathbb{R}}^{3})$. For $\rho\in L^{\frac{6}{5}}({\mathbb{R}}^{3})$, $D(\rho,\rho)$ is equal to $\int_{{\mathbb{R}}^{3}}V^{\rm H}(\rho)\rho$, where $V^{\rm H}$ is the Coulomb, also called Hartree, potential generated by $\rho$:

[TABLE]

Recall that $V^{\rm H}$ can be seen as a unitary operator from the Coulomb space ${\cal C}$ to its dual ${\cal C}^{\prime}$, where

[TABLE]

and

[TABLE]

The term $E_{\rm xc}$ is the exchange-correlation energy. We will restrict ourselves to two kinds of Kohn-Sham models: the rHF model, for which the exchange-correlation energy is taken equal to zero

[TABLE]

and the Kohn-Sham LDA (local density approximation) model, for which the exchange-correlation energy has the form

[TABLE]

where $\epsilon_{\rm xc}$ is the sum of the exchange and correlation energy densities of the homogeneous electron gas. As the function $\epsilon_{\rm xc}:{\mathbb{R}}_{+}\rightarrow{\mathbb{R}}$ is not explicitly known, it is approximated in practice by an explicit function, still denoted by $\epsilon_{\rm xc}$ for simplicity. We assume here that the approximate function $\epsilon_{\rm xc}$ is a $C^{1}$ function from ${\mathbb{R}}_{+}$ into ${\mathbb{R}}_{-}$, twice differentiable on ${\mathbb{R}}^{\ast}_{+}$ and obeying the following conditions

[TABLE]

Note that these properties are satisfied by the exact function $\epsilon_{\rm xc}$. They are also satisfied by Slater’s X$\alpha$ model for which $\epsilon_{\rm xc}(\rho)=-C_{\rm D}\rho^{\frac{1}{3}}$, where $C_{\rm D}=\frac{3}{4}\left(\frac{3}{\pi}\right)^{\frac{1}{3}}$ is the Dirac constant. This model is used in the simulations reported in Section 4.

Remark 1**.**

The minimization set ${\mathcal{K}}_{N}$ defined in (2) is the set of real spin-unpolarized first-order reduced density matrices. We will call its elements non-magnetic states. The general (complex non-collinear spin-polarized, see e.g. [9]) rHF model being convex in the density matrix, and strictly convex in the density, the general rHF ground state density of a given molecular system in the absence of magnetic field, if it exists, is unique, and one of the minimizers is a non-magnetic state. Indeed, using the notation of [9], if $\gamma$ is a complex non-collinear spin-polarized ground state, the non-magnetic state

[TABLE]

where $\overline{\gamma^{\sigma,\sigma}}$ is the complex conjugate (not the adjoint) of the operator $\gamma^{\sigma,\sigma}$, is a non-magnetic ground state. The general rHF ground state energy and density of a molecular system in the absence of magnetic field can therefore be determined by minimizing the rHF energy functional over the set ${\mathcal{K}}_{N}$. The LDA model is not a priori strictly convex in the density, but it is convex over the set of complex non-collinear spin-polarized density matrices having a given density $\rho$. Therefore, the general LDA ground state energy and densities can be obtained by minimizing the LDA energy functional over the set ${\mathcal{K}}_{N}$. In contrast, this argument does not apply to the local spin density approximation (LSDA) model, whose ground states are, in general, spin-polarized.

To avoid ambiguity, for any $z$ and $N$ in ${\mathbb{R}}_{+}^{\ast}$, we denote by

[TABLE]

where

[TABLE]

and

[TABLE]

where

[TABLE]

We recall the following two theorems which ensure the existence of ground states for neutral atoms and positive ions.

Theorem 1** (ground state for the rHF model [5, 17]).**

Let $z\in{\mathbb{R}}_{+}^{\ast}$ and $N\leq z$. Then the minimization problem (9) has a ground state $\gamma^{0,\rm rHF}_{z,N}$, and all the ground states share the same density $\rho^{0,\rm rHF}_{z,N}$. The mean-field Hamiltonian

[TABLE]

is a bounded below self-adjoint operator on $L^{2}({\mathbb{R}}^{3})$, $\sigma_{\rm ess}(H^{0,\rm rHF}_{z,N})={\mathbb{R}}_{+}$, and the ground state $\gamma^{0,\rm rHF}_{z,N}$ is of the form

[TABLE]

where $\epsilon_{z,N,\rm F}^{0,\rm rHF}\leq 0$ is the Fermi level, $0\leq\delta^{0,\rm rHF}_{z,N}\leq 2$ and $\mbox{\rm Ran}(\delta^{0,\rm rHF}_{z,N})\subset\mbox{\rm Ker}(H^{0,\rm rHF}_{z,N}-\epsilon_{z,N,\rm F}^{0,\rm rHF})$. If $\epsilon_{z,N,\rm F}^{0,\rm rHF}$ is negative and is not an accidentally degenerate eigenvalue of $H^{0,\rm rHF}_{z,N}$, then the non-magnetic ground state $\gamma^{0,\rm rHF}_{z,N}$ is unique.

The numerical results presented in Section 4.1.1 indicate that, for neutral atoms, the assumption

$\epsilon_{z,z,\rm F}^{0,\rm rHF}$ is negative and is not an accidentally degenerate eigenvalue of $H^{0,\rm rHF}_{z,z}$,

which guarantees the uniqueness of the non-magnetic rHF ground state density matrix, is satisfied for all the chemical elements of the first two rows of the periodic table, and for most of the elements of the third and four rows. Surprisingly, we observe accidental degeneracies at the Fermi level for Sc and Ti (4p and 3d shells), for V, Cr, Mn and Fe (5s and 3d shells), for Zr (5p and 4d shells), Nb and Mo (6s and 4d shells), and for Pd and Ag (5s and 4d shells). For some of these elements, the Fermi level is clearly negative, and we can conclude the following:

• •

if the Fermi level contains an s and a d shell, then the non-magnetic rHF ground state is unique;

• •

if the Fermi level contains a p and a d shell, and if both shells are partially occupied (which is suggested by our numerical simulations), then the non-magnetic rHF ground state is not unique.

Indeed, in the former case, any rHF ground state is of the form

[TABLE]

where $\alpha\in{\mathbb{R}}$, $\beta\in{\mathbb{R}}^{5\times 5}_{\rm sym}$ and $\gamma\in{\mathbb{R}}^{5}$ are matrices such that 0\leq\left(\begin{array}[]{cc}\alpha&\gamma^{T}\\ \gamma&\beta\end{array}\right)\leq 2 and where

[TABLE]

Here, the $\widetilde{Y}_{l}^{m}$’s are the real spherical harmonics, and $f_{ns}$ and $f_{n^{\prime}d}$ are radial functions with respectively $(n-1)$ and $(n^{\prime}-3)$ nodes in $(0,+\infty)$. Since all the ground state density matrices share the same density, the function

[TABLE]

where $r=(x^{2}+y^{2}+z^{2})^{1/2}$, must be a function of $r$, independent of the chosen ground state density matrix. Since $f_{ns}$ has more nodes than $f_{n^{\prime}d}$ (we have seen above that $n=5$ or $6$ and $n^{\prime}=3$ or $4$), this implies that $\beta$ is a scalar matrix, that $\gamma=0$, and that only one value for the pair $(\alpha,\beta)$ is possible. This demonstrates the uniqueness of the non-magnetic ground state when the Fermi level is negative and contains a pair of accidentally degenerate s and d shells.

In the case when the Fermi level is negative and contains a pair of accidentally degenerate p and d shells, any non-magnetic ground state density matrix is of the form

[TABLE]

where $\alpha\in{\mathbb{R}}^{3\times 3}_{\rm sym}$, $\beta\in{\mathbb{R}}^{5\times 5}_{\rm sym}$ and $\gamma\in{\mathbb{R}}^{3\times 5}$ are matrices such that 0\leq\left(\begin{array}[]{cc}\alpha&\gamma\\ \gamma^{T}&\beta\end{array}\right)\leq 2 and where

[TABLE]

Here, $f_{np}$ and $f_{n^{\prime}d}$ are radial functions with respectively $(n-2)$ and $(n^{\prime}-3)$ nodes in $(0,+\infty)$. Since all the ground state density matrices share the same density, the function

[TABLE]

where $r=(x^{2}+y^{2}+z^{2})^{1/2}$, must be a function of $r$, independent of the chosen ground state density matrix. Since $f_{np}$ has more nodes than $f_{n^{\prime}d}$, this implies that $\alpha$ and $\beta$ are scalar matrices and that, for $\alpha$ and $\beta$ given, the function

[TABLE]

is a given function of $r$. The vector spaces of homogeneous polynomials in $x,y,z$ of total degree equal to $3$ is of dimension 10, and the matrix $\gamma$ has 15 independent entries. Provided $\alpha$ and $\beta$ are not equal to [math] (which is suggested by our numerical simulations), an infinity of density matrices of the form (11) satisfy the rHF equations, and are therefore admissible non-magnetic ground states.

For other chemical elements, such as iron ($z=26)$, the Fermi level is so close to zero that the numerical accuracy of our numerical method does not allow us to know whether it is slightly negative or equal to zero.

Theorem 2** (ground state for the LDA model [1]).**

Let $z\in{\mathbb{R}}_{+}^{\ast}$ and $N\leq z$. Suppose that (5)-(7) hold. Then the minimization problem (10) has a ground state $\gamma^{0,\rm LDA}_{z,N}$. In addition, $\gamma^{0,\rm LDA}_{z,N}$ satisfies the self-consistent field equation

[TABLE]

where $\epsilon_{z,N,\rm F}^{0,\rm LDA}\leq 0$ is the Fermi level, $0\leq\delta^{0,\rm LDA}_{z,N}\leq 2$, $\mbox{\rm Ran}(\delta^{0,\rm LDA}_{z,N})\subset\mbox{\rm Ker}(H^{0,\rm LDA}_{z,N}-\epsilon_{z,N,\rm F}^{0,\rm LDA})$ and the mean-field Hamiltonian

[TABLE]

where $\rho^{0,\rm LDA}_{z,N}=\rho_{\gamma^{0,\rm LDA}_{z,N}}$ and $v_{\rm xc}(\rho)=\frac{d\,\epsilon_{\rm xc}}{d\,\rho}(\rho)$, is a bounded below self-adjoint operator on $L^{2}({\mathbb{R}}^{3})$ and $\sigma_{\rm ess}(H^{0,\rm LDA}_{z,N})={\mathbb{R}}_{+}$.

2.2 Density functional perturbation theory

We now examine the response of the ground state density matrix when an additional external potential $\beta W$ is turned on. The energy functional to be minimized over ${\mathcal{K}}_{N}$ now reads

[TABLE]

and is well-defined for any $\gamma\in{\mathcal{K}}_{N}$, $W\in{\mathcal{C}}^{\prime}$ and $\beta\in{\mathbb{R}}$. The parameter $\beta$ is called the coupling constant in quantum mechanics. Denote by

[TABLE]

The following theorem insures the existence of a perturbed ground state density matrix for perturbation potentials in ${\mathcal{C}}^{\prime}$.

Theorem 3** (existence of a perturbed minimizer [5]).**

Let $z\in{\mathbb{R}}_{+}^{\ast}$, $N\leq z$ and $W\in\cal C^{\prime}$. Assume that the Fermi level $\epsilon_{z,N,\rm F}^{0,\rm rHF}$ is negative and is not an accidentally degenerate eigenvalue of $H^{0,\rm rHF}_{z,N}$. Then the non-magnetic unperturbed rHF ground state, that is the minimizer of (9), is unique, and the perturbed problem (14) has a unique non-magnetic ground state $\gamma_{z,N,\beta W}^{\rm rHF}$, for $\beta\in{\mathbb{R}}$ small enough. The Hamiltonian

[TABLE]

where $\rho_{z,N,\beta W}^{\rm rHF}=\rho_{\gamma_{z,N,\beta W}^{\rm rHF}}$, is a bounded below self-adjoint operator on $L^{2}({\mathbb{R}}^{3})$ with form domain $H^{1}({\mathbb{R}}^{3})$ and $\sigma_{\rm ess}(H^{0,\rm rHF}_{z,N,\beta W})={\mathbb{R}}_{+}$. Moreover, $\gamma_{z,N,\beta W}^{\rm rHF}$ and $\rho_{z,N,\beta W}^{\rm rHF}$ are analytic in $\beta$, that is

[TABLE]

the above series being normally convergent in ${\mathfrak{S}}_{1,1}$ and ${\mathcal{C}}$ respectively.

In the sequel, we will refer to $\gamma_{z,N,W}^{(k)}$ as the $k$-th order perturbation of the density matrix.

Although we focus here on non-magnetic states, it is convenient to consider $H_{z,N}^{0,\rm rHF}$ as an operator on $L^{2}({\mathbb{R}}^{3},{\mathbb{C}})$ in order to expand the angular part of the atomic orbitals on the usual complex spherical harmonics. It would of course have been possible to avoid considering complex wave functions by expanding on real spherical harmonics. However, we have chosen to work with complex wave function to prepare the ground for future works on magnetic systems.

The unperturbed Hamiltonian $H_{z,N}^{0,\rm rHF}$ is a self-adjoint operator on $L^{2}({\mathbb{R}}^{3},{\mathbb{C}})$ invariant with respect to rotations around the nucleus (assumed located at the origin). This operator is therefore block-diagonal in the decomposition of $L^{2}({\mathbb{R}}^{3},{\mathbb{C}})$ as the direct sum of the pairwise orthogonal subspaces ${\mathcal{H}}_{l}:=\mbox{Ker}({\mathbf{L}}^{2}-l(l+1))$:

[TABLE]

where ${\mathbf{L}}=\mathbf{r}\times(-i\nabla)$ is the angular momentum operator. Since we are going to consider perturbation potentials which are not spherically symmetric, but only cylindrically symmetric, or in other words independent of the azimuthal angle $\varphi$ in spherical coordinates, the ${\mathcal{H}}_{l}$’s are no longer invariant subspaces of the perturbed Hamiltonians. The appropriate decomposition of $L^{2}({\mathbb{R}}^{3},{\mathbb{C}})$ into invariant subspaces for Hamiltonians $H_{z,N,\beta W}^{\rm rHF}$ with $W$ cylindrically symmetric, is the following: for $m\in{\mathbb{Z}}$, we set

[TABLE]

where $L_{\mathbf{z}}$ is the $\mathbf{z}$-component of the angular momentum operator ${\mathbf{L}}$ ($L_{\mathbf{z}}={\mathbf{L}}.\mathbf{e}_{\mathbf{z}}$).

Note that

[TABLE]

and

[TABLE]

where $Y_{l}^{m}$ are the spherical harmonics, i.e. the joint eigenfunctions of $\Delta_{S}$, the Laplace-Beltrami operator on the unit sphere ${\mathbb{S}}^{2}$ of ${\mathbb{R}}^{3}$, and ${\mathcal{L}}_{\mathbf{z}}$, the generator of rotations about the azimuthal axis of ${\mathbb{S}}^{2}$. More precisely, we have

[TABLE]

where, in spherical coordinates,

[TABLE]

These functions are orthonormal, in the following sense:

[TABLE]

where $\delta_{ij}$ is the Kronecker symbol and $(Y_{l}^{m})^{*}=(-1)^{m}Y_{l}^{-m}$ is the complex conjugate of $Y_{l}^{m}$.

We also define

[TABLE]

so that $L^{2}({\mathbb{R}}^{3},{\mathbb{C}})$ and $H^{1}({\mathbb{R}}^{3},{\mathbb{C}})$ are decomposed as the following direct sums:

[TABLE]

each ${\mathcal{H}}^{m}$ being $H_{z,N,\beta W}^{\rm rHF}$-stable (in the sense of unbounded operators) for $W$ cylindrically symmetric. This is due to the fact that, for $W$ cylindrically symmetric, the operator $H_{z,N,\beta W}^{\rm rHF}$ commutes with $L_{\mathbf{z}}$. Note that $\sigma(H_{z,N,\beta W}^{\rm rHF})=\overline{\underset{m\in{\mathbb{Z}}}{\cup}\sigma\left(H_{z,N,\beta W}^{\rm rHF}|_{{\cal H}^{m}}\right)}$. Same arguments hold true for $H_{z,N,\beta W}^{\rm LDA}$ under the assumption that the ground state density $\rho_{z,N,\beta W}^{0,\rm LDA}$ is cylindrically symmetric (which is the case whenever it is unique).

We are interested in the Stark potential

[TABLE]

which does not belong to ${\mathcal{C}}^{\prime}$, and thus does not fall into the scope of Theorem 3. We therefore introduce the classes of perturbation potentials

[TABLE]

where ${\mathcal{H}}^{0}_{\rm loc}:={\mathcal{H}}^{0}\cap L^{2}_{\rm loc}({\mathbb{R}}^{3})$, which contain the Stark potential $W_{\rm Stark}$ whenever $s>5/2$. For $W\in{\mathcal{W}}_{s}\setminus{\mathcal{C}}^{\prime}$, the energy functional (13) is not necessarily bounded below on ${\mathcal{K}}_{N}$ for $\beta\neq 0$. Thus the solution of (14) may not exist. This is the case for the Stark potential $W_{\rm Stark}$. However, the $k$-th order perturbation of the ground state may exist, as this is the case when the linear Schrödinger operator of the hydrogen atom is perturbed by the Stark potential $W_{\rm Stark}$ (see e.g [14]). The following theorem ensures the existence of the first order perturbation of the density matrix.

Theorem 4** (first order density functional perturbation theory [6]).**

Let $z\in{\mathbb{R}}_{+}^{\ast}$, $0<N\leq z$, such that $\epsilon_{z,N,\rm F}^{0,\rm rHF}$ is negative111Note that, $\epsilon_{z,N,\rm F}^{0,\rm rHF}<0$ whenever $0<N<z$ (see e.g. [17]). and is not an accidentally degenerate eigenvalue of $H_{z,N}^{0,\rm rHF}$, $s\in{\mathbb{R}}$ and $W\in{\mathcal{W}}_{s}$. In the rHF framework, the first order perturbation of the density matrix $\gamma_{z,N,W}^{(1),\rm rHF}$ is well defined in ${\mathfrak{S}}_{1,1}$.

Note that assumption (8) is used to establish the existence and uniqueness of the first order perturbation of the density matrix $\gamma_{z,N,W}^{(1),\rm LDA}$ in ${\mathfrak{S}}_{1,1}$.

3 Numerical method

In this section, we present the discretization method and the algorithms we used to calculate numerically the ground state density matrices for (9), (10) and (14) for cylindrically symmetric perturbation potentials $W$, together with the ground state energy and the lowest eigenvalues of the associated Kohn-Sham operator. From now on, we make the assumption that the ground state density of (14), if it exists, is cylindrically symmetric which is always the case for the rHF model. Using spherical coordinates, we can write

[TABLE]

(since $Y_{l}^{0}$ is independent of $\varphi$, we use the notation $Y_{l}^{0}(\theta)$ instead of $Y_{l}^{0}(\theta,\varphi)$). As the ground state density $\rho_{z,N,\beta W}$ is assumed to be cylindrically symmetric as well, one has

[TABLE]

The Hartree and the exchange-correlation potentials also have the same symmetry. For $\rho\in L^{1}({\mathbb{R}}^{3})\cap L^{3}({\mathbb{R}}^{3})\cap{\mathcal{H}}^{0}$, we have

[TABLE]

where, for each $l\geq 0$, $V^{\rm H}_{\rho_{l}}(r)$ solves the following differential equation

[TABLE]

with boundary conditions

[TABLE]

while $(v_{\rho}^{\rm xc})_{l}$ can be computed by projection on the spherical harmonics $Y_{l}^{0}$:

[TABLE]

3.1 Discretisation of the Kohn-Sham model

Recall that for $W\in{\mathcal{W}}_{s}$ and $\beta\neq 0$, the energy functional defined by (13) is not necessarily bounded below on ${\mathcal{K}}_{N}$, which implies in particular that (14) may have no ground state. Nevertheless, one can compute approximations of (14) in finite-dimensional spaces, provided that the basis functions decay fast enough at infinity. Let $N_{h}\in{\mathbb{N}}^{\ast}$ and $m_{h}\geq m_{z}^{*}:=\max\{m|\,\exists k>0;\,\epsilon_{m,k}^{0}\leq\epsilon_{z,N,\rm F}^{0}\}$, and let $\{{\mathcal{X}}_{i}\}_{1\leq i\leq N_{h}}\in\left(H^{1}_{0}(0,+\infty)\right)^{N_{h}}$ be a free family of real-valued basis functions. We then introduce the finite-dimensional spaces

[TABLE]

and

[TABLE]

and the set

[TABLE]

Note that since our goal is to compute non-magnetic ground states, we are allowed to limit ourselves to real linear combinations in (19) and (20).

3.1.1 Variational approximation

A variational approximation of (14) is obtained by minimizing the energy functional (13) over the approximation set ${\mathcal{K}}_{N,h}$:

[TABLE]

Any $\gamma_{h}\in{\mathcal{K}}_{N,h}$ can be written as

[TABLE]

with

[TABLE]

[TABLE]

The functions $\Phi_{m,k,h}$ being in ${\mathcal{V}}^{m,h}$, they are of the form

[TABLE]

where for each $-m_{h}\leq m\leq m_{h}$, $1\leq k\leq(m_{h}-|m|+1)N_{h}$ and $|m|\leq l\leq m_{h}$, $u_{l}^{m,k,h}\in{\mathcal{X}}^{h}$. Note that $u_{l}^{-m,k,h}=u_{l}^{m,k,h}$. Expanding the functions $u_{l}^{m,k,h}$ in the basis $({\mathcal{X}}_{i})_{1\leq i\leq N_{h}}$ as

[TABLE]

and gathering the coefficients $U_{i,l}^{m,k}$ for fixed $m$ and $k$ in a rectangular matrix $U^{m,k}\in{\mathbb{R}}^{N_{h}\times(m_{h}-|m|+1)}$, any $\gamma_{h}\in{\mathcal{K}}_{N,h}$ can be represented via (22)-(24) by at least one element of the set

[TABLE]

where

[TABLE]

and

[TABLE]

The matrix $M_{0}$ appearing in the definition of ${\mathcal{U}}_{h}$ is the mass matrix defined by

[TABLE]

and the constraints ${\rm Tr\,}([U^{m,k}]^{T}M_{0}U^{m,k^{\prime}})=\delta_{kk^{\prime}}$ come from the fact that

[TABLE]

Remark 2**.**

An interesting observation is that, if there is no accidental degeneracy in the set of the occupied energy levels of $H_{z,N}^{0,\rm rHF/LDA}$, and if the occupied orbitals are well enough approximated in the space ${\mathcal{V}}^{m,h}$, then the approximate ground state density matrix $\gamma^{0,\rm rHF/LDA}_{z,N,h}$ has a unique representation of the form (22)-(24), up to the signs and the numbering of the functions $u_{l}^{m,k,h}$, that is up to the signs and numbering of the column vectors of the matrices $U^{m,k}$. By continuity, this uniqueness of the representation will survive if a small-enough cylindrically-symmetric perturbation is turned on. This is the reason why this representation is well-suited to our study.

Let us now express each component of the energy functional $\widetilde{E}_{z,N}^{\rm rHF,LDA}(\gamma_{h},\beta W)$ using the representation (22)-(24) of the elements of ${\mathcal{K}}_{N,h}$. For this purpose, we introduce the $N_{h}\times N_{h}$ real symmetric matrices $A$ and $M_{n}$, $n=-2,-1,0,1$ with entries

[TABLE]

The weighted mass matrices $M_{-2}$ and $M_{-1}$ are well-defined in view of the Hardy inequality

[TABLE]

We assume from now on that the basis functions ${\mathcal{X}}_{i}$ decay fast enough at infinity for the weighted mass matrix $M_{1}$ to be well-defined.

In the representation (22)-(24), the kinetic energy is equal to

[TABLE]

where $D_{m}\in{\mathbb{R}}^{(m_{h}-|m|+1)\times(m_{h}-|m|+1)}$ is the diagonal matrix defined by

[TABLE]

All the other terms in the energy functional depending on the density

[TABLE]

we first need to express this quantity as a function of the matrices $U^{m,k}$ and the occupation numbers $n_{m,k}$. As the function $\rho_{h}$ is in ${\mathcal{H}}^{0}$, we have

[TABLE]

Inserting (23) in (29), we get

[TABLE]

We recall the following equality [15]

[TABLE]

with

[TABLE]

where $\displaystyle\left(\begin{matrix}l_{1}&l_{2}&l_{3}\\ m_{1}&m_{2}&m_{3}\end{matrix}\right)$ denote the Wigner 3j-symbols. Inserting the expansion (24) in (31) and using (32) and the fact that

[TABLE]

we obtain

[TABLE]

from which we conclude that

[TABLE]

For $0\leq l\leq 2m_{h}$, we introduce the matrix $R_{l}\in{\mathbb{R}}^{N_{h}\times N_{h}}$ defined by

[TABLE]

where $C^{l,m}\in{\mathbb{R}}^{(m_{h}-|m|+1)\times(m_{h}-|m|+1)}$ is the symmetric matrix222The symmetry of the matrix $C^{lm}$ comes from the following symmetry properties of the 3j-symbols:

$\left(\begin{matrix}l_{1}&l_{2}&l_{3}\\ m_{1}&m_{2}&m_{3}\end{matrix}\right)=(-1)^{l_{1}+l_{2}+l_{3}}\left(\begin{matrix}l_{2}&l_{1}&l_{3}\\ m_{2}&m_{1}&m_{3}\end{matrix}\right)=(-1)^{l_{1}+l_{2}+l_{3}}\left(\begin{matrix}l_{2}&l_{1}&l_{3}\\ -m_{2}&-m_{1}&-m_{3}\end{matrix}\right).$

defined by

[TABLE]

so that

[TABLE]

Note that $C^{0,m}$ is the identity matrix, so that

[TABLE]

and

[TABLE]

and that $C^{1,m}$ is a symmetric tridiagonal matrix whose diagonal elements all are equal to zero.

The Coulomb attraction energy between the nucleus and the electrons then is equal to

[TABLE]

where we have used the orthonormality condition (16) and the fact that $Y_{0}^{0}=\frac{1}{\sqrt{4\pi}}$.

Likewise, since $Y_{1}^{0}(\theta)=\sqrt{\frac{3}{4\pi}}\cos(\theta)$, the Stark potential (18) can be written in spherical coordinates as

[TABLE]

and the potential energy due to the external electric field is then equal to

[TABLE]

Let $\mu$ be a radial, continuous function from ${\mathbb{R}}^{3}$ to ${\mathbb{R}}$ vanishing at infinity and such that $\int_{{\mathbb{R}}^{3}}\mu=1$. The Coulomb interaction energy can be rewritten as follows:

[TABLE]

The reason why we introduce the charge distribution $\mu$ is to make neutral the charge distributions $\rho_{h}-\left(\int_{{\mathbb{R}}^{3}}\rho_{h}\right)\mu$ in the first term of the right-hand side of (36), in such a way that the physical solution $Q_{0,R_{0}}$ to the equation (39) below for $l=0$ is in $H^{1}_{0}(0,+\infty)$.

Introducing the real symmetric matrix $V_{\mu}\in{\mathbb{R}}^{N_{h}\times N_{h}}$ with entries

[TABLE]

where $\mathbf{e}$ is any unit vector of ${\mathbb{R}}^{3}$ (the value of $V^{\rm H}(\mu)(r\mathbf{e})$ is independent of $\mathbf{e}$ since $V^{\rm H}(\mu)$ is radial) the sum of the last two terms of the right-hand side of (36) can be rewritten as

[TABLE]

Denoting by

[TABLE]

we have by symmetry $\widetilde{V}^{\rm H}(\rho_{h})\in{\mathcal{H}}^{0}$ and

[TABLE]

where $Q_{l,R_{l}}$ is the unique solution in $H^{1}_{0}(0,+\infty)$ to the differential equation

[TABLE]

Note that the mappings $R_{l}\mapsto Q_{l,R_{l}}$ are linear. We therefore obtain

[TABLE]

Finally, the exchange-correlation energy is

[TABLE]

3.1.2 Approximation of the Hartree term

Except for very specific basis functions (such as Gaussian atomic orbitals), it is not possible to evaluate exactly the first contribution to the Coulomb energy (39). It is therefore necessary to approximate it. For this purpose, we use a variational approximation of (38)-(39) in an auxiliary basis set $\left\{\zeta_{p}\right\}_{1\leq p\leq N_{h,\rm a}}\in(H^{1}_{0}(0,+\infty))^{N_{h,\rm a}}$, which amounts to replacing $\frac{1}{2}D(\rho_{h},\rho_{h})$ by its lower bound

[TABLE]

where $Q_{l,R_{l}}^{h}$ is the unique solution in $\zeta^{h}=\mbox{span}(\zeta_{1},\cdots,\zeta_{N_{h,\rm a}})$ to the problem

[TABLE]

which is nothing but the variational approximation of (38) in the finite dimensional space $\zeta^{h}$. Expanding the functions $Q_{l,R_{l}}^{h}$ in the basis set $\left\{\zeta_{k}\right\}_{1\leq k\leq N_{h,\rm a}}$ as

[TABLE]

and collecting the coefficients $Q_{p,l}$, $1\leq k\leq N_{h,\rm a}$ in a vector $Q_{l}\in{\mathbb{R}}^{N_{h,\rm a}}$, we obtain that the vector $Q_{l}$ is solution to the linear system

[TABLE]

where the $N_{h,\rm a}\times N_{h,\rm a}$ real symmetric matrices $A^{\rm a}$ and $M_{-2}^{\rm a}$ are defined by

[TABLE]

where $F\in{\mathbb{R}}^{N_{h,\rm a}\times N_{h}\times N_{h}}$ is the three-index tensor with entries

[TABLE]

and where $G\in{\mathbb{R}}^{N_{h,\rm a}}$ is the vector with entries

[TABLE]

Note that since $N={\rm Tr\,}(M_{0}R_{0})$, the mappings $R_{l}\mapsto Q_{l}$ are in fact linear. We finally get

[TABLE]

where $Q_{l}$ is the solution to (42).

3.1.3 Final form of the discretized problem and Euler-Lagrange equations

We therefore end up with the following approximation of problem (14):

[TABLE]

where

[TABLE]

where for each $l$, the matrix $R_{l}$ and the vector $Q_{l}$ are respectively defined by (33) and (42), and where the last but one term in the right-hand side is given by (40).

The gradient of ${\mathcal{E}}_{z,N,\beta}^{\rm rHF/LDA}$ with respect to $U^{m,k}$ is

[TABLE]

where for each $0\leq l\leq 2m_{h}$, the $N_{h}\times N_{h}$ real matrix $V_{\rm xc}^{l}$ is defined by

[TABLE]

where $v_{\rm xc}(\rho):=\frac{d\epsilon_{\rm xc}}{d\rho}(\rho)$ is the exchange-correlation potential.

Diagonalizing simultaneously the Kohn-Sham Hamiltonian and the ground state density matrix in an orthonormal basis, we obtain that the ground state can be obtained by solving the following system of first-order optimality conditions, which is nothing but a reformulation of the discretized extended Kohn-Sham equations exploiting the cylindrical symmetry of the problem:

[TABLE]

where the matrices $A$, $M_{-2}$, $M_{-1}$, $M_{0}$, $M_{1}$, $D_{m}$, $V_{\mu}$, $A^{\rm a}$, $M_{-2}^{\rm a}$, $C^{l,m}$, the 3-index tensor $F$ and the vector $G$ are defined by (26), (28), (34), (37), (43), (44), (45).

3.1.4 ${\mathbb{P}}_{4}$-finite element method

In our calculations, we use the same approximation space to discretize the radial components of the Kohn-Sham orbitals and the radial Poisson equations (38), so that, in our implementation of the method, $N_{h,\rm a}=N_{h}$ and ${\mathcal{X}}^{h}=\zeta^{h}$. We choose a cut-off radius $L_{e}>0$ large enough and discretize the interval $[0,L_{e}]$ using a non-uniform grid with $N_{I}+1$ points $0=r_{1}<r_{2}<\cdots<r_{N_{I}}<r_{N_{I}+1}=L_{e}$. The positions of the points are chosen according to the following rule:

[TABLE]

where $0<s<1$ is a scaling parameter leading to a progressive refinement of the mesh when one gets closer to the nucleus ($r=0)$. To achieve the desired accuracy, we use the ${\mathbb{P}}_{4}$-finite element method.

All the terms in the variational discretization of the energy and of the constraints can be computed exactly (up to finite arithmetics errors), except the exchange-correlation terms (40) and (48), which requires a numerical quadrature method. In our calculation, we use Gaussian quadrature formulas [18] of the form

[TABLE]

where the $0<t_{1,r}<\cdots<t_{N_{{\rm g},r},r}<1$ (resp. $-1<t_{1,\theta}<\cdots<t_{N_{{\rm g},\theta},\theta}<1$) are Gauss points for the $r$-variable (resp. for the $t_{\theta}$-variable) with associated weights $w_{1,r},\cdots,w_{N_{{\rm g},r},r}$ (resp. $w_{1,\theta},\cdots,w_{N_{{\rm g},\theta},\theta}$).

More details about the practical implementation of the method are provided in Appendix.

3.2 Description of the algorithm

In order to solve the self-consistent equations (49)-(55), we use an iterative algorithm. For clarity, we first present this algorithm within the continuous setting. Its formulation in the discretized setting considered here is detailed below. The iterations are defined as follows: an Ansatz of the ground state density $\rho^{[n]}$ being known,

construct the Kohn-Sham operator

[TABLE]

where $v_{\rm xc}=0$ for the rHF model and $v_{\rm xc}=v_{\rm xc}^{\rm LDA}$ for the Kohn-Sham LDA model; 2. 2.

for each $m\in{\mathbb{Z}}$, compute the negative eigenvalues of $H^{[n]}_{m}:=\Pi_{m}H^{[n]}\Pi_{m}$, where $\Pi_{m}$ is the orthogonal projector on the space ${\mathcal{H}}^{m}$:

[TABLE] 3. 3.

construct a new density

[TABLE]

where

[TABLE] 4. 4.

update the density:

[TABLE]

where $t_{n}\in[0,1]$ either is a fixed parameter independent of $n$ and chosen a priori, or is optimized using the Optimal Damping Algorithm (ODA), see below; 5. 5.

if some convergence criterion is satisfied, then stop; else, replace $n$ with $n+1$ and go to step 1.

In the non-degenerate case, that is when $\epsilon_{\rm F}^{[n+1]}$ is not an eigenvalue of the Hamiltonian $H^{[n]}$, the occupation numbers $n_{m,k}^{[n+1]}$ are equal to either [math] (unoccupied) or $2$ (fully occupied), while in the degenerate case the occupation numbers at the Fermi level have to be determined. We distinguish two cases: if $W=0$, or more generally if $W$ is spherically symmetric, and if $\epsilon_{\rm F}^{[n+1]}$ is not an accidentally degenerate eigenvalue of $H^{[n]}$, then the occupation numbers at the Fermi level are all equal; otherwise, the occupation numbers are not known a priori. In our approach we select the occupation numbers at the Fermi level which provide the lowest Kohn-Sham energy. When the degenerate eigenspace at the Fermi level is of dimension $3$, that is when the highest energy partially occupied orbitals are perturbations of a three-fold degenerate p-orbital, the optimal occupation numbers can be found by using the golden search or bisection method [13, Chapter 10] since, in this case, the search space can be parametrized by a single real-valued parameter (this is due to the fact that the sum of the three occupation numbers is fixed and that two of them are equal by cylindrical symmetry). In the general case, more generic optimization methods have to be resorted to.

In the discretization framework we have chosen, the algorithm can be formulated as follows.

Initialization.

Choose the numerical parameters $m_{h}$ (cut-off in the spherical harmonics expansion), $L_{e}$ (size of the simulation domain for the radial components of the Kohn-Sham orbitals and the electrostatic potential), $N_{I}$ (size of the mesh for solving the radial equations), $N_{{\rm g},r}$ (number of Gauss points for the radial quadrature formula), $N_{{\rm g},\theta}$ (number of Gauss points for the angular quadrature formula), and $\varepsilon>0$ (convergence threshold), 2. 2.

assemble the matrices $A=A^{\rm a}$, $M_{-2}=M_{-2}^{\rm a}$, $M_{-1}$, $M_{0}$, $M_{1}$, $C^{l,m}$, $V_{\mu}$ and the vector $G$. The tensor $F$ can be either computed once and for all, or the contractions $F:R_{l}^{[n]}$ can be computed on the fly, depending on the size of the discretization parameters and the computational means available; 3. 3.

choose an initial guess $(R_{l}^{[0]})_{0\leq l\leq 2m_{h}}$ for the matrices representing the discretized ground state density at iteration [math] (it is possible to take $R_{l}=0$ for all $l$ if no other better guess is known).

Iterations. The matrices $(R_{l}^{[n]})_{0\leq l\leq 2m_{h}}$ at iteration $n$ being known,

construct the building blocks of the discretized analogues of the operators $H^{[n]}_{m}$. For this purpose,

1. (a)

solve, for each $l=0,\cdots,2m_{h}$, the linear equation

[TABLE] 2. (b)

assemble, for each $l=0,\cdots,2m_{h}$, the matrix $V_{l}^{{\rm xc},[n]}$ by means of use Gauss quadrature rules

[TABLE]

where

[TABLE] 2. 2.

solve, for each ${0}\leq m\leq m_{h}$, the generalized eigenvalue problem

[TABLE] 3. 3.

build the matrices $R_{l,*}^{[n+1]}$ using the Aufbau principle and, if necessary, optimizing the occupation numbers $n_{m,k}^{[n+1]}$, by selecting the occupation numbers at the Fermi level leading to the lowest Kohn-Sham energy333In practice, this optimization problem is low-dimensional. Indeed, the degeneracy of the Fermi level is typically 3 (perturbation of p-orbitals) or 5 (perturbation of d-orbitals) for most atoms of the first four rows of the periodic table, and some of the occupation numbers are known to be equal for symmetric reasons.:

[TABLE]

where

[TABLE] 4. 4.

update the density:

[TABLE]

where $t_{n}\in[0,1]$ either is a fixed parameter independent of $n$ and chosen a priori, or is optimized using the ODA, see below; 5. 5.

if (for instance) $\max_{0\leq l\leq 2m_{h}}\|R_{l}^{[n+1]}-R_{l}^{[n]}\|\leq\epsilon$ or $|E^{[n+1]}-E^{[n]}|\leq\varepsilon$ then stop; else go to step one.

Note that the generalized eigenvalue problem (2)-(57) can be rewritten as a standard generalized eigenvalue problem of the form

[TABLE]

where the unknowns are vectors (and not matrices) by introducing the column vectors ${\mathbb{V}}_{k}\in{\mathbb{R}}^{(m_{h}+1-|m|)N_{h}}$ and the block matrices

[TABLE]

defined as

[TABLE]

and

[TABLE]

where each of the $(m_{h}-|m|+1)$ block ${\mathbb{H}}^{m}_{l,l^{\prime}}$ is of size $N_{h}\times N_{h}$ with

[TABLE]

[TABLE]

If $\beta=0$ and if the density $\rho_{h}^{[n]}$ is radial, then $R_{l}^{[n]}=0$ for all $l\in{\mathbb{N}}^{\ast}$, and the matrix ${\mathbb{H}}^{m}$ is block diagonal. The generalized eigenvalue problem (58) can then be decoupled in $(m_{h}-|m|+1)$ independent generalized eigenvalue problems of size $N_{h}$. This comes from the fact that the problem being spherically symmetric, the Kohn-Sham Hamiltonian is block diagonal in the two decompositions

[TABLE]

Let us conclude this section with some remarks on the Optimal Damping Algorithm (ODA) [3, 4], used to find an optimal step-length $t_{n}$ to mix the matrices $R_{l,*}^{[n+1]}$ and $R_{l}^{[n]}$ in Step 4 of the iterative algorithm. This step-length is obtained by minimizing on the range $t\in[0,1]$ the one-dimensional function

[TABLE]

where $\gamma^{[n]}$ is the current approximation of the ground state density matrix at iteration $n$ and

[TABLE]

with

[TABLE]

A key observation is that this optimization problem can be solved without storing density matrices, but only the two sets of matrices $R^{[n]}:=(R_{l}^{[n]})_{0\leq l\leq 2m_{h}}$ and $R_{*}^{[n+1]}:=(R_{l,*}^{[n+1]})_{0\leq l\leq 2m_{h}}$, and the scalars

[TABLE]

and

[TABLE]

Indeed, we have for all $t\in[0,1]$,

[TABLE]

where the functional ${\cal F}^{\rm rHF/LDA}$ collects all the terms of the Kohn-Sham functional depending on the density only. When $E_{\rm xc}=0$ (rHF model), the function

[TABLE]

is a convex polynomial of degree two, and its minimizer on $[0,1]$ can therefore be easily computed explicitly. In the LDA case, the minimum on $[0,1]$ of the above function of $t$ can be obtained using any line search method. We use here the golden search method. Once the minimizer $t_{n}$ is found, the quantity $E^{[n]}_{\rm kin}$ is updated using the relation

[TABLE]

4 Numerical results

As previously mentioned, we use in our code, written in Fortran 95 language, the same basis to discretize the radial components of the Kohn-Sham orbitals and of the Hartree potential, that is $({\mathcal{X}}_{i})_{1\leq i\leq N_{h}}=(\zeta_{i})_{1\leq i\leq N_{h}}$, and the ${\mathbb{P}}_{4}$ finite elements method to construct the discretization basis.

In order to test our methodology on LDA-type models, we have chosen to work with the X$\alpha$ model [16], which has a simple analytic expression:

[TABLE]

The exchange-correlation contributions must be computed by numerical quadratures. We use here the Gauss quadrature method with $N_{{\rm g},r}=15$ and $N_{{\rm g},\theta}=30$ (see Section 3.1.4).

We start this section by studying the convergence rate of the ground state energy and of the occupied energy levels of the carbon atom ($z=6$) as functions of the cut-off radius $L_{e}$ and the mesh size $N_{I}$ (see Section 3.1.4). The errors on the total energy and on the occupied energy levels for the rHF and ${\rm X}\alpha$ models are plotted in Fig. 1 (for $L_{e}=50$ and different values of $N_{I}$) and Fig. 2 (for $N_{I}=50$ and different values of $L_{e}$), the reference calculation corresponding to $L_{e}=100$ and $N_{I}=100$. We can see that the choice $L_{e}=50$ and $N_{I}=50$ provide accuracies of about $1\;\mu$Ha (recall that chemical accuracy corresponds to 1 mHa).

4.1 Electronic structures of isolated atoms

We report here calculations on all the atoms of the first four rows of the periodic table obtained with the rHF (Section 4.1.1) and X$\alpha$ (Section 4.1.2) models respectively.

4.1.1 Occupied energy levels in the rHF model

The negative eigenvalues of $H_{\rho^{0}}^{\rm rHF}$ for all $1\leq z\leq 54$ (first four rows of the periodic table) are listed in the tables below. The results for $1\leq z\leq 20$, $27\leq z\leq 39$, $43\leq z\leq 45$ and $48\leq z\leq 54$ correspond to $N_{I}$ increasing from $35$ to $75$ as $z$ increases and $L_{e}$ increasing from $30$ to $100$ as $|\epsilon_{z,z,\rm F}^{0,\rm rHF}|$ decreases, which were sufficient to obtain an accuracy of $1\;\mu$Ha. The remaining atoms are more difficult to deal with because the Fermi level seems to be an accidentally degenerate eigenvalue of $H_{\rho^{0}}^{\rm rHF}$ associated with

• •

the 4p and 3d shells for $z=21$ and $z=22$;

• •

the 5s and 3d shells for $23\leq z\leq 26$, with a Fermi level very close (or possibly equal) to zero;

• •

the 5p and 4d shells for $z=40$, with a Fermi level very close (or possibly equal) to zero;

• •

the 6s and 4d shells for $z=41$ and $z=42$, with a Fermi level very close (or possibly equal) to zero;

• •

the 5s and 4d shells for $z=46$ and $z=47$.

Since the radial component of the highest occupied orbital typically vanishes as $e^{-\sqrt{2|\epsilon_{z,z,\rm F}^{0,\rm rHF}|}r}$ if $\epsilon_{z,N,\rm F}^{0,\rm rHF}<0$ and algebraically if $\epsilon_{z,z,\rm F}^{0,\rm rHF}=0$, a very large value of $L_{e}$ is needed for the atoms for which the Fermi level is very close or possibly equal to zero. For that case, we use a non-uniform grid with $N_{I}^{\prime}=80$ and $L_{e}^{\prime}=100$ as explained in Section 3.1.4 and glue to it a uniform one with $10$ points and length $L_{e}-L_{e}^{\prime}$ varying from $70$ to $700$. Lastly, we add to the basis a function with an unbounded support, equal to $L_{e}/r$ on $[L_{e},+\infty)$ (see Appendix for details). This was sufficient to obtain an accuracy of $10\;\mu$Ha.

When the accidental degeneracy involves an $s$-shell and since the density is radial, the problem of finding the occupation numbers at the Fermi level reduces to finding a single parameter $t_{0}\in[0,1]$, which encodes the amount of electrons on the upper $s$-shell. In other words, one can write

[TABLE]

where $\rho_{\rm f}$ is the density corresponding to the fully occupied shells, and where $\rho_{\rm s}$ and $\rho_{\rm d}$ are densities corresponding to the accidentally degenerate s and d shells. Using the same trick for accidentally degenerate p and d shells, we manage to obtain a self-consistent solution to the rHF equations, which is necessarily a ground state since the rHF model is convex in the density matrix.

In the following tables, we report the rHF occupied energy levels (in Ha) of all the atoms of the first four rows of the periodic table. In some cases, the Fermi level seems to be an accidentally degenerate eigenvalue:

• •

the 4p and 3d orbitals have the same energy for $z=21,22$;

• •

the 5s and 3d orbitals have the same energy for $23\leq z\leq 26$;

• •

the 5p and 4d orbitals have the same energy for $z=40$;

• •

the 6s and 4d orbitals have the same energy for $z=41,42$;

• •

the 5s and 4d orbitals have the same energy for $z=46,47$.

In all these cases, the occupation number $0\leq n\leq 2$ of the partially occupied d orbitals is also given.

Hydrogen and helium:

[TABLE]

First row:

[TABLE]

Second row:

[TABLE]

Third row:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Fourth row:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Remark 3**.**

Our numerical simulations seem to show that for all $1\leq z\leq 54$, there are no unoccupied negative eigenvalues in the rHF ground states of neutral atoms.

We end this section by the following figures, which back up the conjecture that rHF atomic densities are decreasing radial functions of the distance to the nucleus.

4.1.2 Occupied energy levels in the ${\rm X}\alpha$ model

The tables below provide the negative eigenvalues of the Kohn-Sham X$\alpha$ Hamiltonian (in Ha) for all the atoms of the first four rows of the periodic table . We observe that atoms $z$, with $23\leq z\leq 28$ and $41\leq z\leq 44$ have accidentally degenerate Fermi levels, the degeneracy occurring in all cases between an s-shell and a d-shell (4s-3d for $23\leq z\leq 28$, 5s-4d for $41\leq z\leq 44$) . All the results of this section are obtained for $L_{e}=30$ and $N_{I}$ increasing from $30$ to $75$ as $z$ increases.

Hydrogen and helium:

[TABLE]

First row:

[TABLE]

Second row:

[TABLE]

Third row:

[TABLE]

[TABLE]

[TABLE]

Fourth row:

[TABLE]

[TABLE]

[TABLE]

We end this section by the following figures, which show that as in the ${\rm rHF}$ case, the ${\rm X}\alpha$ atomic densities seem to be decreasing radial functions of the distance to the nucleus.

4.2 Perturbation by a uniform electric field (Stark effect)

In this section, we consider atoms subjected to a uniform electric field, that is to an external potential $\beta W_{\rm Stark}$ with

[TABLE]

or, in spherical coordinates,

[TABLE]

As already mentioned in Section 2.2, $\widetilde{\mathcal{I}}_{z,N}^{\rm rHF/LDA}(\beta W_{\rm Stark})=-\infty$ whenever $\beta\neq 0$, and the corresponding variational problem has no minimizer. However, one can find a minimizer $\gamma_{h}\in{\mathcal{K}}_{N,h}$ to the approximated problem $\widetilde{\mathcal{I}}_{z,N,h}^{\rm rHF/LDA}(\beta W_{\rm Stark})$. Hereafter we consider the carbon atom ($z=6$). Even though the cutoff $m_{h}$ is set equal to $6$, all the terms corresponding to a magnetic number $m>1$ are in fact equal to zero.

The following figures are the plots in the $xy$-plane of the densities $\rho_{h}$ multiplied by $|\mathbf{r}|^{2}$ for the carbon atom ($z=6$) obtained for different values of $\beta$:

For $\beta=10^{-2}$ and $\beta=10^{-1}$, we clearly see spurious boundary effects: part of the electronic cloud is localized in the region where the external potential takes highly negative values. This result is obviously not physical. On the other hand, for the X$\alpha$ model and for $\beta=10^{-3}$ we simply observe a polarization of the electronic cloud. The perturbation potential being not spherically symmetric, it breaks the symmetry of the density. This numerical solution can probably be interpreted as a (nonlinear) resonant state. We will come back to the analysis of this interesting case in a following work.

Fig. 7 shows the amount of electrons of the carbon atom which escape to infinity as a function of the coupling constant $\beta$ (for $L_{e}=100$ and $N_{I}=50$), in the rHF and ${\rm X}\alpha$ case.

In general, the standard ODA is used to achieve convergence (see Section 3). However, for $\beta$ small (resp. large) enough, the occupation numbers are selected as follows: $n_{0,1}^{[n]}=n_{0,2}^{[n]}=2$, $n_{0,3}=2(1-t_{0})$ and $n_{1,1}^{[n]}=2-n_{0,3}=2t_{0}$, $t_{0}$ being the minimizer of

[TABLE]

where

[TABLE]

This modification of ODA significantly increases the rate of convergence for $\beta$ small or large, but does not converge for all intermediate values of $\beta$.

While $\widetilde{\mathcal{I}}_{z,N}^{\rm rHF/LDA}(\beta W_{\rm Stark})=-\infty$ and the corresponding variational problem has no minimizer, the first-order perturbation $\gamma_{z,N,W_{\rm Stark}}^{(1),\rm rHF}$ of the ground state density matrix does exist (see Theorem 4). If we consider the carbon atom, it can be expressed as a function of the unperturbed occupied Kohn-Sham orbitals and of their first-order perturbations. We indeed have

[TABLE]

where ${\mathcal{O}}_{6,6}=\{(0,1),(0,2),(0,3),(1,1)\}$, where $\epsilon_{m,k}^{(0)}$ is the $k$-th lowest eigenvalue of $H_{6,6}^{0,\rm rHF}$ in the subspace ${\mathcal{H}}^{m}$, $\Phi_{m,k}^{(0)}$ an associated normalized eigenfunction and

[TABLE]

while $\epsilon_{m,k}^{(1)}$, $\Phi_{m,k}^{(1)}$ and $n_{m,k}^{(1)}$ satisfy the following self-consistent equation

[TABLE]

We denote by $\epsilon_{m,k,h}^{(0)}$, $\epsilon_{m,k,h}^{(1)}$, $\Phi_{m,k,h}^{(0)}$, $\Phi_{m,k,h}^{(1)}$ and $n_{m,k,h}^{(1)}$ the approximations of $\epsilon_{m,k}^{(0)}$, $\epsilon_{m,k}^{(1)}$, $\Phi_{m,k}^{(0)}$, $\Phi_{m,k}^{(1)}$ and $n_{m,k}^{(1)}$, respectively. For each $(m,k)\in{\mathcal{O}}_{6,6}$, define

[TABLE]

Recall that, $\left(\Phi_{m,k,h}(\beta)\right)_{(m,k)\in{\mathcal{O}}_{6,6}}$ (resp. $n_{m,k,h}^{(1)}(\beta)$) are the eigenfunctions (resp. eigenvalues) of the density matrix $\gamma_{h}$, the minimizer of the approximated problem $\widetilde{\mathcal{I}}_{z,N,h}^{\rm rHF}(\beta W_{\rm Stark})$.

Let $U^{m,k}$ and $\widetilde{U}^{m,k}(\beta)$ be such that

[TABLE]

To show that $\widetilde{\Phi}_{m,k,h}^{(1)}(\beta)\to\Phi_{m,k,h}^{(1)}$ when $\beta\to 0$, it is enough to show that for each $l\geq 0$

[TABLE]

The index $(m,k)$ is omitted for simplicity and the vector $\widetilde{Q}_{l}(\beta)$ is the solution to the linear system

[TABLE]

with

[TABLE]

Our numerical results show that, as expected by symmetry, $n_{m,k,h}^{(1)}=\epsilon_{m,k,h}^{(1)}=0$ for all $(m,k)\in{\mathcal{O}}_{6,6}$, and that the left-hand side of (59) converges to zero linearly in $\beta$ (see Fig.8).

5 Acknowledgements

The authors are grateful to Carlos García-Cervera and Vikram Gavini for valuable discussions.

Appendix: ${\mathbb{P}}4$ radial finite elements

In this appendix, we elaborate on the details of the calculation.

A1. Basis functions

We have chosen the following form functions to build the finite element matrices and tensors:

[TABLE]

[TABLE]

[TABLE]

Their derivatives are given by:

[TABLE]

[TABLE]

Finite element basis:

• •

the 1D Schrödinger equation is solved on the finite interval $[0,L_{e}]$ with zero Dirichlet boundary conditions

• •

the interval $[0,L_{e}]$ is decomposed in $N_{I}$ intervals of positive lengths $h_{1},\cdots h_{N_{I}}$. Let $0=r_{1}<r_{2}<\cdots<r_{N_{I}}<r_{N_{I}+1}=L_{e}$ be such that $h_{k}=r_{k+1}-r_{k}$;

• •

we denote by

[TABLE]

the ${\mathbb{P}}_{4}$ finite element space associated with the so-defined mesh. We have

[TABLE]

• •

we then set for all $1\leq k\leq N_{I}$ and $1\leq j\leq 5$,

[TABLE]

so that $p^{k}_{j}(r_{k}+th_{k})=z_{j}(t)$, and define the basis $(\chi_{1},\cdots,\chi_{4N_{I}-1})$ of $V_{h}$ as follows:

[TABLE]

and for all $2\leq k\leq N_{I}$,

[TABLE]

• •

when considering an atom with Fermi level very close or possible equal to zero, an extra basis function of the form

[TABLE]

is added to the space $V_{h}$. Its derivative is equal to

[TABLE]

A2. Assembling the matrices

Let $\Lambda$ be the bijective mapping from $\{0,1,2,3,4\}$ to $\{1,2,3,4,5\}$ defined by

[TABLE]

Recall that the density is equal to

[TABLE]

Using the finite element basis defined above, one gets that $\rho_{h}(r,\theta)$ is equal to

[TABLE]

where $\kappa$ equal to $1$ (resp. [math]) if the discretization space $V_{h}$ (resp. $V_{h}\cup\{\chi_{4N_{I}}\}$) is considered. In particular, for $0<t_{p,\rm r}<1$ and $-1<t_{q,\theta}<1$, we have that $(t_{p,\rm r}h_{k}+r_{k})^{2}\rho(t_{p,\rm r}h_{k}+r_{k},\arccos(t_{q,\theta}))$ is equal to

[TABLE]

where $P_{l}$ are the Legendre polynomials, which can be calculated using the recurrence relation

[TABLE]

with $P_{0}(x)=1$ and $P_{1}(x)=x$. Note that system (65) is used only to calculate the exchange-correlation potential, thus there is no need to define the density at radius greater than $r_{N_{I+1}}$.

For $\mu(r\mathbf{e})=\frac{\eta^{2}}{4\pi}\frac{e^{-\eta r}}{r}$, then

[TABLE]

Thus the vector $G$ in (45) has one of the following form depending on the discretization space considered

[TABLE]

or

[TABLE]

where

[TABLE]

and

[TABLE]

We denote by

[TABLE]

where $C^{l,m}$, $Q_{l}$ and $F$ are defined by (34), (42) and (44), respectively.

All the matrices $A$, $M_{-2}$, $M_{-1}$, $M_{0}$, $M_{1}$, $V_{\mu}$, $[V_{\rm xc}^{l}]$ and $[\hat{\mathcal{H}}_{l,l^{\prime}}]$ defined in (26), (37), (48) and (66), when considering the discretization space $V_{h}$, are symmetric and have the same pattern:

Their entries can be computed using elementary assembling matrices:

• •

diagonal blocks: for any $1\leq k\leq N_{I}$,

[TABLE]

• •

off-diagonal blocks: for $1\leq k\leq N_{I}-1$

[TABLE]

[TABLE]

When $\chi_{4N_{I}}$ is added to the discretization space, an extra row and an extra column must be added to each of the above matrices. The non-zero additional entries are:

[TABLE]

[TABLE]

[TABLE]

The $y^{k}_{ij}$’s are the entries of the elementary assembling matrices. The latter are defined for the matrices $A$, $M_{-2}$, $M_{-1}$, $M_{0}$, $M_{1}$, $V_{\mu}$, $[V_{\rm xc}^{l}]$ and $[\hat{\mathcal{H}}_{l,l^{\prime}}]$ as follows:

[TABLE]

[TABLE]

in the ${\rm X}\alpha$-case, that is for $v_{\rm xc}(\rho)=-\left(\frac{3}{\pi}\right)^{\frac{1}{3}}\rho^{\frac{1}{3}}$,

[TABLE]

where

[TABLE]

and

[TABLE]

Note that $\rho(t_{p,\rm r}h_{k}+r_{k},\arccos(t_{q,\theta}))$ is calculated with the help of (65). The $y_{\infty}$ are defined for the matrices $A$, $M_{-2}$, $M_{-1}$, $M_{0}$, $V_{\mu}$ and $[\hat{H}_{l,l^{\prime}]}$ as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where $f_{\infty}=\int_{r_{N_{I+1}}}^{\infty}\frac{\chi_{4N_{I}}^{3}(r)}{r}dr=\int_{r_{N_{I+1}}}^{\infty}\frac{(r_{N_{I+1}})^{3}}{r^{4}}dr=\frac{1}{3}$.

In addition to assembling the matrices, we need to deal with the following term

[TABLE]

in order to calculate the right-hand side of (42). Let $k_{n}=1+\mbox{int}(\frac{n}{4})$ and $q_{n}=4-(n\,\mbox{mod}\,4)$, so that $n=4k_{n}-q_{n}$. Then the entries $\displaystyle\sum_{i,j=1}^{N_{h}}F_{ijn}[R_{l}]_{ij}$ of the vector $F:R_{l}$ are computed as follows

[TABLE]

When the base $\chi_{4N_{I}}$ is added, the last entry of the vector $F:R_{l}$ is equal to

[TABLE]

We end this section by providing the values of $\alpha_{ij}$, $\beta_{ij}$ and $\nu_{ij}$, for $1\leq i,j\leq 5$,

[TABLE]

[TABLE]

[TABLE]