Lower Semicontinuity of the Universal Functional in Paramagnetic Current-Density Functional Theory.

A cornerstone of current-density functional theory (CDFT) in its paramagnetic formulation is proven. After a brief outline of the mathematical structure of CDFT, the lower semicontinuity and expectation-valuedness of the CDFT constrained-search functional is proven, meaning that there is always a minimizing density matrix in the CDFT constrained-search universal density functional. These results place the mathematical framework of CDFT on the same footing as that of standard DFT.

Density functional theory (DFT) is at present the most widely used tool for first-principles electronic structure calculations in solid-state physics and quantum chemistry. DFT was put on a solid mathematical ground by Lieb in a landmark paper[1] from 1983, where he introduced the universal density functional F(ρ) as the convex conjugate to the concave ground-state energy E(v) for an electronic system in the external scalar potential v.

For electronic systems under the influence of a classical external magnetic potential A, current–density functional theory (CDFT) was introduced by Vignale and Rasolt in 1987.[2] In addition to the density ρ, the paramagnetic current density jp becomes a basic variable. The mathematical foundation of CDFT was put in place by Tellgren et al.[3] and Laestadius[4,5] in the 2010s on the basis of Lieb’s treatment of the field-free standard case. However, a central piece of the puzzle has been missing—namely, whether the CDFT constrained-search functional F(ρ, jp) is lower-semicontinuous and expectation-valued,[6] i.e., whether the infimum in its definition (see eq below) is in fact attained. These foundational issues are important because CDFT is the natural extension of DFT to treat general magnetic systems and several numerical implementations have been reported, although the development of practical functionals lags behind standard DFT.[7−13]

In this Letter, we provide proofs of the above assertions. The CDFT constrained-search functional is indeed convex lower-semicontinuous and can therefore be identified with the CDFT Lieb functional—that is, the Legendre–Fenchel transform of the energy. Without this fact, the ground-state energy functional E(v, A) and the constrained-search functional F(ρ, jp) contain different information. If F(ρ, jp) were not expectation-valued, one would lose the interpretation of the universal functional as intrinsic energy, which is very useful in standard DFT. For the interested reader, suggested further reading for convex analysis are van Tiel’s excellent introductory text[14] and the monograph by Ekeland and Témam.[15] Also, the monograph by Barbu and Precupanu,[16] which treats convex analysis in Banach spaces, and the one by Bauschke and Combettes,[17] which focuses on the Hilbert space formulation, are highly recommended. For more details on trace-class operators, the monograph by Weidmann is an accessible starting point,[18] as well as the now classic volume by Reed and Simon.[19]

For an N-electron system in sufficiently regular external potentials v and A, the ground-state energy is given by the Rayleigh–Ritz variation principle aswhere is the electronic Hamiltonian with kinetic energy operator T(A) = ∑[−∇ + A(r)]2 and two-electron repulsion operator W. The minimization is over all N-electron density matrices Γ of finite kinetic energy, for which the one-electron density is and an element jp ∈ Xp = L1 ∩ L3/2.[20] (The boldface notation indicates a space of vector fields.) The external potential energy (v|ρ) = , the paramagnetic and diamagnetic terms and and thus the Hamiltonian H(v, A) are well-defined for any v ∈ = L3/2 + and A ∈ = L3 + L∞, where XL′ and Xp′ are the dual spaces of XL and Xp, respectively. Examples of such potentials are the nuclear Coulomb potentials and uniform magnetic fields inside bounded domains. The symbol XL for the space of densities is so chosen to indicate that it is the density space of Lieb’s analysis, while Xp indicates “paramagnetic” current densities.

By a well-known reformulation of eq , we obtain the CDFT Hohenberg–Kohn variation principle:Here the Vignale–Rasolt constrained-search density functional F: XL × Xp → [0, +∞] is defined bywhere H0 = T(0) + W is the intrinsic electronic Hamiltonian, and Γ → (ρ, jp) means that the infimum is taken over all N-electron density matrices Γ with density–current pair . Thus, if (ρ, jp) is not N-representable, we have F(ρ, jp) = +∞. The universal density functional F is the central quantity in any flavor of DFT, and its mathematical properties and approximation are of utmost importance to the field.

Although E in eq is not concave, it is readily seen that the reparametrized energyis concave. This reparametrization relies on a technical notion of compatibility of function spaces for the scalar and vector potentials,[20] which is satisfied for the potentials we consider here.

From the concavity and upper semicontinuity of the modified ground-state energy Ẽ, one can deduce the existence of an alternative universal density functional F̂: XL × Xp → [0, +∞] that is related to the ground-state energy by Legendre–Fenchel transformations in the mannerwhere the optimizations are over the space XL × Xp and its dual XL′ × Xp′, respectively. As a Legendre–Fenchel transform, the functional F̂ is convex and lower-semicontinuous. In this formulation of CDFT, the ground-state energy Ẽ and the universal density functional F̂ contain precisely the same information: each functional can be obtained from the other and therefore contains the same of the information about ground-state electronic systems in external scalar and vector fields.

From a comparison of the Hohenberg–Kohn variation principles in eqs and 5, it is tempting to conclude that F̂ and F are the same functional (i.e., F̂ = F), producing the same ground-state energy for each (v, A). However, there exist infinitely many functionals F̃: XL × Xp → [0, +∞] that give the correct ground-state energy E(v, A) (but not necessarily the same minimizing densities, if any) for each (v, A) in the Hohenberg–Kohn variation principle. Each such F̃ is said to be an admissible density functional.[6] Among these, the functional F̂ stands out as being the only lower-semicontinuous and convex universal density functional and a lower bound to all other admissible density functionals (i.e., F̂ ≤ F̃). The functional F̂, the closed convex hull of all admissible density functionals, is thus the most well-behaved admissible density functional. Indeed, we may view it as a regularization of all admissible density functionals, known as the Γ regularization in convex analysis. (This name is unrelated to our notation of density matrices.)

A fundamental result of Lieb’s analysis of DFT is the identification of the transparent constrained-search density functional with the mathematically well-behaved closed convex hull F̂. The identification follows since F is convex and lower-semicontinuous. Whereas convexity follows easily for the CDFT Vignale–Rasolt functional F, the proof of lower semicontinuity is nontrivial. For standard DFT it is given in ref [1], and for CDFT it is provided in the present Letter.

We simplify our analysis by merely assuming that the density–current pairs are (ρ,jp) ∈ L1 × L1 = which we denote as X. With this topology, the potentials must be taken to be bounded functions, (v, A) ∈ X′ = L∞ × L∞ = [L∞ This simplification is irrelevant in this context: if F can be shown to be lower-semicontinuous in the topology, it will be lower-semicontinuous in any stronger topology, as required if we enlarge the potential space to include more singular functions such as those in XL′ × Xp′. Indeed, the original proof of lower semicontinuity of the standard DFT Levy–Lieb functional (eq ) was with respect to the topology, from which the same property with respect to the XL topology immediately follows.

Theorem and Proof. The intrinsic Hamiltonian H0 is self-adjoint (H0 = H0†) over L2, the Hilbert space of square-integrable N-electron wave functions (with spin and permutational antisymmetry built in). The expectation values of H0 and H(v, A) are well-defined on the Sobolev space H1, the subset of L2 with finite kinetic energy.

We denote by the convex set of N-electron mixed states with finite kinetic energy. We have the mathematical characterization[21]where TC(L2) is the set of trace-class operators over L2, the largest set of operators to which a basis-independent trace can be assigned. An operator A is trace-class if and only if the positive square root is trace-class.[18,19] A self-adjoint operator A is trace-class if and only if it has a spectral decomposition of the form , where {ϕ} forms an orthonormal basis and ∑λ is absolutely convergent. Now if and only if λ ≥ 0, ∑λ = 1, and {ϕ} ⊂ H1 and if the total kinetic energy is finite (i.e., ∑λ⟨ϕ|T|ϕ⟩ < +∞).

For any ψ ∈ H1, the density–current pair (ρ,jp) ∈ L1 × L1 is defined bywhere we integrate over all spin variables and over N – 1 spatial coordinates, τ–1 = (σ1, x2, ..., x). For , we can for instance compute ρ = ρΓ from ∑λρ, where ρ is obtained from eq with ψ = ϕ (and similarly for jp).

The theorem involves the weak topology on . Weak convergence of a sequence {x}⊂ X, written as x ⇀ x ∈ X, means that for any bounded linear functional ω ∈ X′, we have ω(x) → ω(x) as a sequence of numbers—that is, weak convergence is the pointwise convergence of all bounded linear functionals. Recall that the dual space of is , so if and only if (f|ρ) → (f|ρ) for every . Likewise, (ρ, jp) ⇀ (ρ, jp) ∈ X if and only if (f|ρ) → (f|ρ) and (a|jp) → (a|jp) for every (f, a) ∈ X′.

The trace-class operators over a separable Hilbert space are examples of compact operators, an infinite-dimensional generalization of finite-rank operators. Indeed, the set of compact operators is the closure of the set of finite-rank operators in the norm topology and thus a Banach space. The dual space of is in fact . For and , the dual pairing is Tr(BA). Similar to the weak topology for a Banach space, the dual of a Banach space can be equipped with the weak-* topology. A sequence of trace-class operators {A} converges weak-* to if, for each , Tr(BA) → Tr(BA).

We now state and prove our main result, from which lower semicontinuity follows in Corollary 1 and expectation-valuedness in Corollary 2. The theorem is the CDFT analogue of Theorem 4.4 in ref (1).

Theorem 1.Suppose that (ρ, jp) ∈ Xand {(ρ, jp)} ⊂ Xare such thatF(ρ, jp) < +∞ andF(ρ, jp) < +∞ for eachand further suppose that (ρ, jp) ⇀ (ρ, jp). Then there existssuch that Γ → (ρ, jp) and Tr(H0Γ) ≤ lim infF(ρ, jp).

Proof of Theorem 1. The initial setup follows ref (1), which we here restate. Without loss of generality, we may replace H0 = T + W by h2 = T + W + 1, which is self-adjoint and positive-definite. The operator h is taken to be the unique positive self-adjoint square root of T + W + 1.

Consider the sequence {g} with elements g ≔ F(ρ, jp). If g → +∞, then the statement of the theorem is trivially true. We therefore assume that {g} is bounded. There then exists a subsequence such that g ≔ lim g exists. Furthermore, for each n there exists such that Γ → (ρ, jp) and Tr(hΓh) = Tr(h2Γ) ≤ g + 1/n. To demonstrate this, we select for each n a density matrix Γ → (ρ, jp) that satisfies Tr(h2Γ) < g + 1/2n and choose m such that |g – g| < 1/2n for each n > m (by taking a subsequence if necessary); for each n > m, we then haveUsing the sequence {hΓh}, we next establish a candidate limit density operator .

The dual-space sequence of (positive-semidefinite) operators y ≔ hΓh ∈ TC(L2) is uniformly bounded in the trace norm: ∥y∥TC ≤ g + 1. By the Banach–Alaoglu theorem, a norm-closed ball of finite radius in the dual space is compact in the weak-* topology. Thus, there exists y ∈ TC(L2) such that for a subsequence, Tr(By) → Tr(By) for each , meaning that y is the (possibly nonunique) weak-* limit of a subsequence of {y}. The limit is positive-definite, since the orthogonal projector PΦ onto is a compact operator, which gives

We now define Γ = h–1yh–1, which fulfills all of the criteria for being an element of , except possibly TrΓ = 1, although TrΓ ≤ 1 is already implied by the weak convergence. (It should be noted that Γ has finite kinetic energy since Tr(h2Γ) < +∞.) If we can show that Γ → (ρ, jp), then we are done with the complete proof since follows from and since

Let (ρ′, jp′) ← Γ be the density associated with Γ. To demonstrate that (ρ′, jp′) = (ρ, jp), we recall that (ρ, jp) ⇀ (ρ, jp) by assumption. Since weak limits are unique, our proof is complete if we can show that (ρ, jp) ⇀ (ρ′, jp′) in . The proof that ρ ⇀ ρ′ is given in ref (1). Here we demonstrate that jp ⇀ jp′ by showing for each that (jp – jp′ | a) → 0.

Let be a bounded domain with characteristic function χ equal to 1 on Ω and 0 elsewhere. Since , for a given ε > 0 we may choose Ω to be sufficiently large that ∫(1 – χ)ρ dr < ε and ∫(1 – χ)ρ′ dr < ε. Since ρ ⇀ ρ, we also have ∫(1 – χ)(ρ – ρ) dr < ε for sufficiently large n. From the triangle inequality, we obtain ∫(1 – χ)ρ dr ≤ ∫ (1 – χ)(ρ – ρ) dr + ∫ (1 – χ)ρ′ dr, implying that ∫(1 – χ)ρ dr < 2ε for sufficiently large n.

In the notation τ = (r1, τ–1) = (x1, x2, ..., x) and τ–1 = (σ1, x2, ..., x) with space–spin coordinates x = (r, σ), we definewhere α denotes a Cartesian component and we have introduced the spectral decomposition Γ = ∑μλμ|ψμ⟩⟨ψμ| ∈ with ψμ ∈ H1. We note that if Γ → (ρΓ, jpΓ), then integration of Uα over τ–1 gives the current component jpαΓ(r) = ∫Uα(r, τ–1) dτ–1.

We now let S = ∏χ(ri) be the characteristic function of . By the definition of Uα, we then haveApplying the Cauchy–Schwarz inequality twice, we obtainNoting that 1 – S ≤ ∑[1 – χ(ri)] and using the symmetry of |ψμ|2, we obtain for the two factorsandWe conclude that I(Uα)2 ≤ 2Ngε. Introducing U = N Im diag ∂1αΓ and proceeding in the same manner, we arrive at the bound I(U)2 ≤ 4Ngε, assuming that n has been chosen to be sufficiently large that ∫(1 – χ)ρ dr < 2ε holds.

We are now ready to consider the weak convergence jp ⇀ jp′ in . For each and for sufficiently large n, using the Cauchy–Schwarz inequality and the Hölder inequality in combination with the bounds I(Uα)2 ≤ 2Ngε and I(U)2 ≤ 4Ngε, we obtain the inequalitySince ε > 0 is arbitrary, it only remains to show that we have ∫(U – Uα)aα(r1)S dτ → 0 as n → ∞.

Let M be the compact multiplication operator associated with aα(r1)S(τ), a bounded function with compact support over . Let Ωσ = {↑, ↓} be the set consisting of the two spin states of the electrons. We note thatviewing Γ as an operator over L2((Ω×Ωσ)) by domain restriction of the spectral decomposition elements—that is, ψμ ∈ H1((Ω×Ωσ)), meaning that the 2 spin components of ψμ are in H1(Ω). For simplicity, the spaces used here are not antisymmetrized.

Our next task is to demonstrate that B = h–1M∂1αh–1 is compact over L2((Ω×Ωσ)). We first show that h–1 is compact with range H1((Ω×Ωσ)). We have with domain H1((Ω×Ωσ)). Now h–1 exists and is bounded since −1 is not in the spectrum of T + W—that is, h–1: L2((Ω×Ωσ)) → H1((Ω×Ωσ)) is bounded. By the Rellich–Kondrachov theorem, H1(Ω) (the standard Sobolev space without spin) is a compact subset of L2(Ω). It follows that H1((Ω×Ωσ)) is a compact subset of L2((Ω×Ωσ)), since the tensor product of compact sets is compact. Hence, h–1 is compact.

Next, the operator ∂1α is, by the definition of the Sobolev space H1(Ω), bounded from H1((Ω×Ωσ)) to L2((Ω×Ωσ)). Thus, ∂1αh–1 is bounded over L2((Ω×Ωσ)). It follows that B ∈ K(L2((Ω×Ωσ))) because it is a product of the compact operator h–1 and the bounded operator M∂1αh–1.

From the compactness of B, it follows thatby the weak-* convergence of y to y. We conclude that jp ⇀ jp′ and hence that (ρ, jp) = (ρ′, jp′), completing the proof.

Corollary 1.is lower-semicontinuous and also weakly lower-semicontinuous.

Proof. Let . From Theorem 1, we then obtainwhere Γ → (ρ, jp). Hence, F is weakly lower-semicontinuous. By Mazur’s lemma,[15] weak lower semicontinuity of a convex function implies strong lower semicontinuity.

Corollary 2. IfF(ρ, jp) < +∞, then the infimum in the CDFT constrained-search functional is a minimum:

Proof. One simply takes (ρ, jp) = (ρ, jp) for all n and applies Theorem 1.

In conclusion, we have extended Theorem 4.4 of ref (1) to CDFT. As immediate corollaries, the constrained-search functional F(ρ, jp) is lower-semicontinuous and expectation-valued, that is, if F(ρ, jp) < +∞, then there exists Γ → (ρ, jp) such that F(ρ, jp) = Tr(H0Γ). These mathematical results are the final pieces in the puzzle of placing CDFT on a solid mathematical ground in a similar manner as done by Lieb for standard DFT.