Existence and nonexistence of HOMO-LUMO excitations in Kohn-Sham density functional theory

TL;DR

This paper investigates the existence of HOMO-LUMO excitations in Kohn-Sham density functional theory, showing they exist for systems with nuclear charge greater than the number of electrons, but not for neutral atoms under certain conditions.

Contribution

The paper provides a rigorous mathematical analysis of HOMO-LUMO excitation existence in KS-DFT, including conditions under which they do or do not exist, with implications for electronic structure calculations.

Findings

HOMO-LUMO excitations exist for Z>N systems.

Excited states do not exist for neutral atoms Z=N under certain approximations.

Implications for ionization potential estimates from HOMO eigenvalues.

Abstract

In numerical computations of response properties of electronic systems, the standard model is Kohn-Sham density functional theory (KS-DFT). Here we investigate the mathematical status of the simplest class of excitations in KS-DFT, HOMO-LUMO excitations. We show using concentration-compactness arguments that such excitations, i.e. excited states of the Kohn-Sham Hamiltonian, exist for $Z>N$, where $Z$ is the total nuclear charge and $N$ is the number of electrons. The result applies under realistic assumptions on the exchange-correlation functional, which we verify explicitly for the widely used PZ81 and PW92 functionals. By contrast, and somewhat surprisingly, we find using a method of Glaser, Martin, Grosse, and Thirring \cite{glaser1976} that in case of the hydrogen and helium atoms, excited states do not exist in the neutral case $Z=N$ when the self-consistent KS ground state density is replaced by a realistic but easier to analyze approximation (in case of hydrogen, the true Schr\"{o}dinger ground state density). Implications for interpreting minus the HOMO eigenvalue as an approximation to the ionization potential are indicated.

Equations

$$-{\epsilon }_H\approx \text{mbox}ionizationpotential$$

$$v_{ext}(x)=-\sum _{\alpha =1}^M{Z}_{\alpha }\frac{1}{|x-R_{\alpha }|}.$$

$${\mathcal{A}}_{nuc}:=\{\mu \in \mathcal{M}({\Omega }_{nuc}):\mu \ge 0,{\int }_{{\Omega }_{nuc}}\mathrm{d}\mu =Z\},$$

$$v_{ext}(x):=-{\int }_{{\Omega }_{nuc}}\frac{1}{|x-y|}\mathrm{d}\mu (y).$$

$$⟨{\phi }_i,{\phi }_j{⟩}_{L^2}={\int }_{{\mathbb{R}}^3}{\phi }_i(x){\phi }_j(x)\mathrm{d}x={\delta }_{ij}\forall i,j\in \{1,\dots ,n\},$$

$${\mathcal{E}}_{\mu }[\Phi ]=\underset{=:T[\Phi ]}{\underbrace{\sum _{k=1}^{n}2{\int }_{{\mathbb{R}}^3}\frac{1}{2}|\nabla {\phi }_k{|}^2(x)\mathrm{d}x}}+\underset{=:V[\rho ]}{\underbrace{{\int }_{{\mathbb{R}}^3}{v}_{ext}(x)\rho (x)\mathrm{d}x}}+\underset{=:J[\rho ]}{\underbrace{\frac{1}{2}{\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}\frac{\rho (x)\rho (y)}{|x-y|}\mathrm{d}x\mathrm{d}y}}+\underset{=:E_{xc}[\rho ]}{\underbrace{{\int }_{{\mathbb{R}}^3}{e}_{xc}(\rho (x))\mathrm{d}x}},$$

$$\rho (x):=2\sum _{k=1}^n|{\phi }_k(x){|}^2.$$

$$e_{xc}(\rho )=-\frac{3}{4}{\left(\frac{3}{\pi }\right)}^{\frac{1}{3}}{\rho }^{\frac{4}{3}}.$$

$$\Phi \in \mathrm{argmin}{\mathcal{E}}_{\mu }\text{ subject to the constraints }\text{eqref}eq:ortho.$$

$$h_{\mu ,\rho }{\phi }_i:=\left(-\frac{1}{2}\Delta +v_{ext}+v_H+v_{xc}\right){\phi }_i=\sum _{j=1}^N{\lambda }_{ij}{\phi }_j,$$

$$v_H(x)={\int }_{{\mathbb{R}}^3}\frac{1}{|x-y|}\mathrm{d}\mu (y),v_{xc}=\frac{\mathrm{d}}{\mathrm{d}\rho }{e}_{xc}.$$

$$h_{\mu ,\rho }{\phi }_i:=\left(-\frac{1}{2}\Delta +v_{ext}+v_H+v_{xc}\right){\phi }_i={\epsilon }_i{\phi }_i.$$

$$\big{(}\varphi_{1},\ldots,\varphi_{n-1},{\color[rgb]{0,0,1}\varphi_{n}}\big{)}\enskip\longrightarrow\enskip\big{(}\varphi_{1},\ldots,\varphi_{n-1},{\color[rgb]{1,0,0}\varphi_{n+1}}\big{)},$$

$$\mathcal{E}_{\mu,\rho}[\chi]=\langle\chi,h_{\mu,\rho}\chi\rangle=\frac{1}{2}\int_{\mathbb{R}^{3}}|\nabla\chi|^{2}+\int_{\mathbb{R}^{3}}\big{(}v_{ext}+v_{H}+v_{xc}(\rho)\big{)}|\chi|^{2}.$$

$${\phi }_H\in \mathrm{argmax}{\mathcal{E}}_{\mu ,\rho }\text{ subject to the constraints }{\phi }_H\in \mathrm{Span}\{{\phi }_1,\dots ,{\phi }_n\},⟨{\phi }_H,{\phi }_H⟩=1$$

$${\phi }_L\in \mathrm{argmin}{\mathcal{E}}_{\mu ,\rho }\text{ subject to the constraints }⟨{\phi }_i,{\phi }_L⟩=0,\forall i\in \{1,\dots ,n\},⟨{\phi }_L,{\phi }_L⟩=1.$$

$$h_{\mu ,\rho }{\phi }_H={\epsilon }_H{\phi }_H,h_{\mu ,\rho }{\phi }_L={\epsilon }_L{\phi }_L,$$

$$e_{xc}(0)=0=v_{xc}(0),v_{xc}\le 0,|v_{xc}(\rho )|\le c_{xc}(1+{\rho }^{p-1})$$

$$\text{There exists }q\in\big{[}1,\tfrac{3}{2}\big{)}\text{ such that }\limsup_{\rho\to 0+}\frac{e_{xc}(\rho)}{\rho^{q}}<0.$$

$$\mathcal{A}=\{({\phi }_1,\dots ,{\phi }_n)\in H^1({\mathbb{R}}^3{)}^n:⟨{\phi }_i,{\phi }_j⟩={\delta }_{ij}\}$$

$$\mathcal{A}^{H}_{\Phi}=\big{\{}\varphi_{H}\in\mathrm{Span}\{\varphi_{1},\ldots,\varphi_{n}\}:||\varphi_{H}||_{2}=1\big{\}}$$

$$\mathcal{A}_{\Phi}^{L}=\big{\{}\varphi_{L}\in H^{1}(\mathbb{R}^{3}):\langle\varphi_{k},\varphi_{L}\rangle=0\text{ for }k\in\{1,\ldots,n\},||\varphi_{L}||_{2}=1\big{\}}.$$

$$\Phi \in \underset{\mathcal{A}}{\mathrm{argmin}}{\mathcal{E}}_{\mu },{\phi }_H\in \underset{{\mathcal{A}}_{\Phi }^H}{\mathrm{argmax}}{\mathcal{E}}_{\mu ,\rho },{\phi }_L\in \underset{{\mathcal{A}}_{\Phi }^L}{\mathrm{argmin}}{\mathcal{E}}_{\mu ,\rho }.$$

$$\Phi \mapsto T[\Phi ]\text{ is continuous and weakly lower semicontinuous on }{H}^1({\mathbb{R}}^3{)}^n.$$

$$(\Phi,\mu)\mapsto V[\rho]\text{ is strong $\times$ weak${}^{*}$ continuous on }\big{(}L^{4}\cap L^{2}\big{)}\times\mathcal{M}.$$

$$\Phi\mapsto J_{H}[\rho]\text{ continuous on }\big{(}L^{\nicefrac{{12}}{{5}}}(\mathbb{R}^{3})\big{)}^{n}.$$

$$\Phi\mapsto E_{xc}[\rho]\text{ is continuous on }\big{(}L^{2p}(\mathbb{R}^{3})\big{)}^{n}.$$

$$V[\rho]=-\int_{\mathbb{R}^{3}}\left(\frac{1}{|\cdot|}\ast\rho\right)\,\mathrm{d}\mu\geq-||\mu||_{\mathcal{M}}||\tfrac{1}{|\cdot|}\ast\rho||_{\infty}\geq-||\mu||_{\mathcal{M}}\big{(}||v_{2}||_{2}||\rho||_{2}+||v_{\infty}||_{\infty}||\rho||_{1}\big{)}$$

$$||\rho |{|}_p\le ||\rho |{|}_q^{\theta }||\rho |{|}_r^{1-\theta }\text{ with }q\le p\le r,\frac{1}{p}=\frac{\theta }{q}+\frac{1-\theta }{r}$$

$$|e_{xc}(\rho)|=\big{|}e_{xc}(0)+\int_{0}^{\rho}v_{xc}(\xi)\,\mathrm{d}\xi\big{|}\leq c_{xc}\big{(}\rho+\tfrac{1}{p-1}\rho^{p}\big{)}$$

Full text

1. 1 Introduction

2. 2 Mathematical setting

3. 2.1 Kohn-Sham equations

4. 2.2 Excitations

5. 3 Existence of HOMO-LUMO excitations

6. 3.1 Assumptions

7. 3.2 Existence proof

8. 3.3 Higher excitations

9. 4 Optimal HOMO-LUMO excitations

10. 5 Nonexistence of HOMO-LUMO excitations in the neutral case $Z=N$

11. 6 Appendix: Analytical properties of PZ81 and PW92

12. Perdew-Zunger (PZ81)

1 Introduction

Electronic excitations play an important role in the description of molecular properties such as absorption spectra, photoexcitation, state-to-state transition probabilities, reactivity, charge transfer processes, and reaction kinetics [5, 6, 12]. In numerical computations of these response properties, the standard model is Kohn-Sham density functional theory (KS-DFT), because of its good compromise between accuracy and feasibility for large systems (see [19] for a textbook account and [3] for a recent review). It is then of interest to investigate the mathematical status of excitations in KS-DFT.

In this paper we analyze the simplest such excitations, HOMO-LUMO transitions, in the setting of the local density approximation (LDA). For a systematic comparison of HOMO-LUMO excitations with experimental data see e.g. [2, 23]. Even in this case we are not aware of previous rigorous results. Our findings are the following (see Figure 1).

For positively charged systems (i.e., total nuclear charge $Z$ greater than the number $N$ of electrons) such excitations – mathematically, excited states of the KS Hamiltonian – are rigorously proven to exist, under realistic assumptions on the exchange-correlation functional which we verify explicitly for the widely used PZ81 and PW92 functionals. See Theorem 1 in Section 3. As a corollary we also establish existence of optimal excitations with respect to suitable control goals recently introduced in [10], without requiring the simplifying assumption in [10] of bounded domains. See Section 4.

By contrast, the neutral case $Z=N$ holds a surprise. In the case of the hydrogen and helium atoms, we prove that excited states do not exist when the self-consistent KS ground state density is replaced by a realistic but easier to analyze closed-form approximation (in case of hydrogen, the true Schrödinger ground state density). See Theorem 4 in Section 5.

Mathematically, the existence result relies on concentration-compactness arguments, and should not come as a surprise to experts. The nonexistence result uses a not widely known method by Glaser, Martin, Grosse, and Thirring (GMGT) [11]. The latter method could, in principle, also be applied to numerical KS ground state densities; we expect that for some atoms and molecules, including hydrogen and helium, the GMGT nonexistence criterion (that a certain integral associated with the effective KS potential lies below a threshold value) would be satisfied.

Physically, these results indicate a significant artefact of KS-DFT. In the full $N$-electron Schrödinger equation, neutral systems (and even systems with $Z>N-1$) are known to possess infinitely many excited states below the bottom of the continuous spectrum. This is a celebrated result by Zhislin [24]; for a modern variational proof see [9]. The analogous result also holds in Hartree-Fock theory: for $Z>N$ the Fock operator associated with the Hartree-Fock ground state density possess infinitely many bound states below the continuous spectrum [18, Lemma II.3], the latter being the interval $[0,\infty)$. (It is also known [15] that the Hartree-Fock enegery functional possess infinitely many critical points below 0.) Our results suggest that in KS-DFT, the threshold for existence of infinitely many excited states is shifted from $Z>N-1$ to $Z>N$. This is a previously unnoticed but important qualitative consequence of the (well known) incomplete cancellation of the self-interaction energy in KS-DFT.

It is interesting to interpret the nonexistence of excitations from the point of view of numerical computations in finite basis sets, or mathematical analysis (as in [10]) in bounded domains. Consider a neutral system for which (exact) excitations do not exist. In a finite basis set, or a bounded domain, the spectrum of the KS Hamiltonian is purely discrete and therefore excited states exist. In the limit as the basis set approaches completeness, or the domain approaches the whole of $\mathbb{R}^{3}$,

(i) the LUMO energy $\varepsilon_{L}$ (i.e., the lowest unoccupied eigenvalue of the KS Hamiltonian) will remain well-defined, and approaches the bottom of the continuous spectrum (which equals [math], see Theorem 2 and Theorem 4)

(ii) the LUMO (i.e., the lowest unoccupied eigenstate) will become more and more delocalized, failing to converge to a bound state.

Thus in contrast to common (explicit or implicit) belief, restriction to finite basis sets or bounded domains may be not just a negligible technicality, but significantly alters the physical nature of LUMO excitations, from stable bound state (i.e., invariant under the dynamics of the KS ground state Hamiltonian) to a delocalized, dispersing state associated with the continuous spectrum.

(ii) makes it very tempting to physically interpret the HOMO-LUMO excitation in the nonexistence case as an (approximation to an) ionization process. This interpretation together with (i) yields ionization potential $\approx\varepsilon_{L}-\varepsilon_{H}=0-\varepsilon_{H}$ (where $\varepsilon_{H}$ is the HOMO energy, i.e. the highest occupied eigenvalue of the KS Hamiltonian), lending new theoretical support to the famous semi-empirical formula

[TABLE]

which often agrees quite well with experimental data [2, 23].

2 Mathematical setting

We start by recalling well-known mathematical facts about Kohn-Sham density functional theory (KS-DFT) [14, 13, 19]. Readers familiar with these facts might want to skip this part. After that we give a variational definition of HOMO-LUMO excitations as introduced recently in [10], which works irrespective of degeneracies and is convenient for the mathematical analysis of excitations.

2.1 Kohn-Sham equations

We consider a system of $N$ non-relativistic electrons in $\mathbb{R}^{3}$ in the electrostatic potential generated by $M$ nuclei of charges $Z_{1},\ldots,Z_{M}$ located at positions $R_{1},\ldots,R_{M}\in\mathbb{R}^{3}$,

[TABLE]

In fact, for our analysis it is not essential that the nuclei are point particles. It suffices to assume the nuclear charge distribution is given by a nonnegative Radon measure $\mu$ with total mass $Z>0$ supported on a compact set $\Omega_{nuc}\subseteq\mathbb{R}^{3}$, i.e. we consider any $\mu$ belonging to

[TABLE]

where $\mathcal{M}(\Omega_{nuc})$ denotes the space of signed Radon measures on $\Omega_{nuc}$, and

[TABLE]

For simplicity we look at a spin-unpolarized system, so the number $N$ of electrons is even, i.e. $N=2n$ for some $n\in\mathbb{N}$. In this case Kohn-Sham DFT describes the electrons by $n$ orbitals $\varphi_{1},\ldots,\varphi_{n}:\mathbb{R}^{3}\to\mathbb{C}$, each occupied by two electrons of opposite spin. They are $L^{2}$-orthonormal, i.e.

[TABLE]

and we denote $\Phi:=(\varphi_{1},\ldots,\varphi_{n})$. Note that in the following $\langle\cdot,\cdot\rangle$ will always denote the $L^{2}$ inner product. Then the corresponding Kohn-Sham energy functional is given by

[TABLE]

where $e_{xc}$ gives the exchange-correlation energy per unit volume and $\rho$ is the total electron density, that is,

[TABLE]

The KS energy hence consists of the following terms: $T$ is the kinetic energy of the electrons, $V$ is the potential energy from the electron-nuclei interaction with $v_{ext}$ being the electrostatic potential of the nuclei (2.3); $J_{H}$ (the Hartree energy) describes the energy corresponding to the interelectron repulsion if the electrons were mutually independent; $E_{xc}$ is the exchange-correlation energy which accounts for correlation effects correcting the simple independent ansatz of $J_{H}$.

Note here that the Coulomb potential over the whole $\mathbb{R}^{3}$ is not in any $L^{p}$-space, but it is in $L^{2}(\mathbb{R}^{3})+L^{\infty}(\mathbb{R}^{3})$ and we will be using the splitting $\tfrac{1}{|\cdot|}=v_{2}+v_{\infty}$, where $v_{2},v_{\infty}$ lie in $L^{2},L^{\infty}$, respectively.

Precise assumptions on $e_{xc}$ which are sufficient for our mathematical results and cover standard local density approximation (LDA) exchange-correlation functionals used in practice are given in Section 3.1. A basic example derived from the homogeneous electron gas is the Dirac exchange energy

[TABLE]

The ground state of the system is given by

[TABLE]

Any ground state $\Phi:=(\varphi_{1},\ldots,\varphi_{N})$ satisfies the Euler-Lagrange equations of the system, the Kohn-Sham equations

[TABLE]

where the Lagrange multipliers $\lambda_{ij}$ arise due to the orthonormality condition (2.4). The Hartree and exchange-correlation potentials are given by

[TABLE]

Since the effective one-body operator $h_{\mu,\rho}$ in (2.9) (the Kohn-Sham Hamiltonian) is invariant under unitary transformations, the KS equations can be brought into their canonical form

[TABLE]

2.2 Excitations

Following [10] we confine ourselves here to the simplest model for electronic excitations, the HOMO-LUMO transition. In this transition an electron pair migrates from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO). For the KS-orbitals $\Phi=(\varphi_{1},\ldots,\varphi_{n})$ ordered by the size of their eigenvalue in (2.11) this means

[TABLE]

where $\varphi_{n}$ is the HOMO and $\varphi_{n+1}$ - the eigenstate corresponding to the next higher eigenvalue of $h_{\mu,\rho}$ - is the LUMO.

To define HOMO and LUMO in a variational way, we consider the excitation energy functional [10] given by the quadratic form associated with KS Hamiltonian $h_{\mu,\rho}$ (2.11),

[TABLE]

Now define a HOMO $\varphi_{H}$ by

[TABLE]

and a LUMO $\varphi_{L}$ by

[TABLE]

If they exist, HOMO and LUMO clearly satisfy the KS equations

[TABLE]

for some eigenvalue $\varepsilon_{H}$ (the HOMO energy) and $\varepsilon_{L}$ (the LUMO energy).

3 Existence of HOMO-LUMO excitations

In this section we show that for positively charged systems $(Z>N)$ there always exist HOMO-LUMO excitations. This generalizes a corresponding result in [10] to unbounded domains, except that in bounded domains no restriction on $Z$ are needed. In the latter case such a result is not straightforward due to the possibility of “mass escaping to infinity”, and requires concentration-compactness arguments [16, 17]. The reader may wonder whether the assumption $Z>N$, which is essential in our proof, is really necessary. The authors of course asked themselves the same question. For counterexamples to existence in the case $Z=N$ see Section 5.

3.1 Assumptions

Assumptions on the exchange-correlation energy

We assume that $e_{xc}:[0,\infty)\to\mathbb{R}\text{ is continuously differentiable }$ with

[TABLE]

for some exponent $p$ with $p\in[1,\frac{5}{3})$ and constant $c_{xc}>0$.

Furthermore we need as in [1]:

[TABLE]

Remark**.**

These assumptions are trivially satisfied by the Dirac exchange (2.7) with $p=q=\tfrac{4}{3}$. In our appendix we check explicitly that these assumptions are also satisfied for the two most popular LDA exchange-correlation functionals: Perdew-Zunger (PZ81) [21] and Perdew-Wang (PW92) [20]

Admissible sets of orbitals In order to write these definitions in a more compact way we introduce the following sets: The KS admissible set

[TABLE]

the HOMO admissible set

[TABLE]

and the LUMO admissible set

[TABLE]

The governing variational principles for the occupied KS orbitals, HOMO and LUMO can now be summarized as

[TABLE]

We start with some estimates for the KS energy functional (2.5).

Lemma 1** (Lower bounds on KS energy functional).**

The terms in the KS energy functionals have the following properties

$T[\Phi]\geq\tfrac{1}{2}T[\Phi]+\tfrac{1}{c}||\rho||_{3}$* and*

[TABLE] 2. 2.

V[\rho]\geq-||\mu||_{\mathcal{M}}\big{(}||v_{\infty}||_{\infty}||\rho||_{1}+||v_{2}||_{2}||\rho||^{\nicefrac{{1}}{{4}}}_{1}||\rho||_{3}^{\nicefrac{{3}}{{4}}}\big{)}* and*

[TABLE] 3. 3.

$J_{H}\geq 0$* and*

[TABLE] 4. 4.

E_{xc}[\rho]\geq-c_{xc}\big{(}||\rho||_{1}+\frac{1}{p-1}||\rho||_{1}^{\nicefrac{{(3-p)}}{{2}}}||\rho||_{3}^{\nicefrac{{3(p-1)}}{{2}}}\big{)},* where $p\in(1,\tfrac{5}{3})$ is the exponent from assumption (3.1), and*

[TABLE]

Proof.

The first inequality is a standard result in DFT, but we include it for the sake of completeness. By a well-known result, see e.g. [4], we have $T[\Phi]\geq||\nabla\sqrt{\rho}||^{2}_{2}$, so the inequality follows by applying the Sobolev embedding $H^{1}\hookrightarrow L^{6}$ to the function $u=\sqrt{\rho}$.

Estimate 2 follows from the duality between $\mathcal{M}(\Omega_{nuc})$ and $C_{b}(\Omega_{nuc})$ and then Cauchy-Schwarz

[TABLE]

and finally bounding the $L^{2}$-norm of $\rho$ by the Hölder interpolation inequality,

[TABLE]

with $p=2,q=1,r=3$.

The positivity of $J_{H}$ is trivial.

Ad 4: By our assumtpion on $v_{xc}$ we have

[TABLE]

The estimate again follows from Hölder interpolation (3.4) with $q=1,r=3,\theta=\tfrac{(3-p)}{2p}$. In all four cases the continuity results follow by pointwise continuity of the integrand and the proven bounds. ∎

The second Lemma considers the excitation functional (2.13). In the following we denote $\rho_{\chi}=|\chi|^{2}$.

Lemma 2** (Lower bounds on excitation functional).**

The terms in the excitation functional have the following properties

$T[\psi]\geq\tfrac{1}{2}T[\psi]+\frac{1}{c}||\rho_{\psi}||_{3}$* and*

[TABLE] 2. 2.

\int_{\mathbb{R}^{3}}v_{ext}\rho_{\psi}\geq-||\mu||_{\mathcal{M}}\big{(}||v_{\infty}||_{\infty}||\rho_{\psi}||_{1}+||v_{2}||_{2}||\rho_{\psi}||^{\nicefrac{{1}}{{4}}}_{1}||\rho_{\psi}||_{3}^{\nicefrac{{3}}{{4}}}\big{)}* and*

[TABLE] 3. 3.

$\int_{\mathbb{R}^{3}}(\tfrac{1}{|\cdot|}\ast\rho)\rho_{\psi}\geq 0$* and*

[TABLE] 4. 4.

\int_{\mathbb{R}^{3}}v_{xc}(\rho)\rho_{\psi}\geq-c_{xc}\big{(}||\rho_{\psi}||_{1}+||\rho||_{p}^{p-1}||\rho_{\psi}||_{1}^{\nicefrac{{(3-p)}}{{(2p)}}}||\rho_{\psi}||_{3}^{\nicefrac{{3(p-1)}}{{(2p)}}}\big{)}* where $p\in[1,\tfrac{5}{3})$ is again the exponent from our assumptions on $v_{xc}$ and*

[TABLE]

In particular, the map $(\Phi,\chi,\mu)\mapsto\mathcal{E}_{\mu,\rho}[\chi]$ is weak $\times$ strong $\times$ weak∗ continuous and weak $\times$ weak $\times$ weak∗ lower semicontinuous on $H^{1}(\mathbb{R}^{3})^{n}\times H^{1}(\mathbb{R}^{3})\times\mathcal{M}$.

Proof.

Statements 1-3 follow by the same line of reasoning as in Lemma 2 above. The fourth assertion follows by the same argument given in [10], but we include it for the sake of completeness. By our assumption on the exchange-correlation potential $v_{xc}$ we have

[TABLE]

The asserted bound now follows from the Hölder interpolation inequality (3.4) with $q=1,r=3,\theta=\frac{3-p}{2p}$. The continuity follows from the pointwise continuity of $\rho\mapsto v_{xc}(\rho)$ together with the bounds (3.1). ∎

Remark**.**

Our existence results do not require $\Phi$ to be the KS ground state, but only to satisfy a fast enough decay property, like

[TABLE]

By the results of [1], whenever $Z\geq N=2n$ and $e_{xc}$ satisfies our assumptions (3.1), then there exists a KS ground state $\Phi=(\varphi_{1},\ldots,\varphi_{n})$ and any such ground state fulfills the additional property (3.5).

Theorem 1** (Existence of HOMO-LUMO excitations).**

For any admissible nuclear charge distribution $\mu\in\mathcal{A}_{nuc}$ and any set of orbitals $\Phi=(\varphi_{1},\ldots,\varphi_{n})\in\mathcal{A}$ the excitation functional possesses a maximizer $\varphi_{H}$ on $\mathcal{A}_{\Phi}^{H}$ (i.e., a HOMO).

If additionally we have $Z>N$, i.e. a positively charged system, and $\Phi\in\mathcal{A}$ satisfies the decay property (3.5), then there exists also a minimizer $\varphi_{L}$ on $\mathcal{A}_{\Phi}^{L}$ (i.e., a LUMO).

Proof.

Existence of a HOMO is elementary since $\chi\mapsto\mathcal{E}_{\mu,\rho}[\chi]$ is continuous on $H^{1}(\mathbb{R}^{3})$ (see Lemma 2) and the admissible set $A_{\Phi}^{H}$ is compact since it is a closed and bounded subset of a finite dimensional space. The case of a LUMO is more difficult and will be treated in the next section. ∎

3.2 Existence proof

For a given $\Phi=(\varphi_{1},\ldots,\varphi_{n})$, define $\gamma_{\Phi}:=\sum_{k=1}^{n}|\varphi_{k}\rangle\langle\varphi_{k}|$, the projection onto the span of the $\varphi_{k}$. Then we define the translation invariant excitation functional by

[TABLE]

Note $\mathcal{E}$ depends on $\mu,\Phi$ but to keep notation simple and distinguish $\mathcal{E}$ from $\mathcal{E}_{\mu,\rho}$ we do not explicitly stress this via the notation.

Now $\mathcal{E}$ is again defined over the whole $H^{1}(\mathbb{R}^{3})$ and if we can prove the existence of a minimizer $\hat{\chi}$ with

[TABLE]

then we have the existence of a LUMO given by $\varphi_{L}:=(Id-\gamma_{\Phi})\hat{\chi}$.

As usual in a concentration compactness argument we embed our problem in a one-parameter family of minimization problems

[TABLE]

parametrized by $\lambda\in\mathbb{R}_{+}$.

We can now state the first central Lemma of this section:

Lemma 3** (Fundamental properties of the infimum).**

The map given by $\lambda\mapsto I_{\lambda}$ has the following properties

1. (i)

$I_{0}=0$* and $-\infty<I_{\lambda}<0$ for every $\lambda>0$,* 2. (ii)

$\lambda\mapsto I_{\lambda}$* is continuous and strictly decreasing, i.e. for $0<\eta<\lambda$ the inequality $I_{\lambda}<I_{\eta}$ holds*

Proof.

The fact that $I_{0}=0$ is clear and since $(Id-\gamma_{\Phi})$ is a bounded operator from $H^{1}$ into itself $I_{\lambda}>-\infty$ follows by the bounds in Lemma 2.

Usually the strict subadditivity condition is the hardest part in a concentration compactness proof, but fortunately here it is quite easy.

The only thing we have to prove is actually $I_{\lambda}<0$ for $\lambda>0$. Since by the structure of our problem we have $\theta I_{\eta}=I_{\theta\eta}$ for all $\theta,\eta\geq 0$. Hence, for $\lambda>\eta$ we choose $\theta>1$ such that $\theta\eta=\lambda$ and obtain

[TABLE]

Also continuity is clear since $\lambda\mapsto I_{\lambda}$ is just a linear function. The inequality $I_{\lambda}<0$ is also important because otherwise we would be considering the constant zero-function.

$\boxed{I_{\lambda}<0}$ Take a radially symmetric testfunction $\psi\in C^{\infty}_{c}(\mathbb{R}^{3})$ with supp$\psi\subseteq B_{K}^{c}$ for some radius $K>0$, $||\psi||=1$. Here $B_{K}^{c}$ denotes the complement of the ball $B_{K}$ around the origin of radius $K>0$. Then define

[TABLE]

Note that

[TABLE]

where we used in the second to last equality that supp$\psi\subseteq B_{K}^{c}$ and in the last step the exponential decay of the KS-orbitals together with the estimate

[TABLE]

Hence $\langle\varphi_{k},\psi_{\sigma}\rangle$ decays exponentially for $\sigma\to 0$, so it can be neglected up to higher orders, i.e.

[TABLE]

With this at hand we can estimate the energy of $\psi_{\sigma}$

[TABLE]

Since we want to estimate the energy from above we do not need to consider the exchange-correlation term since it only gives a negative contribution. The kinetic energy is easily estimated

[TABLE]

Our next step will be to estimate the Hartree and the nuclei term and bring them in a similar form

[TABLE]

where we used the radial symmetry of $\psi$ in the last step.

Exactly the same steps transform the external potential term into

[TABLE]

Since the support of $\mu$ is compact and again due to supp$\psi\subseteq B_{K}^{c}$, we have for $\sigma$ small enough in the last integral $\max\{\sigma|x|,|z+\sqrt{\sigma}\hat{\mathrm{e}}|\}=|z+\sqrt{\sigma}\hat{\mathrm{e}}|$.

Putting everything together we obtain for $\sigma$ small enough

[TABLE]

In this last inequality the assumption of a positively charged system $Z>N$ is crucial.

∎

Our last Lemma is just a standard lower semicontinuity result.

Lemma 4**.**

Let $(\chi_{n})_{n\geq 1}$ be bounded in $H^{1}(\mathbb{R}^{3})$ converging weakly to $\chi$. If $||\chi||_{2}=\lim\limits_{n\to\infty}||\chi_{n}||_{2},$ then $\chi_{n}\to\chi$ strongly in $L^{p}(\mathbb{R}^{3})$ for $p\in[2,6)$ and furthermore

[TABLE]

Proof.

The weak convergence in $H^{1}$ together with the convergence of the $L^{2}$-norms gives strong convergence in $L^{2}$ which in turn gives strong convergence in $L^{p}$, $p\in[2,6)$ due to Hölder interpolation inequality and the fact that $(\chi_{n})_{n\geq 1}$ is bounded in $L^{q}$ for $q\in[2,6]$ by Sobolev embedding.

The inequality follows by weak lower semicontinuity of $\chi\mapsto\mathcal{E}_{\mu,\rho}[\chi]$ established in Lemma 2 and the fact that $(Id-\gamma_{\Phi})$ is a continuous operator from $H^{1}$ into itself. ∎

With all these Lemmata at our disposal, we are now able to prove the existence of a minimizer to (3.7) and hence the existence of a LUMO.

Take a minimizing sequence $(\chi_{n})_{n}$ for $I_{\lambda}$. Then w.l.o.g. $\chi_{n}=(Id-\gamma_{\Phi}){\chi}_{n}$, i.e. $\chi_{n}$ lies in the orthorgonal complement of $\mathrm{span}\{\varphi_{1},\ldots,\varphi_{n}\}$. Therefore the sequence is bounded in $H^{1}(\mathbb{R}^{3})$ by Lemma 2 and hence (up to taking a subsequence) $\chi_{n}\rightharpoonup\chi$ weakly in $H^{1}$ for some $\chi\in H^{1}$.

If $||\chi||_{2}=\lim\limits_{n\to\infty}||\chi_{n}||_{2}$, then by Lemma 4 we get

[TABLE]

and since $\chi$ also lies in the (convex and closed, hence weakly closed) orthorgonal complement of $\mathrm{Span}\{\varphi_{1},\ldots,\varphi_{n}\}$ we have $||(Id-\gamma_{\Phi})\chi||=||\chi||=\lambda$ so $\chi$ is a minimizer.

The rest of the proof consists in ruling out $||\chi||<\lim\limits_{n\to\infty}||\chi_{n}||=:\Lambda$.

For this purpose, assume first $||\chi||=0$. Then $\chi_{n}\to 0$ and $(Id-\gamma_{\Phi})\chi_{n}\to 0$ in $L^{p}_{loc}(\mathbb{R}^{3})$ for $p\in[2,6)$. In particular,

[TABLE]

for any fixed $K>0$.

The first term converges to 0 and the second one can be made arbitrarily small because $v_{ext}$ vanishes at infinity and $v_{xc}(\rho)$ also vanishes since $\rho$ decays exponentially.

So in conclusion we would get $I_{\lambda}\geq 0$, a contradiction to Lemma 3.

So we only have to rule out $\alpha:=||\chi||\in(0,\Lambda)$. Here we use a standard argument. Following e.g. [8], define a quadratic partition of unity $\xi^{2}+\zeta^{2}=1$ with $\xi$ being smooth, radially non-increasing and with

[TABLE]

Define $\xi_{R}:=\xi(\tfrac{\cdot}{R})$. Then for all $n\in\mathbb{N},$ the map $R\mapsto||\chi_{n}\xi_{R}||$ is continuous, non-decreaing, vanishes for $R\to 0$, and converges to $\lambda$ for $R\to\infty$.

Hence we can choose a sequence $(R_{n})_{n\geq 1}$ such that $||\chi_{n}\xi_{R_{n}}||=\alpha$. Furthermore we claim that this sequence goes to infinity as $R_{n}\to\infty$.

Assume for contradiction that this were not the case. Then up to a subsequence we have $R_{n}\to R_{*}$ for some $R_{*}\geq 0$. But this would give us

[TABLE]

which is not possible since $\xi^{2}_{R_{*}}<1$ for $x\neq 0$. Hence $R_{n}\to\infty$.

Next we define $\chi_{1,n}:=\chi_{n}\xi_{R_{n}}$ and $\chi_{2,n}:=\chi_{n}\zeta_{R_{n}}$ such that $|\chi|^{2}=|\chi_{1,n}|^{2}+|\chi_{2,n}|^{2}$ and (by definition of the $R_{n}$) $||\chi_{1,n}||=\alpha,$ $||\chi_{2,n}||=||\chi_{n}||^{2}-\alpha$.

Note that those two sequences are also bounded in $H^{1}(\mathbb{R}^{3})$, since $\chi_{n}$ is bounded in $H^{1}$ and $\xi_{R_{n}}$ belongs to $W^{1,\infty}$.

Furthermore the cutoff gives us that $(\chi_{1,n})_{n}$ converges weakly to $\chi$ and hence as a by-product $\chi_{2,n}$ converges to 0 in $L^{p}_{loc}$. In order to see this, we take a $\Psi\in C^{\infty}_{c}$ and consider

[TABLE]

where we used that $\chi_{n}$ converges weakly to $\chi$ and the fact that $\xi_{R_{n}}$ converges pointwise to 1 and stays uniformly bounded. So we proved weak convergence on a dense subset and by boundedness we conclude.

Now we can apply Lemma 4 to the sequence $\chi_{1,n}$ since by construction $||\chi_{1,n}||=\alpha=||\chi||$, and obtain $\mathcal{E}[\chi]\leq\liminf\limits_{n\to\infty}\mathcal{E}[\chi_{1,n}]$ and $\chi_{1,n}\to\chi$ in $L^{p}$ for $p\in[2,6)$.

Our next step is proving

[TABLE]

where $o(1)$ is the usual Landau notation meaning terms which vanish for $n\to\infty$.

First note that $\chi_{1,n}\rightharpoonup\chi$ and also $\langle\varphi_{k},\chi\rangle=0$, hence

[TABLE]

Even more, remembering $R_{n}\to\infty$ for $n\to\infty$, we have $|\nabla\xi_{R_{n}}|\leq\tfrac{c}{R_{n}}=o(1)$. Therefore we obtain

[TABLE]

This already gives us $T[\chi_{1.n}]\leq T[\chi_{n}]+o(1)$.

For the potential energy we start by proving $V[\chi_{n}]=V[\chi_{1,n}]+V[\chi_{2,n}]+o(1)$.

Due to the fact that $(\chi_{2,n})_{n}$ is bounded and converges pointwise to zero we have $\chi_{2,n}\rightharpoonup 0$ in $H^{1}(\mathbb{R}^{3})$ and hence $\langle\varphi_{k},\chi_{2,n}\rangle=o(1)$, which gives us $(Id-\gamma_{\Phi})\chi_{2,n}=\chi_{2,n}+o(1)$. Consequently

[TABLE]

So up until now we have

[TABLE]

The only fact that remains to be proven is $V[\chi_{2,n}]\geq o(1)$. This follows from the standard argument that the only negative contributions come from $v_{ext}$ and $v_{xc}$, i.e.

[TABLE]

where the first term vanishes because $\chi_{2,n}\to 0$ in $L^{p}_{loc}$, and the second one also tends to [math] because $K>0$ was arbitrary and $v_{xc}$ and $v_{ext}$ decay at infinity.

In conclusion we proved $\mathcal{E}[\chi_{n}]\geq\mathcal{E}[\chi_{1,n}]+o(1)$ which gives us after passing to the limit

[TABLE]

where we used in the last inequality that $\chi$ lies in the orthogonal complement of the occupied orbitals and hence $\alpha=||\chi||=||(Id-\gamma_{\Phi})\chi||$ so $\chi\in K_{\alpha}$.

But since we assumed $\alpha<\lambda$, this contradicts statement $(ii)$ in Lemma 3. This completes the proof of Theorem 1.

3.3 Higher excitations

We now complete the rigorous justification of the picture on the left-hand side of Figure 1, i.e. we fully characterize the spectrum of the KS Hamiltonian in the positively charged case.

Theorem 2** (Spectrum of KS Hamiltonian in the case $Z>N$).**

Consider $Z>N$, i.e. a positively charged system. Then the KS Hamiltonian $h_{\mu,\rho}$ given in (2.11) has infinitely many negative eigenvalues of finite multiplicity below the bottom of its essential spectrum,

[TABLE]

Proof.

Due to assumption (3.1) the potential $v_{ext}+v_{H}+v_{xc}$ belongs to the space $L^{2}(\mathbb{R}^{3})+L^{\infty}_{\varepsilon}(\mathbb{R}^{3})$, consisting of potentials $v$ which for any given $\varepsilon>0$, can be decomposed as $v=v_{1}+v_{2}$ with $||v_{1}||_{2}<\infty$, $||v_{2}||_{\infty}<\varepsilon$. Therefore is relatively compact w.r.t. the Laplacian. Hence by Weyl’s Theorem (see e.g. Chapter XIII.4 of [22]) $h_{\mu,\rho}$ is self-adjoint and we have \sigma_{ess}\big{(}h_{\mu,\rho}\big{)}=\sigma_{ess}\big{(}-\Delta\big{)}=[0,\infty). Furthermore, $h_{\mu,\rho}$ is bounded from below. To see this, write $h_{\mu,\rho}=-\Delta+U_{2}+U_{\infty}$ with $U_{2}\in L^{2}(\mathbb{R}^{3})$ and $U_{\infty}\in L^{\infty}_{\varepsilon}(\mathbb{R}^{3})$. Then for every $\psi\in H^{1}(\mathbb{R}^{3})$ with $||\psi||_{2}=1$ we have

[TABLE]

Now the Sobolev embedding $||\varphi||_{4}^{2}\leq C||\nabla\varphi||_{2}^{\nicefrac{{3}}{{2}}}$ and Young’s inequality give the asserted lower bound. Now note that

[TABLE]

and since by assumption $Z>N$ we can apply Lemma II.1 of [18] which gives us that the right-hand side operator of (3.10) is negative on an infinite-dimensional subspace. Hence so is $h_{\mu,\rho}$. The min-max principle now gives that it has infinitely many negative eigenvalues. This completes the proof. ∎

4 Optimal HOMO-LUMO excitations

In [10], motivated by the design of photovoltaic materials, results are given for optimal HOMO-LUMO excitations with respect to some control goals. The control is the nuclear charge distribution (e.g., the doping profile or the heteroatom substitutions); typical control goals are the spatial electron-hole charge separation or the size of the HOMO-LUMO gap. The analysis in [10] relies on the simplyfing assumption of bounded domains. Here we generalize this analysis to unbounded domains. The main difficulty again consists in handling loss of mass at infinity.

Lemma 5** (Analytic properties of the set of HOMO-LUMO excitations).**

Consider a positively charged system, i.e. $Z>N$, then the joint solution set to the governing variational pronciples for occupied KS orbitals, HOMO and LUMO parametrized by the set of nuclear charge distributions $\mu$,

[TABLE]

has the following properties:

1. (a)

$B$* is weak $\times$ weak $\times$ weak $\times$ weak*∗-closed in $H^{1}(\mathbb{R}^{3})^{n}\times H^{1}(\mathbb{R}^{3})\times H^{1}(\mathbb{R}^{3})\times\mathcal{M}$ 2. (b)

$B$* is strong $\times$ strong $\times$ strong $\times$ weak*∗-compact in $H^{1}(\mathbb{R}^{3})^{n}\times H^{1}(\mathbb{R}^{3})\times H^{1}(\mathbb{R}^{3})\times\mathcal{M}$

Proof.

Part (a) Let $(\Phi^{(\nu)},\varphi_{H}^{(\nu)},\varphi_{L}^{(\nu)})\rightharpoonup(\Phi,\varphi_{H},\varphi_{L})$ in $H^{1}(\mathbb{R}^{3})^{(n+2)}$ and $\mu^{(\nu)}\rightharpoonup^{*}\mu$ in $\mathcal{M}$, then we need to prove

[TABLE]

Ad (i): Since all measures $\mu^{(\nu)}$ are supported on the compact set $\Omega_{nuc}$, the constant functions are in the predual of $\mathcal{M}$, that is the space of the continuous functions on $\Omega_{nuc}$, and hence $Z=\mu^{(\nu)}\left(\Omega_{nuc}\right)\to\mu\left(\Omega_{nuc}\right)$, so $\mu\in A_{nuc}.$

**Ad (ii): ** For any admissible $\Psi\in\mathcal{A}$ we have by the variational definition of $\Phi^{(\nu)}$ and the weak∗ continuity of $\mu\mapsto\mathcal{E}_{\mu}[\Psi]$

[TABLE]

which implies

[TABLE]

since $\Psi\in\mathcal{A}$ was arbitrary.

Using the weak $\times$ weak∗ lower semicontinuity of $(\Phi,\mu)\mapsto\mathcal{E}_{\mu}[\Phi]$ on $H^{1}(\mathbb{R}^{3})^{n}\times\mathcal{M}$ gives

[TABLE]

but unlike in the bounded domain case this does not directly give us the result since $\mathcal{A}$ is not weakly closed.

Fortunately, the $\Phi^{(\nu)}$ are not arbitrary elements of $\mathcal{A}$. Since the energy functional is invariant under unitary transformations, we can always assume the orbitals $\Phi=(\varphi_{1},\ldots,\varphi_{n})$ to be orthogonal. So the only thing that could go wrong is loss of mass due to weak convergence.

Assume for contradiction $\alpha=||\Phi||<n$ so we lose mass in at least one orbital $\varphi_{k}$. In order to see that this is not possible, we will place some small mass at a large but finite distance, and obtain a state with lower energy.

We take $\eta\in C^{\infty}_{c}(\mathbb{R}^{3})$ with $||\eta||_{2}=1$ and consider the test function $\eta_{\lambda,\sigma}:=\sigma^{\nicefrac{{1}}{{2}}}\lambda^{\nicefrac{{3}}{{2}}}\eta(\lambda\cdot)$.

Then

[TABLE]

Using the assumption (3.2) and the fact that $V[|\eta|^{2}]\leq 0$, we obtain that for $\sigma$ small enough there exists a constant $c>0$ such that

[TABLE]

Since $q<\tfrac{3}{2}$, the negative term dominates in the limit $\sigma,\lambda\to 0$ (if we let both go to [math] at the same speed). Hence by choosing the parameters $\lambda$ and $\sigma$ small enough we ensure $\mathcal{E}_{\mu}[\eta_{\lambda,\sigma}]<0$.

Now since we assume loss of mass in the k-th orbital $\varphi_{k}$, let us consider $\tilde{\varphi}^{(n)}_{k}(\cdot)=\varphi_{k}(\cdot)+\eta_{\lambda,\sigma}(\cdot-n\vec{e})$, where $\vec{e}$ is some unit vector in $\mathbb{R}^{3}$. Denoting the orbitals with $\varphi_{k}$ replaced by $\tilde{\varphi}_{k}$ by $\tilde{\Phi}$, we get

[TABLE]

So for $n$ large enough we obtain $\mathcal{E}_{\mu}[\tilde{\Phi}]<\mathcal{E}_{\mu}[\Phi]$ and $||\tilde{\Phi}||>||\Phi||$.

We can now repeat these steps until we have constructed a $\Psi$ with $\Psi\in\mathcal{A}$ (after a suitable unitary transformation), and arrive at the contradiction

[TABLE]

Hence there is no loss of mass, $\Phi\in\mathcal{A}$, and (ii) holds.

Before we move on to (iii) and (iv) we mention that due to the upper and lower bound above, $\mathcal{E}_{\mu^{(\nu)}}[\Phi^{(\nu)}]\to\mathcal{E}_{\mu}[\Phi]$. But this gives us $T[\Phi^{(\nu)}]\to T[\Phi]$ since the other terms in the energy functional, namely $V,J_{h},E_{xc}$, are continuous on $L^{2}\cap L^{4}$, and similarly to the proof of Lemma 4 we infer strong convergence in $L^{p}$ for $p\in[2,6)$, since there is no loss of mass.

Hence the kinetic energy converges as well, which means $||\nabla\Phi^{(\nu)}||_{2}\to||\nabla\Phi||_{2}$, giving us $\nabla\Phi^{(\nu)}\to\nabla\Phi$ in $L^{2}$, so $\Phi^{(\nu)}\to\Phi$ in $H^{1}(\mathbb{R}^{3})^{n}$.

Ad (iii): As in [10] the statements (iii) and (iv) are more difficult since the HOMO and LUMO orbitals $\varphi^{(\nu)}_{H}$ and $\varphi^{(\nu)}_{L}$ are not defined via universal but $\varphi^{(\nu)}$-dependent sets and hence an admissible trial function $\psi$ for the limiting HOMO and LUMO orbitals $\varphi_{H}$ and $\varphi_{L}$ may not be admissible for the variational principle for the approximating orbitals. In short, our argument in (4.1) is not valid anymore.

Hence we need to look at the $L^{2}$-projector of the sequence $\Phi^{(\nu)}=(\varphi_{1}^{(\nu)},\ldots,\varphi_{n}^{(\nu)})$, i.e. $\gamma_{\Phi^{(\nu)}}\chi:=\sum_{k=1}^{n}\langle\varphi_{k}^{(\nu)},\chi\rangle\varphi_{k}^{(\nu)}$. For any given $\chi\in L^{2}(\mathbb{R}^{3})$ the mapping $\Phi\mapsto\gamma_{\Phi}\chi$ is strongly continuous from $H^{1}(\mathbb{R}^{3})^{n}$ to $H^{1}(\mathbb{R}^{3})$, i.e.

[TABLE]

Furthermore by definition we have $\gamma_{\Phi}\chi=\chi$ and $||\chi||_{2}=1$ for any $\chi\in\mathcal{A}^{H}_{\Phi}$. So by (4.4) we have $||\gamma_{\Phi}^{(\nu)}\chi||_{2}>0$ for all $\nu$ large enough and therefore by the variational principle for the HOMO (3.3)

[TABLE]

Here we used the strong convergence of (4.4) and the fact that the map $(\Phi,\chi,\mu)\mapsto\mathcal{E}_{\mu,\rho}[\chi]$ is weak $\times$ strong $\times$ weak∗ continuous on $H^{1}(\mathbb{R}^{3})^{n}\times H^{1}(\mathbb{R}^{3})\times\mathcal{M}$ due to Lemma 2.

Since (4.5) holds for any $\chi\in\mathcal{A}^{H}_{\Phi}$, we have

[TABLE]

Next we prove that $\varphi_{H}\in A^{H}_{\Phi}$. Since by definition the $\varphi_{H}^{(\nu)}$ lie in the span of $\Phi^{(\nu)}$, we obtain by the weak convergence of $\varphi_{H}^{(\nu)}$ in $H^{1}(\mathbb{R}^{3})$ and the strong convergence of $\Phi^{(\nu)}$ in $H^{1}(\mathbb{R}^{3})^{n}$ that

[TABLE]

So taking the limit yields $\varphi_{H}=\gamma_{\Phi}\varphi_{H}$ and hence $\varphi_{H}^{(\nu)}\to\varphi_{H}$ strongly, so $\varphi_{H}\in\mathcal{A}^{H}_{\Phi}$.

Furthermore by the continuity properties proven in Lemma 2 we obtain

[TABLE]

Combining (4.6) and (4.7) yields (iii).

Ad (iv):

The corresponding proof for the LUMO starts similarly. Take any $\chi\in\mathcal{A}^{L}_{\Phi}$, i.e. $\gamma_{\Phi}\chi=0$ and $||\chi||_{2}=1$. Then again by (4.4) and the variational principle for the LUMO (3.3) we obtain

[TABLE]

Minimizing over $\chi\in A^{H}_{L}$ gives

[TABLE]

For the lower bound we use the weak $\times$ weak $\times$ weak∗ lower semicontinuity of the map $(\Phi,\mu,\chi)\mapsto\mathcal{E}_{\mu,\Phi}[\chi]$ by Lemma 2

[TABLE]

Unfortunately we are not done since the $\varphi_{L}^{(\nu)}$ are not known to converge strongly but just weakly, so we still need to prove that $\varphi_{L}$ is admissible. The limit $\varphi_{L}$ lies in the orthogonal complement of the $(\varphi_{k})_{k=1}^{n}$, since $\Phi^{(\nu)}$ converges strongly and therefore

[TABLE]

So we only need to prove $||\varphi_{L}||_{2}=1$. Assume $||\varphi_{L}||_{2}<1$. Then by (4.9) and (4.10)

[TABLE]

But we proved in Lemma 3 that the mass–to–LUMO–energy map $\lambda\mapsto I_{\lambda}$ is strictly decreasing. Hence $||\varphi_{L}||_{2}=1$, so $\varphi_{L}$ is admissible and together with (4.9) and (4.10) this establishes (iv).

Before we prove part b) let us make some remarks. Since there is no loss of mass, $\varphi_{L}^{(\nu)}$ converges strongly in $L^{2}(\mathbb{R}^{3})$ and hence strongly in $L^{p}(\mathbb{R}^{3})$ for $p\in[2,6)$. But since $\mathcal{E}_{\mu^{(\nu)},\Phi^{(\nu)}}[\varphi_{L}^{(\nu)}]\to\mathcal{E}_{\mu,\Phi}[\varphi_{L}]$ and all terms except the kinetic energy converge due to the continuity results in Lemma 2, we must also have $T[\varphi_{L}^{(\nu)}]\to T[\varphi_{L}]$, i.e. $\varphi_{L}^{(\nu)}$ converges strongly in $H^{1}(\mathbb{R}^{3})$.

In conclusion, we know that \big{(}\Phi^{(\nu)},\varphi_{H}^{(\nu)},\varphi_{L}^{(\nu)}\big{)} converges strongly in $H^{1}(\mathbb{R}^{3})^{(n+2)}$.

Part (b)

To prove sequential compactness of the set $\mathcal{B}$ is now quite easy since by part a) and Banach-Alaoglu we just need to prove that any sequence \big{(}\Phi^{(\nu)},\varphi_{H}^{(\nu)},\varphi_{L}^{(\nu)},\mu^{(\nu)}\big{)}\in\mathcal{B} is bounded in $H^{1}(\mathbb{R}^{3})^{(n+2)}\times\mathcal{M}$. For $\Phi^{(\nu)}$ and $\varphi_{L}^{(\nu)}$ this follows from the bounds in Lemma 1 and 2, respectively, noting that the exponent $p\in[1,\tfrac{5}{3})$ in the assumption on $v_{xc}$ has the property that the exponents $\tfrac{3(p-1)}{2}$ and $\tfrac{3(p-1)}{2p}$ of $||\rho||_{3}$ are strictly less than 1. Since the $\Phi^{(\nu)}$ stay bounded, so do the $\varphi_{H}^{(\nu)}$. Lastly, since $\mu^{(\nu)}\geq 0$ we have $||\mu^{(\nu)}||_{\mathcal{M}}=\int_{\mathbb{R}^{3}}\,\mathrm{d}\mu^{(\nu)}=Z$, which concludes the proof. ∎

As an example of a control goal we consider bandgap tuning as introduced in [10]. Here the quantity which one wants to influence by a suitable choice of the nuclear charge distribution $\mu$ is the HOMO-LUMO bandgap $\varepsilon_{H}-\varepsilon_{L}$, where $\varepsilon_{H}$ and $\varepsilon_{L}$ stand for the HOMO and LUMO eigenvalues of the KS Hamiltonian (2.11). Any bandgap tuning functional promoting a desired target value $\varepsilon_{\ast}$ has to reach its minimum when $\varepsilon_{H}-\varepsilon_{L}=\varepsilon_{\ast}$. A simple choice suggested in [10] is

[TABLE]

Theorem 3**.**

For any $\varepsilon_{\ast}>0$ and for $Z>N$, there exists a nuclear charge distribution $\mu\in\mathcal{A}_{nuc}$ which minimizes the bandgap tuning functional (4.12) over $\mathcal{A}_{nuc}$ subject to the constraints (3.3).

Proof.

The bandgap functional $J$ in (4.12) is, due to Lemma 2, weak $\times$ strong $\times$ strong $\times$ weak*∗* continuous on $(H^{1})^{n}\times H^{1}\times H^{1}\times\mathcal{M}$. Hence by the compactness of the set $\mathcal{B}$ proven in Lemma 5 it attains its minimum over this set. ∎

5 Nonexistence of HOMO-LUMO excitations in the neutral case $Z=N$

We now introduce carefully chosen and realistic model densities $\rho$ and prove that the excitation functional $\mathcal{E}_{\mu,\rho}$ admits no excited states, i.e. no bound states other than the ground state. See the picture on the right in Figure 1.

This finding suggests that also for the true KS ground state density $\rho$, it may happen that there are no exact HOMO-LUMO excitations. Of course, the Hamiltonian possesses continuous spectrum above the ground state energy and therefore “metastable” excitations (suitable square-integrable superpositions of continuous eigenstates) still exist.

From a mathematical point of view, the results in this section show that the assumption $Z>N$ in our existence result (Theorem 1) was in fact sharp and cannot be weakend to $Z\geq N$. Note that the model densities considered here satisfy the assumptions of Theorem 1, in particular (3.5).

H-Atom ground state density

Let $\mu=\delta_{0}$ and let $\varphi_{H}(x)=\tfrac{1}{\sqrt{\pi}}e^{-|x|}$ be the hydrogen atom ground state for the Schrödinger equation, i.e. the lowest eigenfunction of $-\tfrac{1}{2}\Delta-\tfrac{1}{|\cdot|}$. Its density is

[TABLE]

We expect this to be a good approximation for the KS density, hence the KS operator $h_{\rho_{H}}$ should be a good approximation to the self-consistent hydrogen KS operator.

In this case the Hartree-potential $v_{H}$ can be explicitly computed. The well known result is

[TABLE]

The first term cancels the external potential $v_{ext}$, hence the excitation functional becomes

[TABLE]

with corresponding Hamiltonian

[TABLE]

Model ground state density for the He-Atom

In order to construct a model density for Helium, we make the following ansatz.

We take a dilated version of the hydrogen orbital, i.e.

[TABLE]

and – following Hans Bethe – determine the parameter $\alpha$ by

[TABLE]

The last term describing the electron-electron interaction comes from

[TABLE]

with the spinors $\vec{\psi}_{1}=\begin{pmatrix}\varphi\\ 0\end{pmatrix}$ and $\vec{\psi}_{2}=\begin{pmatrix}0\\ \varphi\end{pmatrix}.$

The energy is easily computed as

[TABLE]

which implies

[TABLE]

In Figure 2 one sees that this value corresponds to the fact that the second electron does not see the full Coulomb potential of the nucleus but a screened one.

With this density the Hartree-potential can again be computed explicitly and plugging this density into our excitation functional gives us the Hamiltonian

[TABLE]

With these two densities at hand we can state the main result of this section and complete the picture given in Figure 1.

Theorem 4** (Spectrum of the KS Hamiltonian in the case $Z=N$).**

Consider either the hydrogen atom ( $N=1,\mu=\delta_{0}$) with the density $\rho=\rho_{H}$ given by (5.1) or the helium atom ($N=2,\mu=2\delta_{0}$) with the density $\rho=\rho_{\alpha}$ given by (5.2). Furthermore let the exchange-correlation energy be given by either Dirac exchange, PW81 or PZ92. Then the spectrum of the KS Hamiltonian has the form

[TABLE]

In particular, the Hamiltonian possesses exactly one bound state (up to spin in the hydrogen case) and no excited states, i.e. no bound states above the ground state.

The result of Theorem 4 is quite significant from a computational point of view. In numerical methods one, of course, obtains excited states, but in the limit of complete basis sets in infinite volume these might dissolve into metastable states associated with the continuous spectrum.

Proof.

As in the proof of Theorem 2 the potential of the KS Hamiltonian is in $L^{2}+L^{\infty}_{\varepsilon}$, hence in both cases \sigma_{ess}\big{(}h_{\mu,\rho}\big{)}=[0,\infty).

Next, we prove that there is at least one bound state with eigenvalue $\varepsilon_{0}<0$. By the Rayleigh–Ritz method it suffices to find a $\psi\in D(h_{\mu,\rho})=H^{1}(\mathbb{R}^{3})$ with $||\psi||_{2}=1$ and $\langle\psi,h_{\mu,\rho}\psi\rangle<0$. Since for any LDA functional we have $e_{xc}\leq e_{x}$, it suffices to prove the inequality for Dirac exchange. As a test function we choose the corresponding orbitals we used in the construction of our densities. These terms are easily computed to give for the hydrogen case

[TABLE]

So for both atoms we have at least one bound state.

Now we use the upper bound on the number of bound states given in [11] which says that the number $N_{\ell}$ of bound states with angular momentum $\ell$ satisfies

[TABLE]

where the functional $I_{p}(V)$ is given by

[TABLE]

For a radially symmetric potential this reduces to

[TABLE]

We call $I_{p}$ the Glaser-Martin-Grosse-Thirring (GMGT) functional. Here $V_{-}$ denotes the negative part of the potential, and the parameter $p\geq 1$ for the radially symmetric case while for the general case we have the restriction $p\geq\tfrac{3}{2}$.

If one calculates this integral numerically for our hydrogen density (5.1) with Dirac exchange, one obtains that the minimum value is attained at $p=1.4$ with $I_{p}(V)\approx 1.61587$. This gives us the upper bound

[TABLE]

which means $N_{0}\leq 1$ and all other $N_{\ell}=0$ for $\ell\geq 1$. So up to the degeneracy with respect to spin there is only the ground state and there are no excited states.

For the helium density (5.2) we computed the GMGT functional $I_{p}(V)$ for the Dirac [7], the Perdew-Zunger [21] and the Perdew-Wang [20] exchange-correlation functional and obtained

[TABLE]

Here we denote by $V_{-}$ the potential V_{-}(x)=2e^{-2\alpha|x|}\big{(}\frac{1}{|x|}+\frac{1}{\alpha}\big{)}+v_{xc}(\rho_{\alpha}(x)), so in e.g. the Dirac case we have

[TABLE]

As before, it now follows from (5.6) and (5.8) that no bound state other than the ground state orbital exists.

Hence for both hydrogen and helium we have exactly one bound state – the ground state itself – which corresponds to an eigenfunction with eigenvalue $\varepsilon_{0}<0$. ∎

In Figure 3 the values of $I_{p}(V)$ as a function of $p$ are given in the case of helium for the three LDAs mentioned above. The figure shows that the minimum value of the GMGT-functional is always attained at some $p\geq 1.5$, i.e. in the interval $p\in[\tfrac{3}{2},\infty)$, where the upper bound is also valid for non-symmetric potentials. Hence, as long as the real KS-density is only a small perturbation of our model density $\rho_{\alpha}$, even if it were not symmetric, our results would hold.

The overall conclusion of this section is that there are no LUMO excitation for the model densities (5.1) and (5.2).

Our findings raise the following interesting questions which lie beyond the scope of the present paper. First, is the GMGT nonexistence criterion satisfied for numerically obtained ground state densities of hydrogen and helium, or more complex atoms and molecules? A second question is whether the absence of exact excitations persists for more advanced excitation models like the Casida ansatz.

Acknowledgements. This project was supported by the DFG through the IGDK 1754 Optimization and Numerical Analysis for Partial Differential Equations with Nonsmooth Structures.

6 Appendix: Analytical properties of PZ81 and PW92

In the following we verify explicitly that the correlation functionals PZ81 and PW92 satisfy the assumptions (3.1) and (3.2) of this paper. In the physics literature the exchange-correlation functional is usually specified by the energy per particle at the density $\rho$, denoted $\varepsilon_{c}(\rho)$. We work with the mathematically convenient energy per unit volume $e_{xc}(\rho)=\rho\varepsilon_{c}(\rho)$. The exchange part was already discussed in the Remark Remark, so we only need to deal with the correlation part.

Furthermore, recall the Wigner-Seitz radius

[TABLE]

which is a standard parameter in physics to describe the local electron density of a system. Lastly, we remark that in the following $C$ will describe a generic constant, which may have different values at each appearance, but is independent of $\rho$ and $r_{s}$.

Perdew-Wang (PW92)

In this paper we consider only spin-unpolarized systems, so we have $\zeta=\frac{n_{\uparrow}-n_{\downarrow}}{n_{\uparrow}+n_{\downarrow}}=0$. So (in the notation of the original paper) we only need to check the assumptions (3.1) and (3.2) for $e_{c}(r_{s},0)$. The PW92 correlation functional is given by

[TABLE]

In order to improve readability of the arguments below, we define the following functions

[TABLE]

So returning to the notation of our paper we need to check (3.1) and (3.2) for the function

[TABLE]

This function is clearly continuously differentiable for $\rho>0$, so we only check the limit $\rho\to 0$. Since $f(r_{s}(\rho))=O(r_{s})=O(\rho^{-\nicefrac{{1}}{{3}}})$ and $g(r_{s})=\Omega(\sqrt{r_{s}})$ for $\rho\to 0$, repectively $r_{s}\to\infty$, we get

[TABLE]

hence with $e_{c}(0)=0$, $e_{c}$ is continuous.

Next we calculate the derivative

[TABLE]

Applying the inequality $\log(1+x)\leq x$ for $x>-1$, using $\tfrac{1}{g(r_{s})}=O(r_{s}^{-2})$ and taking the limit $\rho\to 0,r_{s}\to\infty$ gives now

[TABLE]

So also $v_{c}$ is continuous and hence it suffices to prove $v_{c}\leq c_{xc}(1+\rho^{p-1})$ for $\rho\to\infty$, i.e. $r_{s}\to 0$. Using $f(r_{s})\to C$,$g(r_{s})\to 0$ and $\tfrac{1}{g(r_{s})}=O(r_{s}^{-\nicefrac{{1}}{{2}}})$ we obtain

[TABLE]

so (3.1)holds with $p=\tfrac{7}{6}$ in.

For (3.2) we can choose $q=\tfrac{17}{12}$ (any value between $\tfrac{4}{3}$ and $\tfrac{3}{2}$ will do). Writing the condition in terms of $r_{s}$ then gives us

[TABLE]

Realizing that $f(r_{s})\sim-r_{s}$ for $r_{s}\to\infty$ and that $\tfrac{x}{1+x}\leq\log(1+x)\leq x$ implies for the $\log$-term $\Theta(\tfrac{1}{1+g(r_{s})})=\log(1+\tfrac{1}{g(r_{s})})$, transforms this into

[TABLE]

which is true if $q>\tfrac{4}{3}$.

Perdew-Zunger (PZ81)

Here we again consider spin-unpolarized systems, i.e. $\zeta=0$, as in the case of PW92. The precise value of the PZ81 constants is important for continuity and continuous differentiability; they were chosen in such a way that $e_{c}$ is continuously differentiable.

The PZ81 correlation functional is given by

[TABLE]

Here we used the notation of the original paper and hence $\gamma,B,D$ are negative. For our analysis we now need to consider $e_{c}=\rho\varepsilon(r_{s}(\rho))$. Hence, we calculate $v_{c}$ to be

[TABLE]

The continuity for $\rho\to 0$ is now checked easily:

[TABLE]

and

[TABLE]

The continuity of $e_{c}$ and $v_{c}$ at the value $r_{s}=1$ follows from the choice of constants in the original paper [21], since

[TABLE]

and

[TABLE]

Since $v_{c}$ is continuous, we only need to check the bound (3.1) for $\rho\to\infty$, i.e. $r_{s}\to 0$. But here the only term which is not bounded is $A\log(r_{s})$ for which we again use a standard $\log$-bound, $0\geq\log(x)\geq 1-\tfrac{1}{x}$ for $0<x\leq 1$, so |v_{c}(r_{s})|\leq C\big{(}1+\tfrac{1}{r_{s}}\big{)} for $r_{s}$ small enough. In terms of $\rho$ this means |v_{c}(\rho)|=C\big{(}1+\rho^{\nicefrac{{1}}{{3}}}\big{)}, so (3.1) holds with $p=\tfrac{4}{3}$.

Also, (3.2) holds with $q=\tfrac{4}{3}$, since then plugging in the relation (6.1) between $r_{s}$ and $\rho$ (6.1) yields

[TABLE]