Kohn-Sham theory with paramagnetic currents: compatibility and functional differentiability

TL;DR

This paper extends Moreau-Yosida regularization to paramagnetic current-density-functional theory, providing a rigorous mathematical framework and a Kohn-Sham iteration scheme for magnetic field effects, with numerical illustrations.

Contribution

It introduces a compatible convex formulation for paramagnetic current-density-functional theory using Moreau-Yosida regularization, including a well-defined Kohn-Sham scheme and numerical results.

Findings

Regularized functionals computed for a quantum ring.

Compatibility conditions for function spaces analyzed.

Numerical performance demonstrated for different regularizations.

Abstract

Recent work has established Moreau-Yosida regularization as a mathematical tool to achieve rigorous functional differentiability in density-functional theory. In this article, we extend this tool to paramagnetic current-density-functional theory, the most common density-functional framework for magnetic field effects. The extension includes a well-defined Kohn-Sham iteration scheme with a partial convergence result. To this end, we rely on a formulation of Moreau-Yosida regularization for reflexive and strictly convex function spaces. The optimal $L^p$-characterization of the paramagnetic current density $L^1\cap L^{3/2}$ is derived from the $N$-representability conditions. A crucial prerequisite for the convex formulation of paramagnetic current-density-functional theory, termed compatibility between function spaces for the particle density and the current density, is pointed out and analyzed. Several results about compatible function spaces are given, including their recursive construction. The regularized, exact functionals are calculated numerically for a Kohn-Sham iteration on a quantum ring, illustrating their performance for different regularization parameters.

Equations

$$T(\vec{A})$$

$$W$$

$$H^{\lambda }(v,\vec{A})=T(\vec{A})+\lambda W+\sum _{k=1}^{N}v({\vec{r}}_k),$$

$$\mathcal{K}(\psi )=\frac{1}{2}\sum _{k=1}^N{\int }_{{\mathbb{R}}^{3N}}|{\nabla }_{{\vec{r}}_k}\psi {|}^2\mathrm{d}{\vec{r}}_1\dots \mathrm{d}{\vec{r}}_N<+\infty .$$

$${\mathcal{W}}_N=\{\psi :\Vert \psi {\Vert }_2$$

$$\text{ ψ is antisymmetric}\}.$$

$${\mathcal{D}}_N=$$

$$\sum _n{p}_n=1,p_n\ge 0,\mathrm{Tr}(\gamma T(\vec{0}))<+\infty \}.$$

$$E^{\lambda }(v,\vec{A})$$

$$=\underset{\gamma \in {\mathcal{P}}_N}{\inf }\mathrm{Tr}(\gamma H^{\lambda }(v,\vec{A}))$$

$$=\underset{\psi \in {\mathcal{W}}_N}{\inf }⟨\psi ,H^{\lambda }(v,\vec{A})\psi ⟩.$$

$${\rho }_{\psi }({\vec{r}}_1)$$

$${\vec{j}}_{\psi }({\vec{r}}_1)$$

$$\frac{1}{2}{\int }_{{\mathbb{R}}^3}|\nabla {\rho }_{\psi }^{1/2}{|}^2\mathrm{d}\vec{r}\le \mathcal{K}(\psi ),$$

$${\mathcal{I}}_N=\left\{\rho \in L^1({\mathbb{R}}^3):\rho \ge 0,\Vert \rho {\Vert }_1=N,{\rho }^{1/2}\in {\mathcal{H}}^1({\mathbb{R}}^3)\right\}$$

$${\left({\int }_{{\mathbb{R}}^3}|f(\vec{r}){|}^6\mathrm{d}\vec{r}\right)}^{1/3}\le S{\int }_{{\mathbb{R}}^3}|\nabla f(\vec{r}){|}^2\mathrm{d}\vec{r},$$

$$\Vert {\rho }_{\psi }{\Vert }_3\le S{\int }_{{\mathbb{R}}^3}|\nabla {\rho }_{\psi }^{1/2}{|}^2\mathrm{d}\vec{r},S=\frac{4}{3(16{\pi }^2{)}^{1/3}}.$$

$$\Vert f{\Vert }_p^p\le \Vert f{\Vert }_r^{r(p-q)/(r-q)}\Vert f{\Vert }_q^{q(r-p)/(r-q)}$$

$${\tau }_{\psi }(\vec{r})=N{\int }_{{\mathbb{R}}^{3(N-1)}}|{\nabla }_{\vec{r}}\psi {|}^2\mathrm{d}{\vec{r}}_2\dots \mathrm{d}{\vec{r}}_N,$$

$${\tau }_{\psi }^k(\vec{r})=N{\int }_{{\mathbb{R}}^{3(N-1)}}|{\nabla }_{r^k}\psi {|}^2\mathrm{d}{\vec{r}}_2\dots \mathrm{d}{\vec{r}}_N$$

$$\Vert \psi {\Vert }_{{\mathcal{H}}^1}^2={\int }_{{\mathbb{R}}^3}\left({\rho }_{\psi }(\vec{r})/N+{\tau }_{\psi }(\vec{r})\right)\mathrm{d}\vec{r}=1+2\mathcal{K}(\psi ).$$

$$\left|N{\int }_{{\mathbb{R}}^{3(N-1)}}{\psi }^{\ast }{\nabla }_{\vec{r}}\psi \mathrm{d}{\vec{r}}_2\dots \mathrm{d}{\vec{r}}_N\right|$$

$${\int }_{{\mathbb{R}}^3}|{\vec{j}}_{\psi }(\vec{r})|\mathrm{d}\vec{r}\le \frac{1}{2}(2N\mathcal{K}(\psi ){)}^{1/2},$$

$$f_{\psi }^k(\vec{r})=N{\int }_{{\mathbb{R}}^{3(N-1)}}{\psi }^{\ast }{\partial }_{r^k}\psi \mathrm{d}{\vec{r}}_2\dots \mathrm{d}{\vec{r}}_N,$$

$$|f_{\psi }^k(\vec{r}){|}^p\le {\left({\rho }_{\psi }(\vec{r}){\tau }_{\psi }^k(\vec{r})\right)}^{p/2}$$

$${\left({\int }_{{\mathbb{R}}^3}|f_{\psi }^k(\vec{r}){|}^p\mathrm{d}\vec{r}\right)}^{1/p}$$

$$\le {\left[{\left({\int }_{{\mathbb{R}}^3}{\rho }_{\psi }(\vec{r}{)}^{pq/2}\mathrm{d}\vec{r}\right)}^{1/q}{\left({\int }_{{\mathbb{R}}^3}{\tau }_{\psi }^k(\vec{r})\mathrm{d}\vec{r}\right)}^{p/2}\right]}^{1/p}$$

$$=\Vert {\rho }_{\psi }{\Vert }_{pq/2}^{1/2}(2{\mathcal{K}}^k(\psi ){)}^{1/2}.$$

$$\Vert j_{\psi }^k{\Vert }_1\le (2N{\mathcal{K}}^k(\psi ){)}^{1/2},\Vert j_{\psi }^k{\Vert }_{3/2}\le S^{1/2}2{\mathcal{K}}^k(\psi ),$$

$$\Vert j_{\psi }^k{\Vert }_1\le {\int }_{{\mathbb{R}}^3}|f_{\psi }^k(\vec{r})|\mathrm{d}\vec{r}\le (2N{\mathcal{K}}^k(\psi ){)}^{1/2}.$$

$$\Vert j_{\psi }^k{\Vert }_{3/2}\le S^{1/2}2{\mathcal{K}}^k(\psi ),$$

$$\Vert j_{\psi }^k{\Vert }_p\le {\left[N^{3-2p}{S}^{3p-3}(2{\mathcal{K}}^k(\psi ){)}^{4p-3}\right]}^{\frac{1}{2p}}$$

$$\frac{1}{2}{\int }_{{\mathbb{R}}^3}|\nabla \sqrt{{\rho }_{\psi }}{|}^2\mathrm{d}\vec{r}+\frac{1}{2}{\int }_{{\mathbb{R}}^3}\frac{|{\vec{j}}_{\psi }{|}^2}{{\rho }_{\psi }}\mathrm{d}\vec{r}\le \mathcal{K}(\psi ).$$

$${\int }_{{\mathbb{R}}^3}$$

Full text

1 Introduction

The theoretical foundation of density-functional theory (DFT) was established in a seminal paper by Hohenberg and Kohn 1. There it was proven that two potentials that differ by more than a constant cannot share the same ground-state particle density $\rho$ (see Eq. (3) for definition). This fact is referred to as the Hohenberg–Kohn (HK) theorem. Using this result, the Schrödinger equation was replaced by a minimization problem involving a universal density functional (HK variational principle). The work by Lieb 2 provided an abstract reformulation of DFT that eliminates some technical difficulties with the HK formulation and constitutes a more tractable framework for rigorous mathematical analysis. Lieb’s formulation relies on Legendre–Fenchel transformations between the ground-state energy and a universal density functional, analogous to the use of Legendre transformations in thermodynamics and classical mechanics. The HK theorem becomes recast into a fact about subgradients of convex functionals that are mapped one-to-one by Legendre–Fenchel transformations 3, 4.

As far as practical purposes are concerned, DFT was first converted into a feasible algorithm for electronic structure calculations by Kohn and Sham 5. Here, both the unknown density of the full system and the effective Kohn–Sham (KS) potential for the non-interacting system are solved for in an iterative manner. Even though the important question of convergence of this procedure has been addressed in several works 6, 7, 8, 9, 10, it has only very recently been answered positively for finite-dimensional settings. 11

The motivation to include current densities and not just the particle density is to obtain a universal functional modelling the internal energy of magnetic systems. In terms of Lieb’s Legendre–Fenchel description, the current couples to the vector potential that now also enters the theory to account for the magnetic field. Recent work in current-density-functional theory (CDFT) has been devoted to the extension of the HK theorem, the HK variational principle, and the KS iteration scheme to include current densities 12, 13, 14, as well as to highlight the complexity of such a generalization 3, 15, 16. Other approaches are feasible as well, e.g., the magnetic-field density-functional theory (BDFT) of Grace and Harris 17, where a semi-universal functional is employed instead. There exists also a convexified formulation, in which BDFT and paramagnetic CDFT are related to each other by partial Legendre–Fenchel transformations 18, 19. Furthermore, the physically important case of linear vector potentials (uniform magnetic fields) has been theoretically studied in linear-vector-potential density-functional theory (LDFT) without the need to include current densities 18. Works beyond the current density generalization exist too, e.g., spin-current density-functional theory, reduced-density-matrix-functional theory, Maxwell–Schrödinger density-functional theory (MDFT), and quantum-electrodynamical density-functional theory (QEDFT) 20, 21, 22, 23, 24. For more generalized density-functional theories see Ref. 25 and the references therein, e.g., the kinetic-density-functional theory of Sim et al. 26, the work of Ayers 27 on $k$-density-functional theory, and Higuchi-Higuchi 28 who explored the use of different physical quantities as variables of the theory.

For mathematical reasons, CDFT is formulated in terms of the paramagnetic current density $\vec{j}$ (see Eq. (4) for definition) rather than with the gauge-invariant total current density. A theoretical foundation in the sense of a HK theorem for the total current density has not yet been proven and its existence remains an open question in the general case 3, 15. However, even if such a result could be shown, a HK variational principle does not exist for the total current density 16. Circumventing these problems may require the Maxwell–Schrödinger variational principle in place of the standard one 23. For the CDFT that makes use of the paramagnetic current density, it is well-known that there are counterexamples that rule out any analogue of the HK theorem 14, 3, 15. Nevertheless, since the particle density and the paramagnetic current density determine the non-degenerate ground state (see Ref. 29 for results in the degenerate case), a universal Levy–Lieb 30, 2 constrained-search functional can be set up, as done by Vignale and Rasolt 12. This functional can be extended to a Lieb functional that in this case also depends on the paramagnetic current density (for a first attempt see Ref. 31 with the choice of domain $(\rho,\vec{j})\in(L^{1}\cap L^{3})\times\vec{L}^{1}$).

Since the Lieb functional within standard DFT suffers from non-differentiability 32, a property that CDFT inherits, we address this particular problem and formulate a regularized theory in a Banach space setting. We here apply our recent work 9 that also extends the mathematical formalism of paramagnetic CDFT in Ref. 31. The need for differentiability—a fact that is usually overlooked in textbooks—is connected to the variational derivation and analysis of the Kohn–Sham scheme. This task, in the setting of paramagnetic CDFT, is the main aim of this work.

To set up a rigorous CDFT including the corresponding Kohn–Sham scheme, which is borrowed from our previous work 9 and here baptized Moreau–Yosida–Kohn–Sham optimal damping algorithm (MYKSODA), we introduce and discuss the condition of compatibility between function spaces for the scalar and vector potentials on the one hand and for the paramagnetic current and the total physical current densities on the other. This condition is necessary both for the convex formulation of CDFT and the subsequent Moreau–Yosida regularization procedure. Moreover, to maintain compatibility the regularization procedure requires a Banach space formulation and we make use of our results employing reflexive Banach spaces 9. In this respect the approach presented here differs from that in standard DFT where a Hilbert space formulation has been previously considered 4, which does not allow the necessary compatibility in CDFT. However, to apply the Banach space formulation outlined in Ref. 9, a suitable function space for the paramagnetic current density first needs to be identified. The choice $\vec{L}^{1}$ from Ref. 31 cannot be used for this purpose since it is not reflexive. It is therefore crucial to first prove that the paramagnetic current density is an element of $L^{p}$ for some $1<p<+\infty$. We prove in Corollary 2 that each component of $\vec{j}$ is an element of $L^{3/2}$ under the assumption of finite kinetic energy.

This article is structured as follows. After introducing the basic quantum-mechanical model for paramagnetic CDFT in Sec. 2.1, we define suitable function spaces for particle and current densities in Sec. 2.2. In such a setting, the usual constrained-search functionals of DFT are defined and the energy functional and generalized Lieb functional are subsequently set up in Sec. 2.3. These functionals serve as the primary objects for a further study of the theory in terms of convex analysis. Here a first problem arises within CDFT: the lack of concavity of the energy functional. As a further ingredient of a well-defined Kohn–Sham iteration scheme, finiteness of the energy functional (and its concave version) is proven in Sec. 2.5. Like the authors recently showed 9, the variational Kohn–Sham construction can only be rigorously set up for a regularized theory. The respective form of Moreau–Yosida regularization is introduced and applied to the setting at hand in Sec. 3.1. Finally, the stage is set for a discussion of the Kohn–Sham iteration scheme in Sec. 3.2 and its precise formulation as MYKSODA in Sec. 3.3. We note possible convergence issues in the particular setting of a two-particle singlet state in Sec. 3.4. We conclude in Sec. 4 with a numerical study of the MYKSODA. For readers that are less familiar with Lebesgue spaces and functional analysis in general we recommend Refs. 33 and 34, while for the tools borrowed from convex analysis we point to Ref. 35.

2 Paramagnetic CDFT

2.1 Ground-state model

In what follows, we consider the Hamiltonian of an $N$-electron system to be specified by an external scalar potential $v:\mathbb{R}^{3}\to\mathbb{R}$ and an external vector potential $\vec{A}:\mathbb{R}^{3}\to\mathbb{R}^{3}$. The components of $\vec{A}$ and other vectors are denoted $A^{k}$ and are not to be confused with the Euclidean norm squared, $|\vec{A}|^{2}=\vec{A}\cdot\vec{A}$. The physical kinetic energy operator and electron-electron repulsion are given by (in SI based atomic units),

[TABLE]

The full Hamiltonian then reads

[TABLE]

where a scale factor $\lambda\geq 0$ is included in front of the electron-electron repulsion term. This means $H^{0}$ is the Hamiltonian of a non-interacting system while the usual interacting system is retrieved with $H^{1}$. This extra parameter is motivated by its usefulness when addressing the KS theory and is standard in the literature.

We consider wavefunctions $\psi=\psi(\vec{x}_{1},\dots,\vec{x}_{N})$, where $\vec{x}_{k}=(\vec{r}_{k},s_{k})$ is the spatial and spin coordinate of the $k$-th particle. The wavefunctions are antisymmetric elements of the $N$-electron space $L^{2}((\mathbb{R}^{3}\times\{\uparrow,\downarrow\})^{N})$, i.e., the usual Hilbert space of quantum mechanics. We will be interested in ground-state CDFT, where several options are available to treat the spin degrees of freedom. Firstly, we could formulate a theory for the global ground state, obtained through minimization over all spin degrees of freedom, in which case a spin-Zeeman term could also be included in the Hamiltonian. Secondly, we could instead formulate a theory for the lowest singlet state ($S^{2}=0$) or the lowest state with some other prescribed value of the spin quantum number $S^{2}$. For simplicity, we formulate a ground-state CDFT for the lowest singlet energy and adapt our notation accordingly. However, our analysis is mostly independent of this choice and applies equally well to a theory for global ground states. Thus without loss of generality the spin coordinate will in the sequel be omitted.

All wavefunctions are assumed to have finite kinetic energy,

[TABLE]

We further assume $L^{2}$ normalization of $\psi$ and denote the $L^{p}$ norm by $\|\cdot\|_{p}$, $1\leq p\leq\infty$. Henceforth the particle number $N$ will be fixed and we define the set of admissible wavefunctions

[TABLE]

Moreover, $\gamma_{\psi}=|\psi\rangle\langle\psi|$ denotes the density matrix of a pure state $\psi$ and $\mathcal{P}_{N}=\{\gamma_{\psi}:\psi\in\mathcal{W}_{N}\}$ is the set of such states. The set of mixed states is given by (where the sum over $n$ can be infinite)

[TABLE]

Note that in case $\gamma=\gamma_{\psi}\in\mathcal{P}_{N}$ we have $\operatorname{Tr}(\gamma T(\vec{0}))=\mathcal{K}(\psi)$.

The energy functional for the ground-state energy can be written in the following alternative forms

[TABLE]

Thus if a minimizer $\gamma\in\mathcal{D}_{N}$ exists one can always also obtain a pure ground state selected from one of the eigenvectors of $\gamma$.

2.2 Function spaces for densities

For any $\psi\in\mathcal{W}_{N}$, we define the particle density and the paramagnetic current density, respectively, according to

[TABLE]

The aim of this section is to extract as much information as possible about the regularity of $\rho_{\psi}$ and $\vec{j}_{\psi}$ in terms of $L^{p}$ spaces from the assumption that $\psi\in\mathcal{W}_{N}$. This will define the sets of admissible densities.

To avoid confusion, a word or two on our notation is appropriate at this point. Since the paramagnetic current density is the main current density of consideration we omit the usual superscript (or subscript) “p” for paramagnetic in $\vec{j}^{\mathrm{p}}$. We write $\underline{\vec{r}}=(\vec{r}_{1},\dots,\vec{r}_{N})\in\mathbb{R}^{3N}$, $\nabla=\nabla_{\underline{\vec{r}}}$, and let $\vec{r}$ denote any $\vec{r}_{i}\in\mathbb{R}^{3}$ but typically $\vec{r}_{1}$. Further $\mathcal{H}^{k}({\mathbb{R}}^{n})$ denotes the Sobolev space that includes all functions in $L^{2}({\mathbb{R}}^{n})$ with weak derivatives up to $k$-th order in $L^{2}({\mathbb{R}}^{n})$. ($\mathcal{H}$ should not be confused with the Hamiltonian $H$.) Finally, $\vec{X}=X\times X\times X$ is the triple copy of a Banach space $X$, here mostly used for $L^{p}$ spaces as $\vec{L}^{p}$.

Hoffmann-Ostenhof and Hoffmann-Ostenhof (see Eq. (3.10) in Ref. 36) and Lieb (Theorem 1.1 in Ref. 2) have shown that the von Weizsäcker term involving $\rho_{\psi}$ is bounded by the kinetic energy of $\psi$, i.e.,

[TABLE]

and therefore $\rho_{\psi}\in\mathcal{I}_{N}$ if $\psi\in\mathcal{W}_{N}$, where

[TABLE]

denotes the set of $N$-representable particle densities. Even though the Hilbert space $\mathcal{H}^{2}({\mathbb{R}}^{3N})$ is the natural domain of the kinetic energy operator, $\psi\in\mathcal{H}^{1}({\mathbb{R}}^{3N})$ is sufficient to guarantee finite $\mathcal{K}(\psi)$. The Sobolev inequality in $\mathbb{R}^{3}$ (see, e.g., Theorem 8.3 in Ref. 33),

[TABLE]

applied to $f=\rho_{\psi}^{1/2}$ yields

[TABLE]

Consequently $\mathcal{I}_{N}\subset L^{1}\cap L^{3}$.

Remark 1*.*

We make a brief comment concerning interpolation. For $1\leq q\leq p\leq r\leq+\infty$, if $f\in L^{q}\cap L^{r}$ then by Hölder’s inequality

[TABLE]

and thus $f\in L^{p}$.

We proceed by summarizing criteria that follow from the works of Lieb and Kato for the space of particle densities in terms of $L^{p}$ spaces.

Proposition 1** (Lieb 2 and Kato 37).**

For $\psi\in\mathcal{W}_{N}$, $\rho_{\psi}\in L^{p}$ with $p\in[1,3]$ and in particular $\rho_{\psi}$ is an element of the Hilbert space $L^{2}$. Moreover, $\psi\in\mathcal{H}^{2}$ implies $\rho_{\psi}\in L^{p}$ for all $p\in[1,\infty]$.

Proof.

The first part follows from Lieb 2 and Remark 1. By Lemma 3 in Ref. 37, the assumption $\psi\in\mathcal{H}^{2}$ gives $\rho_{\psi}\in L^{\infty}$. (Note that $A_{i}f(\vec{r}_{i})$ in Ref. 37 corresponds to $(\rho(\vec{r}_{i})/N)^{1/2}$ in our notation.) ∎

An obvious limitation for the current density is that every component of $\vec{j}_{\psi}$ is in $L^{1}$ 31. Here, by better exploiting the properties of $\psi\in\mathcal{W}_{N}$, we will be able to further characterize the set of current densities.

We start by giving some definitions. The kinetic-energy density $\tau_{\psi}:\mathbb{R}^{3}\rightarrow\mathbb{R}_{+}$ is given by

[TABLE]

and relates to the already defined kinetic energy by $\mathcal{K}(\psi)=\frac{1}{2}\|\nabla\psi\|_{2}^{2}=\frac{1}{2}\|\tau_{\psi}\|_{1}$. We see that $\psi$ being an element of $\mathcal{W}_{N}$ guarantees finite $\mathcal{H}^{1}$ norm and thus $\mathcal{W}_{N}\subset\mathcal{H}^{1}$. Furthermore, let $\vec{r}\in\mathbb{R}^{3}$ be written $\vec{r}=(r^{1},r^{2},r^{3})$ and define the component-wise kinetic energy density

[TABLE]

and $\mathcal{K}^{k}(\psi)=\frac{1}{2}\|\tau_{\psi}^{k}\|_{1}$. A direct computation gives that the usual Sobolev norm satisfies

[TABLE]

An important bound that will be used subsequently is the following one

[TABLE]

which is a consequence of the Cauchy–Schwarz inequality. Since both $\tau_{\psi}^{1/2}$ and $\rho_{\psi}^{1/2}$ are $L^{2}$ functions, integration and the Cauchy–Schwarz inequality give $\vec{j}_{\psi}\in\vec{L}^{1}$. Indeed

[TABLE]

see Proposition 3 in Ref. 31. The idea is now to further extend such $\vec{L}^{p}$-characterizations. With $\vec{j}_{\psi}=(j_{\psi}^{1},j_{\psi}^{2},j_{\psi}^{3})$ we now state and prove

Lemma 2**.**

Set $f_{\psi}^{k}=\partial_{k}\rho_{\psi}+2\mathrm{i}j_{\psi}^{k}$ for $k=1,2,3$ and $\psi\in\mathcal{W}_{N}$. Then $f_{\psi}^{k}\in L^{p}(\mathbb{R}^{3})$ for $1\leq p\leq 3/2$.

Proof.

Since

[TABLE]

we have similar to Eq. (8) for $1\leq p\leq 3/2$

[TABLE]

by the Cauchy–Schwarz inequality. Next, we use Hölder’s inequality with $q$ defined by $p/2+1/q=1$ such that

[TABLE]

To conclude, we note that, by the assumption on $p$, we have $1\leq pq/2\leq 3$ and recall that $\rho$ is in $L^{1}\cap L^{3}$. ∎

Note that $\vec{f}_{\psi}$ in Lemma 2 can be seen as a generalized complex current, similar to the one considered by Tokatly 38 in a lattice version of time-dependent CDFT. It has also been considered before as “momentum density” 39.

As a direct consequence of Lemma 2, we have our main result about function spaces for the current density

Theorem 3**.**

For $\psi\in\mathcal{W}_{N}$, each component $j_{\psi}^{k}$ of the paramagnetic current density $\vec{j}_{\psi}$ is in $L^{p}$ for any $1\leq p\leq 3/2$ and we write $\vec{j}_{\psi}\in\vec{L}^{p}$. In particular, we have

[TABLE]

with the constant $S$ given by Eq. (7).

Proof.

Set $p=1$ and apply Lemma 2. Since $pq/2=1$ and $\|\rho_{\psi}\|_{1}=N$, Eq. (10) gives

[TABLE]

With the choice $p=3/2$ instead, we have $pq/2=3$. Eq. (10) then reduces to

[TABLE]

where we have also used the Sobolev inequality (6) and Eq. (5). Using interpolation (see Remark 1), it follows

[TABLE]

and thus $j_{\psi}^{k}\in L^{p}$ for $1\leq p\leq 3/2$. ∎

From the proof of Lemma 2 (the Cauchy–Schwarz inequality applied to Eq. (9) with $p=2$), we obtain the current correction to the von Weizsäcker kinetic energy (see Ref. 40 and note that this sharpens Theorem 14 in Ref. 31)

[TABLE]

The well-known inequality $\frac{1}{2}\int_{{\mathbb{R}}^{3}}|\vec{j}_{\psi}|^{2}\rho_{\psi}^{-1}\,\mathrm{d}\vec{r}\leq\mathcal{K}(\psi)$ is immediate from Eq. (11).

We next note that the space for the current density cannot be further restricted since $j_{\psi}^{k}\notin L^{p}$, $p>3/2$, for some $\psi\in\mathcal{W}_{N}$. Before proving this fact, a further characterization of $\vec{j}_{\psi}$ using $\psi\in\mathcal{H}^{2}$ is given.

Proposition 4**.**

$\psi\in\mathcal{H}^{2}$* implies $\vec{j}_{\psi}\in\vec{L}^{2}$.*

Proof.

Suppose $\psi\in\mathcal{H}^{2}$. By Proposition 1 it follows that $\rho_{\psi}\in L^{\infty}$. Then from Eq. (8)

[TABLE]

which gives $\vec{j}_{\psi}\in\vec{L}^{2}$. (Note that we could have argued by means of the Cauchy–Schwarz inequality to obtain $\int_{\mathbb{R}^{3}}|\vec{j}|^{2}\,\mathrm{d}\vec{r}\leq\|\rho\|_{\infty}\int_{\mathbb{R}^{3}}|\vec{j}|^{2}\rho^{-1}\,\mathrm{d}\vec{r}$, where the integral over $|\vec{j}|^{2}\rho^{-1}$ is bounded in terms of $\mathcal{K}(\psi)$.) ∎

Proposition 5**.**

For $N=1$, there exists $\phi\in\mathcal{W}_{N}$ such that $\vec{j}_{\phi}\notin\vec{L}^{p}$, for every $p>3/2$.

Proof.

Consider the function

[TABLE]

where $0\leq r\leq 1$, $\alpha>-1$, and $2\beta>-\alpha-1$. It then holds that $\phi\in\mathcal{H}^{1}(\mathbb{R}^{3})$, since

[TABLE]

We wish to show that $\vec{j}_{\phi}\notin\vec{L}^{p}$, for all $p>3/2$ for some choice of $(\alpha,\beta)$ in the set

[TABLE]

Note that $\vec{j}_{\phi}=\text{Im}\,\phi^{*}\nabla\phi=r^{\alpha}\nabla r^{\beta}$. Let $\delta$ be an arbitrarily small positive number and set $\alpha=-1+3\delta$. Then with $\beta=-\delta>-3\delta/2$ we have $(\alpha,\beta)\in\Theta$. For all $p>3/(2(1-\delta))$ we obtain

[TABLE]

Since any $\delta>0$ is allowed, we conclude $\vec{j}_{\phi}\notin\vec{L}^{p}$, $p>3/2$. ∎

Remark 2*.*

Note that the same counterexample also shows that $\rho_{\phi}\notin L^{p}$, for $p>3$.

From $\psi\in\mathcal{W}_{N}$, with Proposition 1 and Remark 2 we have thus arrived at the well-known optimal choice $L^{1}\cap L^{3}$ of $L^{p}$ spaces for the particle density 2. Similarly, by Theorem 3 and Proposition 5 the optimal choice for the paramagnetic current density is $\vec{L}^{1}\cap\vec{L}^{3/2}$. Note that densities and currents that are not from these spaces cannot be represented by admissible wavefunctions $\psi\in\mathcal{W}_{N}$. Later this choice will be further limited by the demands coming from compatibility (see Sec. 2.4), reflexivity, and strict convexity (in connection with the regularized KS iteration scheme, Sec. 3.3).

We summarize this section with some definitions and a concluding corollary. We also refer to Refs. 41, 2, 42, 39, 43, 44 for further discussions on this topic (not only confined to CDFT).

Definition 6**.**

We say that a density pair $(\rho,\vec{j})$ is $N$-representable if there exists a $\psi\in\mathcal{W}_{N}$ such that $\rho_{\psi}=\rho$ and $\vec{j}_{\psi}=\vec{j}$. If such a $\psi$ is the ground state of some $H^{\lambda}(v,\vec{A})$, then $(\rho,\vec{j})$ is ($\lambda$) $v$-representable. Furthermore, we distinguish between fully interacting ($\lambda=1$) and non-interacting ($\lambda=0$) $v$-representability. If $\psi\in\mathcal{W}_{N}$ is replaced by $\gamma\in\mathcal{D}_{N}$ we call the above property ensemble $v$-representability.

Corollary 7**.**

The set of $N$-representable density pairs $(\rho,\vec{j})$ is a subset of $(L^{1}\cap L^{3})\times(\vec{L}^{1}\cap\vec{L}^{3/2})$.

The above set is given in terms of $L^{p}$ spaces only. There are other well established constraints such as $\rho(\vec{r})\geq 0$, $\int_{{\mathbb{R}}^{3}}\rho\,\mathrm{d}\vec{r}=N$, and $\int_{{\mathbb{R}}^{3}}|\vec{j}|^{2}\rho^{-1}\,\mathrm{d}\vec{r}<+\infty$ 2, 42, 31.

It will later be important to impose a reflexive (R) and strictly convex Banach space setting, and we therefore define such a space that includes the set of $N$-representable density pairs. Hanner’s inequality 45 shows that $L^{p}$, $1<p<\infty$, is even a uniformly convex space which implies both strict convexity and reflexivity.

Definition 8**.**

$X_{\mathrm{R}}\times Y_{\mathrm{R}}=L^{3}\times\vec{L}^{3/2}$.

2.3 Constrained-search functionals

To formulate a rigorous CDFT, several requirements need to be placed on densities and potentials. Some of these requirements are related to $N$-representability and thus do not amount to any restriction, but merely exclude irrelevant densities that are invalid in the sense that they cannot arise from any quantum-mechanical state $\psi\in\mathcal{W}_{N}$. Other requirements amount to assumptions about the ground-state densities or restrictions on the external potentials that can be considered.

The universal part of the Hamiltonian $H^{\lambda}(v,\vec{A})$ in Eq. (1), independent of the potential pair $(v,\vec{A})$, is $H^{\lambda}(0,\vec{0})=T(\vec{0})+\lambda W$. On a Banach space $X\times Y\subset X_{\mathrm{R}}\times Y_{\mathrm{R}}$ of measurable functions for particle densities ($X$) and current densities ($Y$), define the universal Levy–Lieb-type functionals $F_{\mathrm{VR}}^{\lambda}$ and $F_{\mathrm{VR,DM}}^{\lambda}$ according to

[TABLE]

These functionals are derived from the parts of the energy expressions in Eq. (2) that are independent of the potential pair $(v,\vec{A})$. In analogy with Eq. (2) we generally have two possibilities, searching either over pure or mixed states. If a given density pair $(\rho,\vec{j})$ cannot be represented by a pure or mixed state, then the value of the functional will just be set to $+\infty$ by definition. Unlike the ground-state energy functional in Eq. (2), the pure and mixed search domains do not yield equivalent results, since the former is subject to more severe representability restrictions 39.

Remark 3*.*

The density functionals are here denoted “VR” which stands for Vignale and Rasolt to credit their work 12. We remark that we could just as well have chosen to credit Levy, Valone, and Lieb due to the obvious counterpart in DFT 30, 2, 46. “DM” refers to the use of density matrices for mixed states and was introduced for the DFT constrained-search functional by Valone in Ref. 46.

Remark 4*.*

At this point we wish to keep the setting general and thus the space $X\times Y$ is not specified any closer. This setting includes the choice $X=L^{1}\cap L^{3}$, $Y=L^{1}\cap L^{3/2}$ such that all information from the previous section on function spaces is used. However, if a regularized theory is to be obtained, we need to have a reflexive Banach space setting and therefore we also consider the less restrictive choice $X=X_{\mathrm{R}}$, $Y=Y_{\mathrm{R}}$.

The Levy–Lieb-type functional $F_{\mathrm{VR}}(\rho,\vec{j})$ is not convex, see Proposition 8 in Ref. 31. Yet by the linearity of the map $\gamma\mapsto(\rho,\vec{j})$ it follows that $F_{\mathrm{VR,DM}}(\rho,\vec{j})$ is convex. Both functionals are admissible 47 in the sense that they can be used to compute the ground-state energy.

Since the energy expression will naturally include integrals over couplings of potentials with densities, it is helpful to introduce the notion of dual pairings (between elements of dual Banach spaces). For measurable functions $f,g$ with domain ${\mathbb{R}}^{3}$ let

[TABLE]

whenever the integral is well-defined in $\mathbb{R}\cup\{\pm\infty\}$, and similarly for vector-valued functions $\vec{f},\vec{g}$, but with the pointwise product replaced by $\vec{f}\cdot\vec{g}$. Then the energy expressions in Eq. (2) can be written as

[TABLE]

and equivalently by employing $F^{\lambda}_{\mathrm{VR,DM}}$ defined with mixed states instead of $F_{\mathrm{VR}}^{\lambda}$. In particular, $\lambda=1$ corresponds to the fully interacting system, and $\lambda=0$ to a non-interacting one.

At the outset, the formulation of paramagnetic CDFT relies on a decomposition of the total kinetic energy into canonical kinetic energy, the paramagnetic term $\langle\vec{A},\vec{j}\rangle$, and the diamagnetic term $\langle|\vec{A}|^{2},\rho\rangle$, with each of the terms separately finite 3, 15, 48. As in standard DFT, the electrostatic interaction with the external potential, $\langle v,\rho\rangle$, needs to be finite too. In the convexified form, a new potential variable is formed by absorbing the diamagnetic term into the scalar potential $u=v+\frac{1}{2}|\vec{A}|^{2}$ 3, 31. Minimally, then, the underlying function spaces should be such that

[TABLE]

A convex formulation achieves this automatically as it requires the stronger condition that densities and potentials are elements of dual Banach spaces:

Definition 9** (Density-potential duality).**

We say that there is density-potential duality, or just duality, when densities and potentials are confined to dual Banach spaces

[TABLE]

Remark 5*.*

At this moment we do not assume any more specific properties for $X\subset X_{\mathrm{R}}$ and $Y\subset Y_{\mathrm{R}}$ besides duality. However, reflexivity and strict convexity of $X$ and $Y$ are additional assumptions that will become important in Sec. 3.

Remark 6*.*

Note that $u\in X^{*}$, $X=L^{1}\cap L^{3}$ and $\vec{A}\in Y^{*}$, $Y=L^{1}\cap L^{3/2}$ imply $u\in L^{3/2}+L^{\infty}$ and $\vec{A}\in\vec{L}^{3}+\vec{L}^{\infty}$. For the condition on the original scalar potential $v$ see the next section on compatibility. As far as the vector potential is concerned, the restrictions on $\vec{A}$ are stronger than the familiar setting of $\vec{A}\in L_{\mathrm{loc}}^{2}$ (see e.g., Ref. 33), which is implied by $\vec{A}\in\vec{L}^{3}+L^{\infty}$. Also the reflexive setting with $X=X_{\mathrm{R}}$ and $Y=Y_{\mathrm{R}}$, where $u\in L^{3/2}$ and $\vec{A}\in\vec{L}^{3}$, implies $\vec{A}\in L_{\mathrm{loc}}^{2}$ again. Moreover, $u\in L_{\mathrm{loc}}^{3/2}$ is a natural assumption 33. We remark that our consideration of dual spaces is mathematically motivated and not a physical necessity. A truncated space domain is in many cases needed to cover the usual potentials of physical systems (see Ref. 4 for a discussion on this topic).

Since the potentials $(v,\vec{A})$ are not paired linearly with the densities $(\rho,\vec{j})$, the functional $E^{\lambda}$ defined in this way is not concave. The change of variables $u=v+\tfrac{1}{2}|\vec{A}|^{2}$ results in a convexification of paramagnetic CDFT, meaning that

[TABLE]

is a jointly concave functional 3. The consequences of this variable change for the choice of function spaces will be discussed in Sec. 2.4. The price to pay for concavity is a convoluted gauge symmetry. For all scalar fields $\chi$ with gradients in the same function space as $\vec{A}$ one has

[TABLE]

But the benefit is much greater, making $\bar{E}^{\lambda}$ jointly concave in both potentials and highlighting the linear structure of coupling between potentials and densities

[TABLE]

The convex formulation of paramagnetic CDFT can be outlined as follows. Let the dual space of $X\times Y$ be given by $X^{*}\times Y^{*}$. We define the generalized Lieb functional $F^{\lambda}(\rho,\vec{j})$ as the supremum over the energy plus the linear coupling between densities and potentials, i.e.,

[TABLE]

Such a functional is by construction convex. To extract more from Eq. (13) we first need

Definition 10**.**

Let $B$ be a Banach space with dual $B^{*}$, $f:B\to{\mathbb{R}}\cup\{\pm\infty\}$, and $g:B^{*}\to{\mathbb{R}}\cup\{\pm\infty\}$.

• (i)

If $f$ is convex, lower semi-continuous, has $f>-\infty$, and is not identically equal to $+\infty$, then it is called closed convex and we write $f\in\Gamma_{0}(B)$. Analogously, with weak-* lower semi-continuity we define $\Gamma_{0}^{*}(B^{*})$. We also introduce the sets $\Gamma(B)=\Gamma_{0}(B)\cup\{\pm\infty\}$ and $\Gamma^{*}(B^{*})=\Gamma^{*}_{0}(B^{*})\cup\{\pm\infty\}$.

• (ii)

Following Refs. 2, 4, we define the (skew) conjugate functionals (Legendre–Fenchel transformations),

[TABLE]

Definition 11**.**

The standard norms for the intersection of two Banach spaces $B,B^{\prime}$ and their set sum are given by the following expressions 49

[TABLE]

Theorem 3.6 in Lieb 2 can be straightforwardly generalized to the statement that $F^{\lambda}$ is lower semi-continuous on the space $X\times Y=(L^{1}\cap L^{3})\times(\vec{L}^{1}\cap\vec{L}^{3/2})$, with topology defined by the norm from the above definition. (See also Proposition 12 in Ref. 31 where this was done for $(L^{1}\cap L^{3})\times\vec{L}^{1}$.) Also, by the same argument, we have for the reflexive setting

Lemma 12**.**

$F^{\lambda}\in\Gamma_{0}(X_{\mathrm{R}}\times Y_{\mathrm{R}})$.

Since $F^{\lambda}$ is convex and lower semi-continuous, this clears the way for the application of powerful tools from convex analysis. Because Eq. (13) is already the Legendre–Fenchel transformation of the energy functional $\bar{E}^{\lambda}$ we are able to switch back from $F^{\lambda}$ to $\bar{E}^{\lambda}$ with the inverse transformation (see Theorem 1 in Ref. 9)

[TABLE]

Using the notation from Definition 10 (ii), we sum up the situation as

[TABLE]

Moreover the closed convex $F^{\lambda}$ is the smallest possible admissible functional,

[TABLE]

Solving the variational problem in Eq. (14) is the general task of CDFT.

Remark 7*.*

It is to the best of our knowledge an open question whether $F_{\mathrm{VR,DM}}^{\lambda}$ is lower semi-continuous and hence equal to $F^{\lambda}$ in the context of paramagnetic CDFT.

2.4 Compatibility of function spaces

Duality as in Definition 9 is not strong enough to guarantee finiteness of the diamagnetic term. It also does not guarantee another natural condition on the function space for the diamagnetic contribution to the current density that we call compatibility.

Definition 13**.**

The density function space $X$ and the current density function space $Y$ are said to be compatible if, for all $\rho\in X$ and all $\vec{A}\in Y^{*}$,

[TABLE]

We emphasize both conditions in the definition above, as they have different physical interpretations, although it will be seen in Theorem 14 that (C1) and (C2) are equivalent. The first condition (C1) requires that the scalar potential $v$ and $|\vec{A}|^{2}$ share the same function space, so that changes of variables between $v=u-\frac{1}{2}|\vec{A}|^{2}$ and $u=v+\frac{1}{2}|\vec{A}|^{2}$ stay within the space $X^{*}$. The second condition (C2) requires that the paramagnetic contribution, $\vec{j}$, and the diamagnetic contribution, $\rho\vec{A}$, to the total current density share the same function space $Y$.

Compatibility also ensures a sensible behavior under gauge transformations. Duality imposes a restriction on the gauge transformations that are allowed within the theory. Any gauge function $\chi$ that is used to transform $\vec{A}$ to $\vec{A}^{\prime}=\vec{A}+\nabla\chi$ must satisfy $\nabla\chi\in Y^{*}$. If $(v,\vec{A})$ has the ground-state density pair $(\rho,\vec{j})$, the gauge transformed potential pair $(v,\vec{A}^{\prime})$ would be expected to have the ground-state density pair $(\rho,\vec{j}^{\prime}=\vec{j}+\rho\nabla\chi)$. The second compatibility condition (C2) ensures that $\vec{j}^{\prime}\in Y$, so that ground-state density pairs are never lost after an allowed gauge transformation.

The following theorem shows that the two compatibility conditions are in fact equivalent. First, we mention a fundamental result that will be used in the sequel. Suppose we are given a general Banach space $B$ of measurable functions and a measurable function $g$. We can then check that $g$ is contained in the dual $B^{*}$ by verifying that the pairing $\langle g,f\rangle=\int\!g\,f\,\mathrm{d}\vec{r}$ is finite for all $f\in B$, see Appendix A for a full proof of this statement.

Theorem 14**.**

If $X,Y$ and their duals are Banach spaces of measurable functions, then (C1) $\iff$ (C2).

Proof.

Part 1 ($\impliedby$): From (C2), we have that

[TABLE]

for all $\rho\in X$ and all $\vec{A},\vec{A}^{\prime}\in Y^{*}$. Specialization to the case $\vec{A}^{\prime}=\vec{A}$ yields (C1).

Part 2 ($\implies$): Suppose (C2) is false, i.e., $\rho\vec{A}\notin Y$. Then there exists an $\vec{A}^{\prime}\in Y^{*}$ such that

[TABLE]

and we obtain

[TABLE]

Thus $|\vec{A}|^{2}+|\vec{A}^{\prime}|^{2}\notin X^{*}$ and either $|\vec{A}|^{2}$ or $|\vec{A}^{\prime}|^{2}$ is not an element of $X^{*}$ (or both). This contradicts (C1). ∎

Compatible function spaces can be built up recursively, by combining different function spaces that are already known to be compatible.

Theorem 15** (Compatibility of intersections and sums).**

Suppose $X_{1}$ and $Y_{1}$ are compatible, and the same holds for $X_{2}$ and $Y_{2}$. Then (a) the intersections

[TABLE]

are compatible. Moreover, (b) the sums

[TABLE]

are compatible. (Norms of intersections and sums are as given in Definition 11.)

Proof.

By Theorem 14 it is sufficient to prove (C1).

Part (a): The dual spaces are $X^{*}=X_{1}^{*}+X_{2}^{*}$ and $Y^{*}=Y_{1}^{*}+Y_{2}^{*}$. Decompose the vector potential as $\vec{A}=\vec{A}_{1}+\vec{A}_{2}$, with $\vec{A}_{1}\in Y_{1}^{*}$ and $\vec{A}_{2}\in Y_{2}^{*}$. Using the inequality

[TABLE]

property (C1) follows immediately, since

[TABLE]

and each of the two terms is finite by hypothesis.

Part (b): The dual spaces in this case are $X^{*}=X_{1}^{*}\cap X_{2}^{*}$ and $Y^{*}=Y_{1}^{*}\cap Y_{2}^{*}$. Trivially, for any $\vec{A}\in Y^{*}$ we have $\vec{A}\in Y_{i}^{*}$ and therefore, by hypothesis, $|\vec{A}|^{2}\in Y_{i}^{*}$ ($i=1,2$). Hence, $|\vec{A}|^{2}\in Y_{1}^{*}\cap Y_{2}^{*}$. ∎

As demonstrated in Ref. 31 (see also Theorem 3 above), the paramagnetic current density satisfies $\vec{j}_{\psi}\in\vec{L}^{1}$ for $\psi\in\mathcal{W}_{N}$. Combined with compatibility, this becomes a substantial condition on the vector potential space.

Theorem 16**.**

Let $X$ and $Y$ be compatible function spaces for the particle density and current density. Suppose furthermore that $Y\subseteq\vec{L}^{1}$. Then it follows that

[TABLE]

and from that automatically $\vec{X}\subseteq Y\subseteq Y^{**}$.

Proof.

Let $\vec{A}=(f,g,h)\in Y^{*}$ be an arbitrary vector potential. From $Y\subseteq\vec{L}^{1}$ we immediately have $Y^{*}\supseteq\vec{L}^{\infty}$ and therefore we can add a constant to one of the components of $\vec{A}$ without leaving the potential space $Y^{*}$,

[TABLE]

Next, as we assume compatibility of $X$ and $Y$, we have

[TABLE]

Furthermore, by assumption $f^{2},g^{2},h^{2}\in X^{*}$. Compatibility combined with the fact that $1\in L^{\infty}$ yields $1\in X^{*}$. Thus, the only remaining term $2f$ must be an element of $X^{*}$ too. Repeating the proof, with trivial changes for the other components, yields $f,g,h\in X^{*}$. ∎

When the preconditions of the above theorem are satisfied and $\vec{A}$ is an allowed vector potential, we thus have the peculiar situation that both $|\vec{A}|^{2}$ and $|\vec{A}|$ are contained in the function space of allowed scalar potentials.

Next, we turn to examples of reasonable choices of function spaces for paramagnetic CDFT that also illustrate the duality and compatibility conditions. Following Lieb 2, we may choose the non-reflexive space $X=L^{1}\cap L^{3}$ for the particle densities and its dual $X^{*}=L^{3/2}+L^{\infty}$ for the scalar potentials. For the current densities, we first consider the choice in the literature 31, where current densities were placed in the non-reflexive space $Y=\vec{L}^{1}$. Compatibility then follows trivially.

Proposition 17**.**

The choice $X=L^{1}\cap L^{3}$ and $Y=\vec{L}^{1}$ is compatible.

As a second example a choice of reflexive, compatible spaces should be given. Reflexivity is imperative for the construction of a well-defined KS scheme like in Sec. 3.3. To achieve this we just drop the non-reflexive $L^{1}$ from $X$ in the example above and switch to $\vec{L}^{3/2}$ instead of $\vec{L}^{1}$ for $Y$.

Proposition 18**.**

The strictly convex and reflexive choice $X=L^{3}$, $Y=\vec{L}^{3/2}$ is compatible.

Finally the admissible spaces $X=L^{1}\cap L^{3}$ and $Y=\vec{L}^{1}\cap\vec{L}^{3/2}$ that were derived from considerations regarding $N$-representability in Sec. 2.2 should be examined with respect to compatibility. It turns out that by applying Theorem 15 we can just build these spaces as intersections of the ones from Propositions 17 and 18.

Corollary 19**.**

The choice $X=L^{1}\cap L^{3}$ and $Y=\vec{L}^{1}\cap\vec{L}^{3/2}$ is compatible.

Many other compatible function spaces can be constructed. The previous examples exclude the common case of uniform magnetic fields as these require linearly growing vector potentials, e.g., $\vec{A}=\frac{1}{2}\vec{B}\times\vec{r}$, which do not belong to any $\vec{L}^{p}(\mathbb{R}^{3})$ space. We refer to Ref. 18 for a treatment of this situation. The next result provides function spaces that allow for inclusion of such vector potentials. It is a little detour to other possible choices for compatible Banach spaces, before we return to the discussion of energy functionals defined on them.

Let $w(\vec{r})$ be a suitable weight function. We use the notation $f\in L^{p}(w)$ for the weighted Lebesgue space defined by

[TABLE]

Note that $wf\in L^{p}$ is equivalent to $f\in L^{p}(w^{p})$. Some care is required however, as the two forms may produce inequivalent results when multiple weighted $L^{p}$-spaces are considered. For example, $f\in L^{p}(w^{\prime})\cap L^{q}(w^{\prime\prime})$ in general cannot be represented as $wf\in L^{p}\cap L^{q}$, for any weight function $w$. In the following examples the weight function is assumed to satisfy $w(\vec{r})\geq 1$ for all $\vec{r}\in\mathbb{R}^{3}$.

Theorem 20**.**

Let $Z$ be a normed space with dual $Z^{*}$. Then each of the following choices of function spaces is compatible,

[TABLE]

or, with $1<p<\infty$,

[TABLE]

Proof.

In both cases, we exploit the equivalence of (C1) and (C2), and only prove one of them. First, for the choice in Eq. (17), condition (C1) is trivial, since $w^{-1}\vec{A}\in\vec{L}^{\infty}$ directly implies $w^{-2}|\vec{A}|^{2}\in L^{\infty}$. Second, for Eq. (18), we similarly have that $w^{-1}\vec{A}\in\vec{L}^{2p}$ directly implies $w^{-2}|\vec{A}|^{2}\in L^{p}$, establishing (C1). ∎

In the first of the above examples, Eq. (17), we can make the trivial choice $w(\vec{r})=1$ and $Z=L^{3}$ to recover the function space analyzed in Ref. 31. Moreover, with suitable, non-trivial weight functions, unbounded vector potentials can be considered. For example, the choice $w(\vec{r})=(1+|\vec{r}|^{2})^{1/2}$ was studied extensively in Ref. 18 as it allows for uniform magnetic fields and always ensures that the angular momentum is well-defined.

2.5 Finiteness of the energy functional

The following general property of the CDFT energy functional $E^{\lambda}(v,\vec{A})$ will be useful later. It says that for both interacting and non-interacting systems the energy is finite for all considered potentials. Furthermore, if the choice of function spaces is compatible, the same boundedness from below holds for $\bar{E}^{\lambda}(u,\vec{A})$.

Lemma 21**.**

$E^{\lambda}(v,\vec{A})$* is finite for*

[TABLE]

corresponding to the density space $(\rho,\vec{j})\in(L^{1}\cap L^{3})\times(\vec{L}^{1}\cap\vec{L}^{3/2})$. By compatibility, $\bar{E}^{\lambda}(u,\vec{A})$ is also finite on the same domain.

Proof.

To prove finiteness (of the infimum) it is enough to prove boundedness from below. By definition, we have for any $(v,\vec{A})\in(L^{3/2}+L^{\infty})\times(\vec{L}^{3}+\vec{L}^{\infty})$

[TABLE]

Moreover, if $(\rho,\vec{j})$ is not $N$-representable then $F^{\lambda}(\rho,\vec{j})=+\infty$ and we can consequently consider only $N$-representable density pairs. By definition $F^{\lambda}\geq F^{0}$ and furthermore $F^{0}(\rho,\vec{j})\geq\mathrm{inf}_{\psi}\,\mathcal{K}(\psi)$ for all $\psi$ with $\rho_{\psi}=\rho$ and $\vec{j}_{\psi}=\vec{j}$. By the von Weizsäcker bound in Eq. (11), we obtain

[TABLE]

for $N$-representable density pairs $(\rho,\vec{j})$. The Sobolev inequality in Eq. (6) further allows for

[TABLE]

Combining Eqs. (19) and (20), it follows

[TABLE]

where the minimization is restricted to $N$-representable density pairs. Using the obvious inequality

[TABLE]

to replace the full square, we have

[TABLE]

Next we bound the r.h.s. of Eq. (21) from below. Since $\mathcal{C}_{0}^{\infty}$ is dense in $L^{3/2}$, also $L^{\infty}$ is. This means that we can split $v=v_{S}+v_{\infty}$ with $v_{S}\in L^{3/2}$, $v_{\infty}\in L^{\infty}$ in such a way that $\|v_{S}\|_{3/2}$ is arbitrarily small. We choose the decomposition such that $\|v_{S}\|_{3/2}\leq(2S)^{-1}$ with the constant $S$ from Eq. (7). Then using Hölder’s inequality Eq. (21) can be estimates as

[TABLE]

This proves the claim for $E^{\lambda}(v,\vec{A})$ and by compatibility also for $\bar{E}^{\lambda}(u,\vec{A})$. ∎

Remark 8*.*

The reflexive space $X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$ is a subset of the domain given in Lemma 21 and corresponds to a compatible density space by Proposition 18. It follows by an adaptation of the proof of Lemma 21 that $E^{\lambda}(v,\vec{A})$ and $\bar{E}^{\lambda}(u,\vec{A})$ are finite on $X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$. In this case when we just assume $(\rho,\vec{j})\in L^{3}\times\vec{L}^{3/2}$, however, it is crucial to exploit that the minimization can be restricted to $N$-representable density pairs, such that

[TABLE]

where $\varphi\in C_{0}^{\infty}$ has been chosen such that $v-\varphi$ has sufficiently small $L^{3/2}$ norm.

3 Regularization and the Kohn–Sham iteration scheme

In our previous work 9 the general theory of a quantum system described by the (density) variable $b\in B$, where $B$ is a reflexive ($B^{**}=B$) and strictly convex Banach space, was presented. In this theory the general problem

[TABLE]

(here $\langle b^{*},b\rangle$ denotes the dual pairing that is not necessary given by an integral) with $f\in\Gamma_{0}(B)$ is studied. We here wish to apply this structure to paramagnetic CDFT, i.e., $b=(\rho,\vec{j})$, $b^{*}=(u,\vec{A})$, and $f(b)=F^{\lambda}(\rho,\vec{j})$. We choose the density space from Definition 8

[TABLE]

that was already discussed in Proposition 18 to meet the requirements of strict convexity, reflexivity, and compatibility. The dual potential space then is

[TABLE]

Another way to gain reflexivity is by limiting the spatial domain to a bounded set $\Omega\subsetneq{\mathbb{R}}^{3}$. Then $X=L^{1}(\Omega)\cap L^{3}(\Omega),Y=\vec{L}^{1}(\Omega)\cap\vec{L}^{3/2}(\Omega)$ automatically collapse to $X=L^{3}(\Omega),Y=\vec{L}^{3/2}(\Omega)$, which are again reflexive. We also want to remark that in this case the usual Coulomb potential is included in $L^{3/2}(\Omega)$ since it is in $L_{\mathrm{loc}}^{3/2}({\mathbb{R}}^{3})$.

Furthermore, compatibility of $X_{\mathrm{R}}\times Y_{\mathrm{R}}$ gives that the concave $\bar{E}^{\lambda}$ (restricted to $X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$) can be defined. By Lemma 21 and Remark 8, this energy is also bounded below. To connect the CDFT functionals $F_{\mathrm{VR}}^{\lambda}(\rho,\vec{j})$ (or $F_{\mathrm{VR,DM}}^{\lambda}(\rho,\vec{j})$) and $\bar{E}^{\lambda}(u,\vec{A})$, we have already introduced the Legendre–Fenchel transformation in Eqs. (13) and (14). It also relates the Lieb functional $F^{\lambda}$ and the concave energy functional $\bar{E}^{\lambda}$ vice versa as a conjugate pair. Also note that $F^{\lambda}\in\Gamma_{0}(X_{\mathrm{R}}\times Y_{\mathrm{R}})$ by Lemma 12. The next step lies in another type of transformation that makes $F^{\lambda}$ functionally differentiable too, which is achieved by the Moreau–Yosida regularization.

3.1 Moreau–Yosida regularization

The original problem given in Eq. (14) of finding a ground-state density pair $(\rho,\vec{j})\in X_{\mathrm{R}}\times Y_{\mathrm{R}}$ by minimizing the convex functional $F^{\lambda}$ plus the potential energy can be restated using sub-/superdifferentials (Definition 2 in Ref. 9), both denoted $\partial$. It means selecting $(\rho,\vec{j})$ from the superdifferential of $\bar{E}^{\lambda}$ at $(u,\vec{A})\in X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$. Through the Legendre–Fenchel transformation Eq. (13), the same is possible for $F^{\lambda}$. With a minus sign in front, the potential pair $(u,\vec{A})$ yielding the ground state lies in the subdifferential of $F^{\lambda}$ at $(\rho,\vec{j})$ (see Lemmas 3 and 4 in Ref. 9),

[TABLE]

This statement can be seen as a more general reformulation of the HK theorem, but only with a $(u,\vec{A})$ potential pair, which is different from the physical $(v,\vec{A})$ setting.

The generalized notions of differentiability for convex/concave functionals involve the difficulty of non-existence or non-uniqueness. Sub- and superdifferentials are set-valued and can thus be empty or contain many elements. It is thus beneficial to “smooth out” the functional $F^{\lambda}$ in such a way that it is differentiable, which implies that the subdifferential contains only one single element. (Note that in infinite dimensions only for a continuous functional a single element in the subdifferential implies differentiability.) This is achieved by the Moreau–Yosida regularization of $F^{\lambda}\in\Gamma_{0}(X_{\mathrm{R}}\times Y_{\mathrm{R}})$, for $\varepsilon>0$ given by

[TABLE]

Since $F^{\lambda}$ is convex, the new functional $F_{\varepsilon}^{\lambda}$ is convex as well and now also functionally differentiable by Theorem 9 in Ref. 9. This regularized functional then serves as the basis for defining a new energy functional $\bar{E}^{\lambda}_{\varepsilon}$ through the Legendre–Fenchel transformation Eq. (14) again

[TABLE]

Note carefully that $\bar{E}_{\varepsilon}^{\lambda}$ is not the Moreau–Yosida regularization of some functional, but instead the Legendre–Fenchel conjugate of the regularized functional $F_{\varepsilon}^{\lambda}$. Then Theorem 10 in Ref. 9 can be used to relate the two energy functionals through

[TABLE]

Since the infimum in the definition Eq. (23) is always uniquely attained at some $(\rho_{\varepsilon},\vec{j}_{\varepsilon})\in X_{\mathrm{R}}\times Y_{\mathrm{R}}$ (see Sec. 2.2.3 in Ref. 50), we can define a mapping $(\rho,\vec{j})\mapsto(\rho_{\varepsilon},\vec{j}_{\varepsilon})$ that is called the proximal mapping, i.e.,

[TABLE]

and

[TABLE]

The proximal mapping maps density pairs that are solutions of the regularized problem $(\rho,\vec{j})\in\partial\bar{E}_{\varepsilon}^{\lambda}(u,\vec{A})$ back to solutions of the corresponding unregularized problem, $\mathrm{prox}_{\varepsilon F}(\rho,\vec{j})\in\partial\bar{E}^{\lambda}(u,\vec{A})$ by Corollary 11 in Ref. 9. Furthermore, the original functional $F^{\lambda}$ is subdifferentiable at $(\rho_{\varepsilon},\vec{j}_{\varepsilon})\in X_{\mathrm{R}}\times Y_{\mathrm{R}}$, which means that the density pair $(\rho_{\varepsilon},\vec{j}_{\varepsilon})$ is $v$-representable. We note that by Theorem 9 in Ref. 9

[TABLE]

where $\mathcal{J}:X_{\mathrm{R}}\times Y_{\mathrm{R}}\to X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$ is the duality map that is always homogeneous (Definition 7 in Ref. 9). By Proposition 1.117 in Ref. 50, it is further bijective in the present setting of reflexive and strictly convex Banach spaces (including their duals). Letting $\mathcal{J}=(\mathcal{J}_{X_{\mathrm{R}}},\mathcal{J}_{Y_{\mathrm{R}}})$, we get from $-(u,\vec{A})=\nabla F_{\varepsilon}^{\lambda}(\rho,\vec{j})$ that

[TABLE]

which straightforwardly transforms to

[TABLE]

Here the compatibility of $X_{\mathrm{R}}$ as given by Proposition 18 again becomes important since it implies that we can decompose $u$ as

[TABLE]

Recall that for all $u\in L^{3/2}$ and $\vec{A}\in\vec{L}^{3}$ we have

[TABLE]

We conclude this section by discussing a Hilbert (H) space formulation. A direct adaptation of the approach taken in Ref. 4 is to choose

[TABLE]

i.e., the Hilbert space built up from four copies of $L^{2}$. The regularization presented in Ref. 4 for standard DFT can then be directly applied to the four-vector $(\rho,\vec{j})$ instead of just the particle density $\rho$. We note that $\psi\in\mathcal{H}^{2}$ is sufficient to obtain $j^{k}\in L^{2}$, see Proposition 4. However, when $X\times Y=X_{\mathrm{H}}\times Y_{\mathrm{H}}$ the density and potential spaces are not compatible in the meaning of Definition 13. This causes a problem for the regularization procedure of the Legendre–Fenchel pair $F^{\lambda}$ and $\bar{E}^{\lambda}$, because we cannot decompose $u$ as in Eq. (27) any more. To see this, note that $|\vec{A}|^{2}$ cannot in general be assumed to satisfy $|\vec{A}|^{2}\in X_{\mathrm{H}}$, which in turn would yield the desired $v\in X_{\mathrm{H}}$. Consequently, we cannot obtain the physical setting of $E^{\lambda}(v,\vec{A})$ from the Moreau–Yosida setting of $\bar{E}^{\lambda}(u,\vec{A})$. This again highlights the usefulness of the more general reflexive Banach space formulation that allows a compatible choice of function spaces and makes a regularized paramagnetic CDFT possible.

3.2 Regularized Kohn–Sham iteration scheme in CDFT

We now revisit the KS iteration scheme that we previously analyzed for generic Banach spaces 9. Due to the similarity to the Optimal Damping Algorithm 51, 52 constructed for an unregularized setting, we baptize this iteration scheme the Moreau–Yosida Kohn–Sham Optimal Damping Algorithm (MYKSODA).

Again, let $X\times Y=X_{\mathrm{R}}\times Y_{\mathrm{R}}$ so that the space of densities is compatible, reflexive, and strictly convex. The latter two properties are also fulfilled by the dual $X^{*}\times Y^{*}=X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$. From Sec. 3.1, the ground-state problem Eq. (14) can be reformulated in terms of sub- and superdifferentials

[TABLE]

The regularized functionals $F_{\varepsilon}^{\lambda}$ and $\bar{E}_{\varepsilon}^{\lambda}$ then allow the same relation with the benefit that $F_{\varepsilon}^{\lambda}$ is now differentiable. This means we can switch from the subdifferential $\partial$ to the gradient $\nabla$ of $F_{\varepsilon}^{\lambda}$, yet this is not permitted for $\bar{E}_{\varepsilon}^{\lambda}$.

We set up two problems side by side, the interacting problem with $\lambda=1$ and the non-interacting reference problem corresponding to $\lambda=0$, i.e.,

[TABLE]

In the setting of (regularized) KS theory the external potential pair $(u_{\mathrm{ext}},\vec{A}_{\mathrm{ext}})$ is fixed and $(u_{\mathrm{KS}},\vec{A}_{\mathrm{KS}})$ is to be determined under the assumption that both problems give the same (regularized) ground-state density pair $(\rho_{\mathrm{reg}},\vec{j}_{\mathrm{reg}})$. We have highlighted the $\lambda$ dependence by using “ext” for $\lambda=1$ and by “KS” for $\lambda=0$ (but for the functionals $\bar{E}$ and $F$ we keep 1 and 0). By combining Eqs. (28) and (29) we arrive at

[TABLE]

from which the iteration scheme will be derived by replacing the unknown variables by sequences. Let $(\rho_{i},\vec{j}_{i})$ be the element of a sequence towards the (regularized) ground-state density pair $(\rho_{\mathrm{reg}},\vec{j}_{\mathrm{reg}})$, and thus the next step towards the KS potential pair $(u_{\mathrm{KS}},\vec{A}_{\mathrm{KS}})$ follows by

[TABLE]

The expression $\nabla(F_{\varepsilon}^{1}-F_{\varepsilon}^{0})$, which can be identified as a “Hartree exchange-correlation potential”, is where the usual approximation techniques of DFT enter. Yet in the domain of CDFT the variety of tried and tested functionals is meager 53, 54, 55 compared to the wealth of options in conventional DFT 56.

The second major step in the iteration scheme is then the solution of the non-interacting reference system (instead of the computationally difficult interacting problem) which selects a ground-state density pair $(\rho_{i+1},\vec{j}_{i+1})$ corresponding to the approximated KS potential pair $(u_{i+1},\vec{A}_{i+1})$ that then serves as the next input in Eq. (30)

[TABLE]

That this (super)differential is indeed always non-empty, meaning that the associated ground-state problem has at least one solution, can be shown by using the result of finiteness of $\bar{E}^{\lambda}$ from Lemma 21 (see the proof of Theorem 12 in Ref. 9).

The iteration stops when

[TABLE]

which means that $(\rho_{i},\vec{j}_{i})$ solves the original interacting ground-state problem with external potential pair $(u_{\mathrm{ext}},\vec{A}_{\mathrm{ext}})$. This is so because this condition is equivalent to

[TABLE]

by Eq. (30), and the next step given by Eq. (31) would just yield the same density pair $(\rho_{i+1},\vec{j}_{i+1})=(\rho_{i},\vec{j}_{i})$ again. In such a case, or if we decide we have converged close enough to the supposed ground-state density pair $(\rho_{\mathrm{reg}},\vec{j}_{\mathrm{reg}})$ of the regularized problem, a fixed relation between this solution and the solution of the unregularized problem is established by Eq. (26).

The question of convergence of the sequences $\{(\rho_{i},\vec{j}_{i})\}$ and $\{(u_{i},\vec{A}_{i})\}$ with respect to the Banach space topologies of $X_{\mathrm{R}}\times Y_{\mathrm{R}}$ and $X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$ is immediately raised. The authors have answered this in Ref. 9, but only in a weak sense. More precisely, the associated energy sequence

[TABLE]

can be guaranteed to converge to some value larger or equal the correct value of the regularized energy functional $\bar{E}^{1}_{\varepsilon}(u_{\mathrm{ext}},\vec{A}_{\mathrm{ext}})$. Arguably this is not what one expects from a well-formed KS iteration, where convergence to the correct ground-state density pair is the obvious aim. Further, such convergence in terms of energy can only be guaranteed if an additional step is inserted into the scheme consisting of Eqs. (30) and (31), coined “optimal damping” 51, 57, 6, that limits the step of the new density in such a way that the energy value assuredly decreases. It should be noted that a very recent development 11 treating the finite dimensional case finally proves full convergence of the exact KS iteration.

3.3 Weak-type convergence of MYKSODA

We now have all the ingredients we need for an application of Theorem 12 in Ref. 9, which constitutes the main theoretical result of this work.

Theorem 22**.**

For the density spaces $X_{\mathrm{R}}\times Y_{\mathrm{R}}=L^{3}\times\vec{L}^{3/2}$, and the corresponding dual for potential pairs $X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}=L^{3/2}\times\vec{L}^{3}$, a well-defined KS iteration can be set up for the energy functional $\bar{E}^{\lambda}$ from Eq. (2). It starts with a fixed potential pair $(u_{\mathrm{ext}},\vec{A}_{\mathrm{ext}})\in X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$ by setting $(u_{1},\vec{A}_{1})=(u_{\mathrm{ext}},\vec{A}_{\mathrm{ext}})$ and selecting $(\rho_{1},\vec{j}_{1})\in\partial\bar{E}^{0}_{\varepsilon}(u_{\mathrm{ext}},\vec{A}_{\mathrm{ext}})$. Then iterate $i=1,2,\ldots$ according to:

1. (a)

Set

[TABLE]

and stop if

[TABLE] 2. (b)

Select $(\rho_{i+1}^{\prime},\vec{j}_{i+1}^{\prime})\in\partial\bar{E}_{\varepsilon}^{0}(u_{i+1},\vec{A}_{i+1})$. 3. (c)

Choose $t_{i}\in(0,1]$ maximally such that for

[TABLE]

one still has

[TABLE]

Then the strictly descending sequence

[TABLE]

converges as a sequence of real numbers to

[TABLE]

Thus,

[TABLE]

is an upper bound for the ground-state energy $\bar{E}^{1}(u_{\mathrm{ext}},\vec{A}_{\mathrm{ext}})$.

Proof.

To be able to apply Theorem 12 from Ref. 9 we have to make sure that $X_{\mathrm{R}}\times Y_{\mathrm{R}}$ and $X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$ are reflexive and strictly convex, the non-interacting energy functional $\bar{E}^{0}$ needs to be finite on all of $X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$, and $F^{\lambda},\bar{E}^{\lambda}$ must form a convex-concave pair linked by the Legendre–Fenchel transformation. Now, Proposition 18 shows that the chosen density spaces are indeed reflexive, all $L^{p}$ with $1<p<\infty$ are strictly convex anyway, and also that they are compatible. Compatibility gives that the energy functional can be transformed to a concave $\bar{E}^{\lambda}$ by Eq. (12) that links to a convex Lieb functional $F^{\lambda}$ by Eq. (13) (Legendre–Fenchel transformation). Then Lemma 21 and Remark 8 prove that $\bar{E}^{0}$ is indeed finite on $X_{\mathrm{R}}^{*}\times Y_{\mathrm{R}}^{*}$. With the results from Theorem 12 in Ref. 9 we get a strictly decreasing and converging sequence

[TABLE]

with the given lower bound. The transformed energy bound follows directly from Eq. (24). ∎

Remark 9*.*

Any candidate for a possible ground-state density pair from the iteration defined in Theorem 22 can be transformed to a solution of the corresponding “physical” unregularized problem with the help of Eq. (26). But it has not been proven that the iteration actually converges in terms of densities and potentials as elements of the given Banach spaces and dual spaces or that if it converges, it actually reaches the ground-state density pair $(\rho_{\mathrm{reg}},\vec{j}_{\mathrm{reg}})$ and the associated KS potential pair $(u_{\mathrm{KS}},\vec{A}_{\mathrm{KS}})$.

Remark 10*.*

We already discussed $(u_{\mathrm{ext}},\vec{A}_{\mathrm{ext}})=-\nabla F_{\varepsilon}^{1}(\rho_{i},\vec{j}_{i})$ as a stopping condition for the iteration before. To have such a stopping condition is important for at least a possible convergence to the correct ground-state density pair that is then a fixed point. Note that we still have the appearance of a whole set of ground-state density pairs in step (b), signifying that degeneracy is admitted. Thus, it is more beneficial to look at the sequence of potential pairs $(u_{i},\vec{A}_{i})$ that, if the stopping condition is eventually reached, gives the correct KS potential for some ground-state densities. This is a reason why it was important to switch to differentiable functionals $F_{\varepsilon}^{\lambda}$ through regularization: To be able to define a unique sequence of potentials that can converge to the KS potential. See also Ref. 10 where the traditional iteration in density space is supplemented with a bivariate formalism.

Remark 11*.*

Note that step (c) in the KS iteration scheme above corresponds to a line-search between the points $(\rho_{i},\vec{j}_{i})$ and $(\rho_{i+1}^{\prime},\vec{j}_{i+1}^{\prime})$. If more solutions of the regularized, non-interacting reference system from step (b) are taken into account (degeneracy), then step (c) gets generalized to a search over a convex polytope.

3.4 Kohn–Sham iteration scheme for two-electron systems

As only a partial convergence result is presently available for the KS algorithm—leaving open the possibility that it does not always converge to the right energy and potential—it is interesting to consider a case where non-interacting $N$- and $v$-representability issues pose a challenge for the algorithm. Consider a formulation of paramagnetic CDFT for singlet ground states and restrict attention to a two-electron system. Then the (unregularized) non-interacting KS system is represented by a single orbital and its vorticity,

[TABLE]

vanishes if differentiability is assumed. For the small set of KS potentials that yield ground-state degeneracies, this can be circumvented by allowing the KS system to be represented by a mixed state. However, for most KS potentials, the KS ground state is unique and has trivial vorticity. On the other hand, for most external potentials, the correlated ground state of the interacting system has a nontrivial vorticity. Hence, most of the ground-state densities are not non-interacting ($\lambda=0$) $N$-representable.

This situation poses an interesting challenge for the KS iteration scheme as the vorticity of the KS system cannot develop gradually. Until the algorithm has constructed potentials that yield an exact ground-state degeneracy, the vorticity is trivial, and it is not clear how “visible” the corresponding degrees of freedom are to the optimization algorithm. Moreau–Yosida regularization alleviates the challenge somewhat, since the relevant densities have contributions from the potentials, complicating the non-interacting $N$-representability conditions. With regularization, the relevant vorticity that should reproduce the interacting system using Eq. (26) is

[TABLE]

If counterexamples that prevent a full convergence proof exist at all, the type of system sketched above is a promising candidate for further analysis.

4 Numerical application to quantum ring

The theory of regularized CDFT can be directly applied to a one-dimensional quantum ring. Although this is a toy model, it is sufficiently rich to contain simple formal analogues of many aspects of CDFT for a three-dimensional spatial domain. We limit attention to two-electron systems in singlet spin states. The Hamiltonian is given by

[TABLE]

where $R$ is the radius of the ring and the potentials $v$ and $A$ as well as the electron-electron repulsion $W$ are considered functions of the angular position along the ring. Note that gradients and the vector potential only have tangential components and may therefore effectively be treated as scalars.

Because of the limitation to singlet states, the spatial wave function $\psi(\theta_{1},\theta_{2})=\psi(\theta_{2},\theta_{1})$ must be symmetric. Any uncorrelated state, e.g., a KS state, takes the form $\psi(\theta_{1},\theta_{2})=\phi(\theta_{1})\phi(\theta_{2})$ and is defined by a single orbital $\phi$. The densities that arise from such an uncorrelated, single-orbital $\phi(\theta)=\sqrt{\rho(\theta)}\,e^{i\chi(\theta)}$ state must satisfy

[TABLE]

with $m$ an integer if $\rho>0$ everywhere. Note that $\chi(\theta+2\pi)=\chi(\theta)+2\pi m$ is in general a multivalued phase function. By contrast, a correlated state can well give rise to a fractional value of $m$. This is the quantum ring analogue of the fact that vorticity is trivial for single-orbital systems in a three-dimensional spatial domain.

Next, in order to study regularized CDFT numerically, we discretize the quantum ring into $N_{G}$ uniformly spaced grid points. The approach described below is implemented in a Matlab program named MYring 58. We replace the Laplacian by the standard second-order finite difference expression

[TABLE]

where $h=2\pi R/N_{G}$ is the grid spacing. The paramagnetic term is discretized using the symmetric first-order expression,

[TABLE]

Defining the particle density and current density at grid point $k$ by

[TABLE]

we can define a constrained-search functional as well as linear pairings between $\rho$ and $u=v+A^{2}/2$ as well as $j$ and $A$,

[TABLE]

The densities and potentials may all be regarded as vectors in $\mathbb{R}^{N_{G}}$. Because all norms in finite dimensions are mathematically equivalent, we can choose to endow all function spaces with the same Euclidean $l^{2}(N_{G})$ norm without losing compatibility. However, the norms may not be numerically equivalent. Moreover, to connect to the continuum limit when $N_{G}\to\infty$, it is likely that one needs more carefully chosen norms. This is left for future studies.

The grid discretization makes it trivial to construct compatible finite-dimensional function spaces. However, this is not true in arbitrary basis expansions of density pairs and potential pairs. Unless the respective basis sets have special properties, compatibility is in general lost.

By solving the discretized Schrödinger equation, we then obtain the ground-state energy and the regularized energy

[TABLE]

The universal density functional can be computed from the Lieb variational principle

[TABLE]

with $G^{\lambda}_{\varepsilon}(u,A;\rho,j)=\bar{E}^{\lambda}_{\varepsilon}(u,A)-\langle u,\rho\rangle-\langle A,j\rangle$. We have found that a cutting-plane bundle method for convex optimization 59, 60 provides robust, though occasionally very slow convergence to the maximum value. In more detail, our implemented method maintains a “bundle” of data $(g_{l},u_{l},A_{l},\sigma_{l},k_{l})$ from previous iterations. The bundle contains the function value $g_{l}=G^{\lambda}_{\varepsilon}(u_{l},A_{l};\rho,j)$ and a supergradient $(\sigma_{l},k_{l})$ evaluated at $(u_{l},A_{l})$. Then a model function is defined by all the tangent planes encoded in the bundle,

[TABLE]

The next sample point $(u_{l+1},A_{l+1})$ is determined by maximizing the model function subject to a trust region constraint (to guard against $Q_{l}$ being unbounded, as may happen in the first iterations).

The stopping criterion requires care. As mentioned above, because the Moreau–Yosida regularization is only applied to $F^{\lambda}$, the energy functional $\bar{E}^{\lambda}_{\varepsilon}$ is not more differentiable than the original $\bar{E}^{\lambda}$. Hence, there is no guarantee that $G^{\lambda}_{\varepsilon}$ is differentiable with respect to the potentials at the maximum. This is particularly true for the KS potentials at $\lambda=0$, where $N$-representability constraints become more visible. Hence, it is not feasible to rely on vanishing (super)gradients as a stopping criterion for the optimization. This is connected to ground-state degeneracy and can be diagnosed by computing the energy gap to the first excited state.

4.1 Kohn–Sham potentials from the Lieb variational principle

We consider a discretization with $N_{G}=30$ grid points and set the electron-electron interaction to

[TABLE]

We choose the external potentials

[TABLE]

in order to obtain a non-trivial example that is nonetheless simple to specify. The external potentials are visualized in Fig. 1. The resulting Hamiltonian $H^{1}(v_{\mathrm{ext}},A_{\mathrm{ext}})$ has a highly correlated ground state with densities $(\rho,j)$ displayed in Fig. 2. Performing maximization in the Lieb variation principle (Eq. (13) or (33)) defining $F_{\varepsilon=0}^{\lambda}(\rho,j)$ yields KS potentials $(u_{\mathrm{KS}},A_{\mathrm{KS}})$ as a by-product, visualized in Fig. 1. The Hamiltonian $H^{0}(v_{\mathrm{KS}},A_{\mathrm{KS}})$, with $v_{\mathrm{KS}}=u_{\mathrm{KS}}-(A_{\mathrm{KS}})^{2}/2$, has a two-fold ground-state degeneracy and one of these ground-state densities is shown in Fig. 2. The vanishing gap is seen in Fig. 3 and results in a non-differentiable kink in the ground-state energy $\bar{E}_{\varepsilon=0}^{0}(u_{\mathrm{KS}},A_{\mathrm{KS}})$. The interacting density pair $(\rho,j)$ is a supergradient at this non-differentiable point, but it is neither a left- nor a right-derivative. Due to the limitation that our implementation is limited to pure states, and furthermore that the the choice of degenerate eigenvector basis is not optimized, it is seen in Fig. 2 that the interacting ground-state density pair $(\rho,j)$ is not reproduced exactly by the KS ground state. In general, exact reproduction requires mixed states.

Next we illustrate the regularized setting by taking $\varepsilon=0.1$. This relatively large regularization parameter is used to make the effects of regularization noticeable. It is now the pair $(\sigma=\rho-\varepsilon u_{\mathrm{ext}},k=j-\varepsilon A_{\mathrm{ext}})$ that takes over the role played by the density pair in the unregularized setting. In particular, the Lieb variation principle now yields a KS potential pair $(u_{\mathrm{KS}},A_{\mathrm{KS}})$ such that $(\sigma_{\mathrm{KS}}=\rho_{\mathrm{KS}}-\varepsilon u_{\mathrm{KS}},k_{\mathrm{KS}}=j_{\mathrm{KS}}-\varepsilon A_{\mathrm{KS}})$ coincides with the density pair $(\sigma,k)$, but $(\rho_{\mathrm{KS}},j_{\mathrm{KS}})\neq(\rho,j)$. Hence, the KS potentials shown in Fig. 4 are different from those in the unregularized setting (Fig. 1). The resulting densities are shown in Fig. 5.

4.2 Kohn–Sham potentials from the iterative algorithm

In the previous section, the KS potentials were determined by first solving for correlated ground-state wave function of the interacting system, then constructing its densities, and finally plugging these densities into the Lieb variation principle. The iterative KS algorithm discussed in Sec. 3.2 above provides an alternative that does not require any a priori information of the correlated ground state or its associated density. For simplicity, we have implemented the pure-state version of this algorithm, enabling us to see the consequences when a density pair is not representable by a pure ground state. The linesearch for the interpolation parameter $t$ was implemented in the following way:

• (i)

Successively try $t=1,\tfrac{1}{2},\tfrac{1}{4},\tfrac{1}{8},\ldots$, until the criterion from the optimal damping step Eq. (32) is fulfilled.

• (ii)

If already $t=1$ fulfills Eq. (32), then use this value. Otherwise, let $t=2^{-k}$ be the first parameter value such that Eq. (32) holds and estimate the critical $t$ value by linear interpolation between $t=2^{-k+1}$ and $2^{-k}$. If the criterion is still not fulfilled at this $t$, perform another linear interpolation and choose the best of the sampled values.

The computation of gradients $\nabla F^{1}_{\varepsilon}(\rho,j)$ is done using the Lieb variation principle, with a maximum of 300 bundle optimization iterations and a convergence criterion of $10^{-5}$ for stopping earlier. When there is a degenerate ground state, the gradient criterion does not apply and we instead test for a small gap and stagnated bundle iterations. In cases of numerically very small, but non-zero gap between the ground state and first excited state, our implementation may fail to obtain an adequate solution from the (pure-state) Lieb variation principle. As the algorithm was not formulated to account for such failures, the energy seen in the KS iterations may not be bounded from below by the true energy $\bar{E}^{1}(u_{\mathrm{ext}},A_{\mathrm{ext}})$, unless we override the reference density pair and instead use the actual density pair returned from the Lieb optimization.

Continuing with the same numerical example as in the previous section, we ran the KS iteration for different values of the regularization parameter. In the unregularized case, the consequences of failure of pure-state representability, both for the KS system and interacting systems corresponding to trial densities encountered in the course of the iterations, prevented a meaningful result. With Moreau–Yosida regularization, we were able to converge within the expected accuracy, given the finite precision of our implementation of the Lieb variation principle. In Fig. 6 the convergence of the energy difference,

[TABLE]

is shown for four different values, $\varepsilon=0.05,0.1,0.2,0.3$, of the regularization parameter. Fig. 7 shows the convergence of the gradient norm,

[TABLE]

which vanishes when the ground-state density of the interacting system has been reproduced. Although not encountered in the example studied here, small numerical inaccuracies especially in the Lieb variation principle lead to occasional small increases of the energy. The convergence is slow compared to experience with standard algorithms, such as Pulay’s DIIS 61, and approximate density functionals, as these result in quadratic convergence in favorable cases. However, most standard algorithms also lack formal convergence guarantees and have, for practical reasons, never been tested with the exact functional. An exception is the work by Wagner et al. that did explore convergence of the exact functional using an algorithm applied to one-dimensional systems 8, 44. In Ref. 44 an adaptive choice of the damping (mixing) parameter was investigated, including discussions on line search and Hermite spline fit to the energy as a function of the damping parameter. (It is interesting to note that they use the curvature of the energy as information. In the regularized setting where derivatives are guaranteed to exist, the curvature is a key ingredient in the convergence proof of Ref. 11.) Furthermore, their study of an optimal damping parameter demonstrated numerically that convergence is more difficult for strongly correlated systems.

The present work is the first time a KS vector potential, corresponding to an exact CDFT functional, is calculated using a KS iteration scheme. As expected, the convergence in Figs. 6 and 7 is faster for larger values of the regularization parameter. This is partly due to the fact that the unregularized case features a KS system with vanishing gap and partly due to the increased regularity of the problem for larger $\varepsilon$.

Finally, Fig. 8 shows the calculated values of the damping parameter $t$ as a function of iteration number. The parameter values vary substantially between the examples with different $\varepsilon$ and also from one iteration to the next. In particular, for $\varepsilon=0.1$ the MYKSODA iterations alternate between smaller $t\sim 0.05$ and larger values in the range $0.1<t<0.65$. A simpler iterative algorithm could use a fixed $t$ in all iterations, as was done in Ref. 8. To explore this possibility, we fixed a conservative value $t=0.05$ for the damping parameter. As seen in Fig. 9 this yields dramatically slower convergence, showing that $t$ in general needs be chosen adaptively.

5 Conclusions

We have given a comprehensive account of the rigorous formulation of Kohn–Sham theory for CDFT. An important point is that textbook treatments of DFT rely on ill-defined functional derivatives 32. However, recent work has demonstrated that functional derivatives can be made well-defined and rigorous using Moreau–Yosida regularization 4, 9. We have extended that approach to functional differentiation in CDFT, enabling us to obtain well-defined Kohn–Sham potentials as well as an iteration scheme (MYKSODA). The presented MYKSODA is an algorithm for practical calculations in the setting of ground-state CDFT within a regularized framework. A toy model in the form of a quantum ring is solved numerically, and allowed a study of MYKSODA for the exact universal density functional. The calculations illustrate the performance of the algorithm and highlight the difference to iteration schemes with a constant damping factor. It is also the first implementation of a Moreau–Yosida regularized Kohn–Sham approach.

While our model was solved numerically with the exact functional, this is of course not feasible for more realistic settings where we must resort to density-functional approximations. This raises the question of how to develop such approximations for the Moreau–Yosida regularized setting, or alternatively, of how to compute the Moreau–Yosida regularization of well-established density-functional approximations. This is an interesting topic for future investigation.

Central to the theory developed here was the concept of compatibility of spaces of densities and current densities. It allows a fully convex formulation of the theory and demands the use of Banach spaces for the basic variables. The respective $\vec{L}^{p}$ constraints for current densities were determined optimally in order to complement knowledge from traditional DFT and previous work on CDFT. This article sets the stage for further inquiries into the field, such as the possible full convergence of the iteration scheme and the study of approximate (regularized) functionals for CDFT.

Acknowledgments

We thank an anonymous referee for improvements on our proof of Lemma 21. This work was supported by the Norwegian Research Council through the CoE Hylleraas Centre for Quantum Molecular Sciences Grant No. 262695. AL is grateful for the hospitality received at the Max Planck Institute for the Structure and Dynamics of Matter in Hamburg, while visiting MP and MR. MP acknowledges support by the Erwin Schrödinger Fellowship J 4107-N27 of the FWF (Austrian Science Fund) and is thankful for an invitation to the Hylleraas Centre just taking place writing this. AL and SK were supported by ERC-STG-2014 under grant agreement No. 639508. EIT was supported by the Norwegian Research Council through Grant No. 240674.

Appendix A A theorem on everywhere defined functionals on spaces of measurable functions

On any infinite-dimensional Banach space (assuming the axiom of choice) there exist everywhere defined linear maps that are unbounded. The following theorem shows that this cannot happen for linear functionals on spaces of measurable functions that are defined as integrals. The proof is based on a construction by D. Fischer 62.

Theorem 23**.**

Let $B$ be a Banach space consisting of measurable functions $f:\mathbb{R}^{n}\to\mathbb{R}$. Let $g$ be a measurable function. Then the functional $T:f\mapsto\int gf\,\mathrm{d}\mu$ is in $B^{*}$ if and only if for all $f\in B$,

[TABLE]

Proof.

Since a bounded linear functional must be everywhere defined, the only if part is trivial. Suppose $g$ is measurable and that the integral $\int fg\,\,\mathrm{d}\mu$ exists for all $f\in B$. For $n\in\mathbb{N}$, define a sequence of bounded functions with bounded support,

[TABLE]

Then $g_{n}$ is measurable for all $n$, and $h_{n}(x)=g_{n}(x)f(x)\to h(x)=g(x)f(x)$ for all $x$. Moreover $|g_{n}(x)f(x)|\leq|h(x)|$ for all $x$, the latter function being integrable by assumption. By the dominated convergence theorem,

[TABLE]

as $n\to+\infty$. Thus, the family of continuous linear functionals $T_{n}\colon f\mapsto\int fg_{n}\,\mathrm{d}\mu$ is pointwise bounded.

The uniform boundedness principle states that a family $\{T_{n}\}$ of pointwise bounded linear functionals is in fact uniformly bounded. Thus, $\sup_{n}\|T_{n}\|_{B^{*}}<+\infty$. It then follows that

[TABLE]

Hence, $T\in B^{*}$. ∎