Uniform magnetic fields in density-functional theory

We construct a density-functional formalism adapted to uniform external magnetic fields that is intermediate between conventional Density Functional Theory and Current-Density Functional Theory (CDFT). In the intermediate theory, which we term LDFT, the basic variables are the density, the canonical momentum, and the paramagnetic contribution to the magnetic moment. Both a constrained-search formulation and a convex formulation in terms of Legendre--Fenchel transformations are constructed. Many theoretical issues in CDFT find simplified analogues in LDFT. We prove results concerning $N$-representability, Hohenberg--Kohn-like mappings, existence of minimizers in the constrained-search expression, and a restricted analogue to gauge invariance. The issue of additivity of the energy over non-interacting subsystems, which is qualitatively different in LDFT and CDFT, is also discussed.

We construct a density-functional formalism adapted to uniform external magnetic fields that is intermediate between conventional Density Functional Theory and Current-Density Functional Theory (CDFT). In the intermediate theory, which we term LDFT, the basic variables are the density, the canonical momentum, and the paramagnetic contribution to the magnetic moment. Both a constrained-search formulation and a convex formulation in terms of Legendre-Fenchel transformations are constructed. Many theoretical issues in CDFT find simplified analogues in LDFT. We prove results concerning N -representability, Hohenberg-Kohn-like mappings, existence of minimizers in the constrained-search expression, and a restricted analogue to gauge invariance. The issue of additivity of the energy over non-interacting subsystems, which is qualitatively different in LDFT and CDFT, is also discussed.

I. INTRODUCTION
Ground-state density functional theory (DFT) is a widely successful electronic structure method due to its good trade off between accuracy and computational cost. Although DFT is commonly used as a pragmatic approach to compute the response of a molecule to an external magnetic field, this application is formally outside the scope of the theory. None of the independent theoretical foundations of DFT-that is, the Hohenberg-Kohn (HK) theorem [1], the constrained-search formulation [2], and Lieb's convex analysis formulation [3,4]-directly allows for a magnetic vector potential in the Hamiltonian. All three of these frameworks establish, firstly, that the density-scalar potential interaction can be isolated from other energy terms. Secondly, the remaining energy terms can be obtained from a universal density functional, which requires the density but not the external potential as input. Some generalizations of DFT that allow for a magnetic vector potential in the Hamiltonian have been proposed. Grayce and Harris dropped the universality property and formulated magnetic-field DFT, or BDFT, which can be viewed as a family of densityfunctional theories parametrized by the external magnetic field [5,6]. For a vanishing magnetic field, the parametrized theory coincides with conventional DFT. The resulting density functional is semi-universal-that is, universal with respect to the scalar potential, but not with respect to the magnetic field. The dominant alternative to BDFT that preserves universality is the paramagnetic current-density functional theory (CDFT) due to Vignale and Rasolt [7,8].
Here a universal functional of the density and paramagnetic current density is established and many results from conventional DFT carry over to this framework. Both a constrained-search formulation and a generalization of Lieb's convex formulation are possible for paramagnetic CDFT. The HK theorem of conventional DFT has only a weaker analogue in CDFT, which is sufficient for defining a ground-state energy CDFT functional but not for generalization of all the stronger statements that can be made based on the HK mapping.
In the present work, we explore a framework intermediate between DFT and CDFT. By restricting attention to uniform magnetic fields, we obtain a theory where universality does not require the full paramagnetic current density as an additional variable besides the density. Instead, the canonical momentum and the paramagnetic moment (equivalent to canonical angular momentum in the absence of spin) are sufficient. The resulting theory is termed LDFT. The reduction from an infinite-dimensional vector field to two three-dimensional quantities represents a substantial simplification. In some respects it is also a qualitative change, since the current density is a local quantity, whereas linear momentum and magnetic moment are global quantities. Both the orbital and spin magnetic moments contribute to the total magnetic moment. LDFT can therefore also be viewed as a minimal framework for incorporating spin dependence in a universal density functional.
The outline of this paper is as follows. In Sec. II, we review paramagnetic CDFT and adapt some of the technical mathematical details to allow for uniform magnetic fields. We continue in Sec. III by detailing the formulation of LDFT. Both constrained-search and convex formulations are given. Some HK-like results, which turn out to be stronger than the CDFT analogues, are also explored along with issues such as expectation-valuedness and Nrepresentability. In Sec. IV we discuss Kohn-Sham theory, with focus on invariance with respect to gauge degrees of freedom and additive separability of the exchange-correlation functional. Next, in Sec. V, we briefly comment on the possibility of using gauge-invariant, physical quantities instead of the paramagnetic current density or canonical momenta. We also comment on a recent hybrid formulation due to Pan and Sahni featuring both the canonical, gauge-dependent magnetic moment and the physical current density [9]. Finally, in Sec. VI, we give some concluding remarks.

II. REVIEW OF PARAMAGNETIC CURRENT-DENSITY FUNCTIONAL THEORY
Current-density functional theory is the natural generalization of DFT to the case when there is an external magnetic vector potential. For an N -electron system, we consider Schrödinger-Pauli Hamiltonians of the type where atomic units are used,Ŝ k = 1 2σ k is a spin operator, B(r k ) := ∇ k × A(r k ) is the magnetic field, and g s is the electron spin g-factor. Although g s = 2 is the empirically relevant value for a non-relativistic theory, the CDFT formalism is mathematically valid for any value of g s . When we are only interested in analysing orbital effects and states that are ground states in the absence of the spin-Zeeman term, we are therefore free to set g s = 0. In what follows, we consider this parameter to be fixed but arbitrary.
We introduce the short-hand notation for the pairing integral of two scalar fields f and g, and similarly for vector fields (where the pointwise product inside the integral is replaced by the scalar product). Our point of departure is the constrained-search expression introduced by Vignale and Rasolt. For given external potentials (v, A), we can write the ground-state energy as where the standard notation ψ → (ρ, j m ) means that the wave function gives rise to these densities-that is, ρ ψ = ρ and j m;ψ = j m . In detail, the densities corresponding to a wave function ψ (and similarly for a mixed state Γ) are and j m;ψ (r 1 ) = N Im Above, we let x k = (r k , ω k ) and the integral sign denotes integration both over spatial coordinates r k and over discrete spin coordinates ω k ∈ S. The magnetization current density j m (r) = j p (r) + g s ∇ × m(r) is the sum of the paramagnetic current density and the spin-current density. Both the magnetization current density and the paramagnetic current density are gauge-dependent quantities. The gauge-invariant, physical current density is given by j(r) = j m (r) + ρ(r) A(r). To obtain the paramagnetic term (j m |A) above, it is necessary to transfer the curl operator from the magnetic vector potential to the spin density-that is, (m|∇ × A) = (∇ × m|A). As is common, we understand the kinetic energy operators 1 2 (−i∇ l + A(r l )) 2 and − 1 2 ∇ 2 l as quadratic forms. When the spin-Zeeman term is understood as arising from the Pauli kinetic energy 1 2 (σ l · (−i∇ l + A(r l ))) 2 , this same point of view means that B(r l ) = ∇ l × A(r l ) is a distributional derivative, to be transferred to the wave function product it is integrated with. Hence, (m|∇ × A) = (∇ × m|A) by definition.
Two remarks can made at this stage about Eq. (3): Firstly, the search domains for the different minimizations should be specified. Most of the appeal of density functional theory comes from the fact that the inner minimization is a universal functional of (ρ, j m ). Hence, not only the expectation value itself in the inner minimization, but also the inner search domain should be free of dependencies on the external potentials. To setup a Lieb formalism, also the outer search domain needs to be independent of the external potentials. Secondly, and even more fundamentally for a density functional theory, the expectation value of the kinetic energy has been split into a diamagnetic, a paramagnetic, and a canonical part that may not be separately finite. This is necessary to express E(v, A) in terms of a current-density functional at all.

A. Choice of function spaces
For the ground-state problem of the Schrödinger-Pauli Hamiltonian, the natural choice of wave-function space is the magnetic Sobolev space defined by L 2 functions with finite physical kinetic energy, where we leave implicit the restriction to properly anti-symmetric and normalized ψ (strictly speaking, H 1 A (R 3N ×S N ) is then a subset of a magnetic Sobolev space) and the index l should be understood as a generic particle index-that is, the stated condition holds for all 1 ≤ l ≤ N . The magnetic Sobolev space H 1 A is natural since it makes the weak formulation of the ground-state problem well defined. In particular, any eigenfunction is in H 1 A . On the other hand, this space depends on the magnetic vector potential. Note also that decomposition of the kinetic energy into terms such as ψ| − 1 2 l ∇ 2 l |ψ , (j p;ψ |A), and (ρ ψ |A 2 ) is not in general possible since finite physical kinetic energy does not guarantee that these terms are separately finite. For the purposes of constructing a CDFT, these two facts are serious obstacles. To ensure that all three terms are separately finite, we may make the additional assumption that the paramagnetic term, (j p;ψ |A), is finite. This suffices because the other two terms are always positive. Alternatively, it suffices to assume that the sum of the canonical kinetic energy and the diamagnetic term is finite as this in turn entails a finite paramagnetic term.
is a condition on the physical kinetic energy, excluding the spin-Zeeman term. To ensure the finiteness of (j m |A), we must verify that (∇ × m|A), in addition to (j p |A), is finite. Fortunately, this follows automatically. The Cauchy-Schwarz inequality directly yields since the squared spin operator for a single particle index is a multiple of the identity operator,Ŝ 2 1 = 3 4 . For the spin-dependent term, we now have the following implication: ψ ∈ H 1 A (R 3N × S N ) and finiteness of (j p |A) entails finiteness of the canonical kinetic energy and the diamagnetic term, which in turn entails finiteness of (∇ × m ψ |A) and (j m |A).
A wave function with finite canonical kinetic energy belongs to the standard Sobolev space, If the magnetic vector potentials are confined to the space of bounded functions, A ∈ L ∞ , no further restrictions beyond Eq. (9) are needed [10]. However, in order to allow for the present CDFT formalism to later be specialized to LDFT in Sec. III, we shall choose a wave-function space that automatically gives rise to a finite physical kinetic energy and a finite paramagnetic term when the external vector potential grows linearly without bound such as A(r) = 1 2 B × r. As a by-product, this also ensures a well-defined magnetic moment. To achieve this, we introduce a weight function g(r) and take the wave-function space to be defined by We have in mind the specific choice which guarantees a finite paramagnetic term (j m |A) and finite magnetic moments. However, much of the discussion remains valid for a generic g(r) bounded away from zero to ensure that 1/g(r) ∈ L ∞ (R 3 ) is a bounded function. This slightly generalizes the previous work by Laestadius [10], where the choice g = 1 was implicit-that is, the unweighted Sobolev space was used. The membership g j p ∈ L 1 is a straightforward application of the Cauchy-Schwarz inequality and the fact that g 2 ρ and the kinetic energy density are L 1 functions: A similar Cauchy-Schwarz estimate shows that also the spin-current belongs to the same space, g∇ × m ∈ L 1 . Hence, the total magnetization current j m is in the same space.
To summarize, Eq. (10) entails that the densities belong to the function spaces This is a slight modification of the result ρ ∈ L 1 ∩ L 3 found by Lieb [3] for a standard, unweighted Sobolev space. The Hölder inequality guarantees finite interactions (ρ|v), (ρ|A 2 ), and (j m |A) when the external potentials belong to the spaces Indeed, Y * is the continuous dual of the weighted Lebesgue space Y , whereas X * is the continuous dual space of X, a Banach spaces under the norm ρ X := g 2 ρ L 1 + ρ L 3 . With this choice, the density-dependent energy term (ρ|v + 1 2 A 2 ) is guaranteed to be finite for all ρ ∈ X, v ∈ X * , and A ∈ Y * . As a consequence, the diamagnetic potential must also belong to the dual space of the densities-that is, A 2 ∈ X * and v + 1 2 A 2 ∈ X * .
We finally remark that finiteness of the canonical kinetic energy, Eq. (10), implies finiteness of the j p -corrected von Weizsäcker energy [10,11]: Recent work has shown this bound to be valid also when j p;ψ is replaced by j m;ψ [12].

B. Constrained search and Lieb formulations of CDFT
Having discussed the search domains, we introduce the notationĤ 0 =Ĥ(0, 0) and define the universal functionals Within the present, mildly restricted choice of search domains, the first functional is simply the pure-state constrainedsearch functional of Vignale and Rasolt [7], while the second functional is its extension to mixed states, with D denoting the convex set of all valid density operators [13]. Because all pure states ψ can be represented by density operators Γ = |ψ ψ|, it follows immediately that F VR,DM (ρ, j m ) ≤ F VR (ρ, j m ). The pure-and mixed-state functionals can be used interchangeably in the CDFT variation principle, The last equality follows because E(v, A) could just as well have been expressed as an infimum over mixed states. Because our choice of wave-function space amounts to a mild restriction to wave functions with well-defined secondorder moments, there is a risk that some potentials (v ′ , A ′ ) ∈ X * × Y * have ground states that decay so slowly that ψ ′ / ∈ H 1 0 and (ρ ′ , j ′ m ) / ∈ X × Y . Such ground states could have an infinite diamagnetic energy, which would be compensated by an infinite, negative paramagnetic energy. For the present purposes, we simply accept the fact that minima attained by slowly decaying ground states without well-defined second-order moments may be lost in this formulation. However, the restrictions on the function space do not affect the value of the infimum defining E(v, A), it only means that there are more instances where the infimum is not a minimum, attained by a wave function-for a careful discussion of this aspects, we refer to Kvaal et al. [14]. We expect molecular systems to feature exponentially decaying ground states, so this mild restriction leaves many areas of application unaffected. Several rigorous results exist that establish, for example, the exponential decay of ground-state wave functions for wide classes of potentials in the absence of magnetic fields [15] and in the presence of constant or asymptotically constant fields [16,17].
We now turn to the convexity properties of the functionals. The pure-state constrained.search functional F VR (ρ, j m ) is not convex-see Proposition 8 in Ref. 10, where we here in addition assume that the ground states in the counterexample are elements of H 1 0 . By contrast, the mixed-state functional F VR,DM (ρ, j m ) is jointly convex in (ρ, j m ) because the map Γ → (ρ, j m ) is linear. If the ground-state energy functional E(v, A) had featured a universal functional and potentials paired linearly with densities, it would have had the form of a Legendre-Fenchel transform and be manifestly concave. However, the diamagnetic term 1 2 (ρ|A 2 ) prevents such an identification and the existence of diamagnetic molecules indeed shows that E(v, A) is not concave in A. The non-concavity can be remedied by a change of variables that absorbs the diamagnetic term into the scalar potential [13], Since u ∈ X * follows from v ∈ X * and A ∈ Y * , this change of variables "stays within" the already specified space of scalar potentials. DefiningH(u, A) := H(u − 1 2 A 2 , A), we may now writē Note thatĒ(u, A) = E(u − 1 2 A 2 , A). It can be remarked that not all potential pairs in X * × Y * correspond to physical systems. For example, some pairs correspond to harmonic oscillator-type potentials with negative sign. A concrete example is a one-electron system subject to u(r) = −1/r and A(r) = 1 2 B × r; in this case the HamiltonianH(u, A) fails to be bounded from below as there is no diamagnetic term to balance the orbital Zeeman effect [18].
The functionalĒ(u, A) has the form of a Legendre-Fenchel transform and is jointly concave in (u, A). Repeated Legendre-Fenchel transformation now yields the conjugate pair Any function given as a Legendre-Fenchel transform automatically has some regularity. Thus,F is convex and lower semicontinuous (l.s.c.), whileĒ is concave and upper semicontinuous (u.s.c.); see, for example, Ref. [19]. Lower semicontinuity of f (x) means that, if x k → x, then lim inf k→∞ f (x k ) ≥ f (x) (and similarly for u.s.c. functions). All three universal functionals F VR , F VR,DM , andF are equivalent in the sense that they can be used interchangeably in the CDFT variation principle in Eqs. (20), (22), and (24), producing the right ground-state energy. In the language of Ref. [20] they are admissible functionals. However,F is the unique admissible functional that is both convex and l.s.c.-that is, it is the unique functional that satisfies the variation principle (23). This is significant, because this functional it the only one that can-in principle, at least-be computed fromĒ(u, A).
Removing the minimization in Eq. (22) and moving the density-potential pairing integrals to the left-hand side, we obtainĒ for all (ρ, j m ) ∈ X × Y and all (u, A) ∈ X * × Y * . Taking the supremum over potentials, The above argument is not specific to F VR,DM , but valid for any admissible functional. Therefore, in general,F ≤ F for any admissible F , such as F VR . For the l.s.c. ofF , in particular, see the straightforward argument adapted to the CDFT setting in the proof of Proposition 12 in Ref. 10.
The non-convexity of the pure-state functional means that, in general,F (ρ, j m ) = F VR (ρ, j m ). The identification of F and F VR,DM hinges on the l.s.c. of the latter functional, which is presently an open question; see Ref. 14. If F VR,DM is l.s.c., it follows that F VR,DM =F sinceF is unique. It is remarkable that, in standard DFT, we have the result F (ρ) = F DM (ρ) [3].

III. LDFT: A DFT FORMALISM FOR AFFINE VECTOR POTENTIALS
We study in this section how the paramagnetic CDFT of Vignale and Rasolt can be specialized and simplified when the vector-potential space is suitably constrained. Let us restrict attention to magnetic vector potentials that can be Taylor expanded around some reference point G, Truncating at some fixed order n is equivalent to restricting the function space of allowed vector potentials to polynomials of degree n in (x, y, z). In what follows, we study the case n = 1 with the additional gauge condition ∇ · A(r) = 0. However, everything can be straightforwardly generalized to an arbitrary finite order n without additional gauge constraints.
We thus focus here on vector potentials of the form A = 1 2 B × r. We find it instructive to not completely eliminate all gauge degrees of freedom, allowing constant shifts a of the vector potentials. Retention of these gauge degrees of freedom makes it easier to discuss, for example, additivity of physical energies over independent subsystems. Hence, we take the scalar potential space to be X * (as in the previous section) and the vector potential space to be a six-dimensional vector space, where the weight function is g(r) = 1 + |r| 2 . This choice of function space guarantees that both the canonical momentum p = j p dr and the paramagnetic moment L C = (r − C) × j m dr are finite. It is also worth being explicit about the role of spin: Hence, if g s = 0, then the spin degrees of freedom contribute to the magnetic moment but not to the linear momentum. The physical momentum can similarly be expressed as π = jdr = p + ρAdr and the physical magnetic moment as J C = (r − C) × jdr = L C + (r − C) × ρAdr. Magnetic moments relative to different reference points are related as Moreover, for any other eigenstate of the Hamiltonian (and, in particular, the ground state), we have where µ G = (r − G) ρdr is the electric dipole moment relative to G. As a consequence, for energy eigenstates, the physical magnetic moment J C is independent of the reference point C.
The only allowed gauge transformations are now gauge shifts-that is, constant shifts of the magnetic vector potential. The canonical momentum and paramagnetic moment transform in a simple manner under such gauge shifts, where N = ρdr is the number of electrons. We also define a gauge-shift invariant quantity that we term the intrinsic magnetic moment, The intrinsic magnetic moment is identical to the paramagnetic moment with respect to the average position of the electrons or centre of charge, Λ = L R with R = G + N −1 µ G . In some sense, the intrinsic magnetic moment extracts the part of L G that is independent of the physically arbitrary reference point G. It is therefore better suited as a parameter in a functional describing a physical energy, which must be independent of arbitrary reference points. We return to this point in Sec. IV, where we discuss the Kohn-Sham decomposition of the total energy.

A. Constrained-search and Lieb formulations of LDFT
Let us now define a ground-state functional on the space X * × Y * lin by where the constrained-search functional can be defined either with pure states or with mixed states, For a magnetization current density consistent with the canonical momentum and magnetic moment-that is, p = j m dr and L G = r G × j m dr-we have the inequalities Next, noting that A 2 ∈ X * , we may introduce a new potential u := v + 1 2 A 2 ∈ X * that remains in the same function space as v ∈ X * . After this reparametrization, we obtain This expression takes the form of a Legendre-Fenchel transformation (concave conjugation). Hence, a Lieb formalism is obtained for affine vector potentials. Repeated conjugation gives the conjugate pair These two transformations constitute a generalization of Lieb's formulation of DFT to include uniform magnetic fields.

B. Expectation-valuedness of the pure state functional
In the terminology of Kvaal and Helgaker, a DFT functional is said to be expectation-valued if its finite function values are always attained by expectation value ψ|Ĥ 0 |ψ for some state ψ [20]. The conventional DFT functionals are expectation valued [3] as is the pure-state CDFT functional F VR [10]. It is an open question whether or not the universal CDFT functional F VR,DM is expectation valued. Below we modify the CDFT argument in Ref. 10 to prove that also the pure state LDFT functional F has this property.
We first introduce the set of density triples corresponding to N electrons, All densities that are pure-state N -representable (with finite kinetic energy) are contained in this set. The LDFT functional attains its minimum whenever it is finite. In the language of Ref. 20, F is expectation valued: Proof. Let {ψ k } ⊂ H 1 0 be a minimizing sequence such that and ψ k → (ρ, p, L G ) for each k. This means that, for any k, we have j m;k dr = p and r G × j m;k dr = L G , where j m;k := j m;ψ k . Furthermore, denoting the volume element of the combined coordinate and spin space which means that {ψ k } is bounded in H 1 0 . By the Banach-Alaoglu theorem ψ k ⇀ ψ weakly in H 1 0 for some ψ ∈ H 1 0 . The argument in the proof of Theorem 3.3 in Lieb's work [3] shows that ψ k → ψ in L 2 -norm and ψ → ρ. We may then conclude that ψ ∈ H 1 0 , because Furthermore, since ψ k ⇀ ψ in H 1 0 and the expectation value ofĤ 0 is weakly l.s.c., it holds that If we can prove j m;k → j m;ψ as well as and the proof would be complete. Equation (53) implies a simple rule for changing the reference point from G to another location in Eq. (54). Hence, without loss of generality, we simplify the notation by moving the reference point to the origin in what follows, writing r instead of r G . To prove convergence of the current-density sequence, we note that φ k := gψ k is a bounded sequence in L 2 . Again, by the Banach-Alaoglu theorem, φ k ⇀ φ weakly in L 2 for some φ ∈ L 2 . The argument used in the proof of Theorem 3.3 in Ref. 3 shows that (for a subsequence) φ k → φ in L 2 -norm. This argument relies on the fact that weak convergence φ k ⇀ φ is given and it only remains to show norm convergence φ k L 2 → φ L 2 to prove that φ k → φ in L 2 . Since lim k→∞ φ k L 2 ≥ φ L 2 is a generic consequence of weak convergence, only the reverse inequality needs to be proven below.
Let ε > 0 be given. Since g 2 ρ ∈ L 1 , we may choose a characteristic function χ of some bounded set in R 3 such that Next, let δ > 0 be given and choose R > 0 such that B c R |φ − gψ| 2 dτ ≤ δ 2 , where we leave the spin degrees of freedom implicit. Let M ⊂ B R be any measurable set and note that M φdτ = lim k M gψ k dτ = M gψdτ , implying that φ = gψ almost everywhere on B R . Choosing N such that k > N implies φ k − φ L 2 ≤ ε, we obtain We now show that, with l, l ′ ∈ {1, 2, 3} and a bounded self-adjoint operatorΩ with operator norm less than 1, as k → ∞. The difference between the indicated limit and the left hand-side vanishes because, using φ k − gψ L 2 → 0, we have and furthermore since r lΩ ψ * ∈ L 2 , Having established Eq. (59), we can setΩ to the identity operator to conclude that R 3 r l (j p;k ) l ′ dr → R 3 r l (j p;ψ ) l ′ dr. We can also setΩ to be any component of the spin operator to conclude that the spin-current contribution in Eq. (5) converges in the same way. Hence, we have An immediate generalization of the above theorem is obtained by noting that the proof remains valid when the electron-electron repulsion operator is scaled by a non-negative parameter λ. The only parts affected by the scaling are Eqs. (48) and (52), which still hold for any non-negative scaling.

Theorem 2. Consider the modified potential-free Hamiltonian
For all λ ≥ 0 and all triples (ρ, p, L G ) ∈ J N , there exists a ψ λ ∈ H 1 0 such that ψ λ → (ρ, p, L G ) and This form of the theorem is relevant for the adiabatic-connection expression for the exchange-correlation energy [21][22][23][24]. The adiabatic-connection expression is essentially the integral where the last equality relies on the Hellmann-Feynman theorem and the adiabaticity assumption that F λ is differentiable with respect to λ-in fact, it is sufficient that the left-or right-derivative with respect to λ exists. This form, combined with approximate Lieb optimization to determine the scalar potential needed to enforce the density constraint, has been used to characterize the exchange-correlation functional in standard DFT [25][26][27][28][29][30][31][32]. Furthermore, we introduce the density-matrix functional Although the above proof does not apply to this functional so that its expectation-valuedness is presently not established, it will be useful in later discussions-see Sec. IV.

C. Mixed-state N -representability
In this section, we focus on states with zero spin S = 0, which is sufficient to establish a general N -representability result for LDFT. Mixed-state N -representability has previously been established very generally for pairs (ρ, j p ) [13].
The N -representability problem for (ρ, p, L G ) is therefore solved if one can always construct a current density j p that is compatible with the given triple. This can be achieved, for example, by minimization of the current-correction to the von Weizäcker energy, |j p | 2 /2ρdr, subject the constraints p = j p dr and Λ = r R × j p dr. The minimizer is a velocity field of the form where the vector field κ (in the notation in Ref. [13]) is twice the paramagnetic velocity, and the constants (Lagrange multipliers) ζ and ν are chosen so as to satisfy the constraints. The latter multiplier is in fact the vorticity of this particular velocity field. The equations to determine the multipliers are straightforward, where the last equality defines the moment-of-inertia tensor, Q. Inserting our specific κ into the general density-matrix construction of Sec. IV.A of Ref. [13], we obtain where λ, µ ≥ 2p π ( 1 4 N ρ q ) 2p/3 and p, q > 1 are Hölder conjugates. Because ρ ∈ L 1 ∩ L 3 , we may for simplicity fix, for example, p = q = 2. Here, we have modified the original form slightly by replacing r and s by r R and s R , respectively, where it makes a difference in the exponentials; note that r − s = r R − s R .
The constant shift of the coordinates does not affect the results in Ref. [13], to which we refer for proofs that D λµ is a valid one-particle reduced density matrix (e.g., it is shown that spatial occupations numbers in the interval [0, 2]). The density matrix D λµ reproduces (ρ, j p ) and therefore also p, Λ, and L G . Given our choice of function spaces, this construction thus yields a strong N -representability result for LDFT: The canonical kinetic energy of this explicit construction is given by where T W = |∇ρ| 2 /8ρdr is the von Weizsäcker kinetic energy, setting an upper limit on the density-matrix functional F λ=0 DM (ρ, p, L G ).

D. Hohenberg-Kohn-like results
Paramagnetic CDFT admits a HK-like result in the form that (v, A)-representable ground-state densities (ρ, j m ) determine the ground-state wave function, provided that it is non-degenerate [7,33]. In the degenerate case, we can consider two HamiltoniansH(u k , A k ), k = 1, 2, with ground-states Γ k . If the ground states share the same densities (ρ 1 , j m;1 ) = (ρ 2 , j m;2 ), then Γ 1 is also a ground state ofH(u 2 , A 2 ) and vice versa (see Sec. III.B in Ref. 34). The literature proofs are given for j p (or g s = 0) but hold also for j m with trivial changes. In fact, HK-like results are possible whenever potentials are paired linearly with densities in the ground-state energy expression. We term this type of results weak HK-like results, as they are weaker than the HK theorem in standard DFT, where also the potentials are determined by the ground-state densities.
Below, we state three weak HK-like results specialized to the LDFT setting, all equally valid with (g s = 2) and without (g s = 0) inclusion of the spin-Zeeman term in the Hamiltonian. Theorem 4. Let two or more potential triples (u k , a k , B k ), with k = 1, . . . , K ≥ 2, be given. Let Γ k denote a pure or mixed ground state ofH(u k , a k , B k ). Suppose that all Γ k give rise to the same density triples (ρ Γ k , p Γ k , L G;Γ k ) = (ρ, p, L G ). Then Γ k is also a ground state ofH(u l , a l , B l ), for all k, l.
Proof. The energy that the lth Hamiltonian assigns to the kth state can be written as Tr(Γ kH (u l , a l , B l )) = (u l |ρ) + a l · p + 1 2 B l · L G + Tr(Γ kĤ0 ).
Adding the same energy expression with the indices interchanged, we obtain Tr(Γ kH (u l , a l , B l )) + Tr(Γ lH (u k , a k , B k )) = Tr(Γ kH (u k , a k , B k )) + Tr(Γ lH (u l , a l , B l )).
Together with the variation principle Tr(Γ kH (u l , a l , B l )) ≥ Tr(Γ lH (u l , a l , B l )), the result now follows.
The above result can be generalized slightly by taking convex combinations of density matrices and potentials. The set of potential triples {(u, a, B)} that give rise to a given density triple is in fact a convex set. Moreover, each such potential triple is a subgradient of the functionalF at (ρ, p, L G ).
Under the additional assumption that the ground state is non-degenerate, we obtain a HK-like result in the more familiar form that ground-state densities (ρ, p, L G ) determine the ground-state wave function.
Theorem 5. Let two or more potential triples (u k , a k , B k ) with k = 1, . . . , K ≥ 2 be given. Suppose that each Hamilto-nianH(u k , a k , B k ) has a non-degenerate ground state ψ k that gives rise to the same density triple (ρ ψ k , p ψ k , L G;ψ k ) = (ρ, p, L G ). Then ψ k and ψ l are equal up to a global phase for all k, l.
A corollary is that non-degenerate ground-state densities (ρ, p, L G ) also determine the paramagnetic current density j p as well as the magnetization current density j m . We now recall that the eigenstates of the Hamiltonian H(v, A) have divergence-free physical current densities, ∇ · j = ∇ · (j p + ρA) = 0 and a vanishing physical momentum. If two or more potential triples (u k , a k , B k ) give rise to the same non-degenerate ground-state density triple (ρ, p, L G ), and consequently the same ground-state wave function ψ and paramagnetic current density j p , then Hence, for all k, l, we have The first term can be eliminated by exploiting the fact that p k = p l = p and the relation Inserting where R is the center of charge. If B k = B l , it follows that the ground-state density ρ must be cylindrically symmetric about an axis passing through R and directed along B k − B l . Hence, for non-cylindrically symmetric densities ρ, LDFT does admit a strong HK-like result in the form that the external potentials are determined: Theorem 6. Let two or more potential triples (u k , a k , B k ) with k = 1, . . . , K ≥ 2 be given. Suppose that each HamiltonianH(u k , a k , B k ) has a non-degenerate ground state ψ k with the same density (ρ ψ k , p ψ k , L G;ψ k ) = (ρ, p, L G ). Suppose further that ρ is not cylindrically symmetric about any axis. Then ψ k and ψ l are equal up to a global phase, u k = u l + const and (a k , B k ) = (a l , B l ), for all k, l.
To establish the HK-like result that (ρ, p, L G ) determines the external potentials, we exploited the fact that this triple determines j p (and the ground state ψ) under a non-degeneracy assumption. Intuitively, we expect (ρ, Λ) in a similar fashion to determine the gauge-shift invariant components of the potentials-that is, B and u to within a constant, leaving a undetermined. To investigate this, we must deal with the fact that a gauge shift in the magnetic vector potential results in a non-trivial change of the phase of the ground-state wave function, thereby affecting the paramagnetic current density. Theorem 7. Let two potential triples (u k , a k , B k ) with k = 1, 2 be given. Suppose that each HamiltonianH(u k , a k , B k ) has a ground state ψ k that gives rise to the same density pair (ρ ψ k , Λ ψ k ) = (ρ, Λ). Then, for a constant c (determined below), the state ψ ′ 1 = ψ 2 e −ic· j rj is a ground state ofH(u 1 , a 1 , B 1 ) and ψ ′ 2 = ψ 1 e ic· j rj is a ground state of H(u 2 , a 2 , B 2 ).
Proof. Inserting the identity L k;G = Λ + (R − G) × p k into the orbital Zeeman term, we obtain Noting that A k (R) = a k + 1 2 B k × (R − G) and adding and subtracting the terms needed to form E 2 := ψ 2 |H(u 2 , a 2 , B 2 ), we next find As in the standard HK proof, we obtain a new equation by interchanging the indices 1 and 2, being careful to also replace occurrences of c by −c, yielding where E 1 := ψ 1 |H(u 1 , a 1 , B 1 )|ψ 1 . Adding the last two equations and setting c = A 1 (R) − A 2 (R), we obtain By the variation principle, we also have E ′ 1 ≥ E 1 and E ′ 2 ≥ E 2 . The only possibility is E ′ 1 = E 1 and E ′ 2 = E 2 . Thus, both ψ 1 and ψ ′ 1 must be ground states ofH(u 1 , a 1 , B 1 ). Likewise, both ψ 2 and ψ ′ 2 must be ground states of H(u 2 , a 2 , B 2 ).
Under the additional non-degeneracy assumption, it follows that two HamiltoniansH(u k , a k , B k ) with the same ground-state density pair (ρ, Λ) have ground states related by ψ 2 = ψ 1 e ic· j rj with c = A 1 (R) − A 2 (R). Consequently, the paramagnetic current densities are related as j p;2 = j p;1 + ρc.
The corresponding physical current densities must be divergence free, Subtracting the first of these equations from the second, we obtain If B 1 = B 2 , it follows that ρ is cylindrically symmetric about an axis passing through the center of charge in the direction B 2 − B 1 . Hence, only components of the magnetic field parallel to a cylindrical symmetry axis are undetermined. To summarize, we have proved the following theorem: Theorem 8. Let two or more potential triples (u k , a k , B k ) with k = 1, . . . , K ≥ 2 be given. Suppose that each HamiltonianH(u k , a k , B k ) has a non-degenerate ground state ψ k with the same density triple (ρ ψ k , Λ ψ k ) = (ρ, Λ). Suppose further that ρ is not cylindrically symmetric about any axis. Then ψ k = ψ l e i(A k (R)−A l (R))· j rj (to within a global phase) and B k = B l , for all k, l.
A brief remark can be made about the case when there is a rotational symmetry axis. From symmetry considerations, we conclude that projection of the rotation axes from B k and Λ results in (anti-)parallel vectors. Let C denote the projection onto the space spanned by all rotational symmetry axes and C ⊥ the projection onto the complement. Then Another triple that determines the external magnetic field is the gauge-shift invariant (ρ, Λ, J R ), where J R is the physical magnetic momentum with respect to the center of charge R, The last term vanishes because µ R = (r − R)ρdr = 0 is the dipole moment with respect to the center of charge.
Rewriting the other integral and noting that L R = Λ, we obtain Recalling the definition of the moment-of-inertia tensor, Q = (|r R | 2 I − r R r T R )ρdr, it follows that The inertia tensor Q is symmetric and it is therefore always possible to rotate the coordinate axes so that Q becomes diagonal. Without loss of generality, we therefore write This matrix is always positive definite and invertible. Hence, we find To determine the magnetic field, it was not even necessary to assume that (ρ, Λ, J R ) comes from a ground state. As long as (ρ, Λ) has been obtained from some state and J R has been calculated from the same state and some additional vector potential A = a + 1 2 B × (r − G), the above expression for B is valid.

E. Constructive relationship between B and ground-state densities (ρ, jm)
The above HK-like results have been established using proofs by contradiction. Going beyond the strict LDFT framework and studying the current density, we may derive a more explicit, constructive relationship. In particular, for the density pair (ρ, j m ), arising from a ground state of some potential (u, A = a + 1 2 B × r G ), a constructive relationship can be established.
First note that the condition of a vanishing physical momentum can be expressed as π = p + N A(R) + 1 2 B × µ R . Because the dipole moment relative to the center of charge vanishes, we find that the canonical momentum determines the vector potential at the center of charge, p = −A(R)/N . The condition that the physical current is divergence-free can then be expressed as where (r R × ∇ρ) T is a row vector. Multiplying by the column vector ρ −1 r R × ∇ρ and integrating over all space, we obtain where the first equality defines M. The rationale for the factor 1/ρ is simply that it makes each element of M bounded by a weighted von Weizäcker energy. Therefore, M is guaranteed to be finite, given our choice of function spaces. As long as the resulting integral is well defined, the factor 1/ρ can be replaced by any function of ρ. An alternative expression for M, which is somewhat more conducive to analysis, is the form It can be verified (by writing out the expression in cylindrical coordinates) that M = 0 for a spherically symmetric density with rank(M) = 2 when there is a single cylindrical symmetry axis. For a non-cylindrically symmetric density, rank(M) = 3. Hence, for non-cylindrically symmetric densities ρ, where j m = j m + ρA(R) = j m − ρ p N is the deviation of j m from its spatially averaged value p = j m dr. An alternative constructive relationship is obtained by multiplying the condition ∇ · j by a coordinate product (r a − R a )(r b − R b ), with a, b ∈ {x, y, z}, and integrating over all space: Noting that the term proportional to A c (R) vanishes by virtue of being a constant times components of the dipole moment relative to the center of charge, and moving the B d -dependent terms to the left-hand side, we obtain When a = b, the right-hand side resembles the two terms in the magnetic-moment component Λ c (with c = a and c = b) but with different signs: this expression is symmetric in a, b, whereas the magnetic moment is anti-symmetric. Contraction by the Levi-Civita tensor ǫ cab yields zero on both sides. We therefore contract by the absolute value of each element, |ǫ cab |, to obtain ρ|ǫ cab | ǫ ade r Rb r Re B d dr = −2 |ǫ cab |r Ra j m;b dr.
We write the right-hand side (excluding sign) more compactly using the notation, L + R;c = (r Rb j m;a + r Ra j m;b )dr = |ǫ cab | r Ra j m;b dr, where |ǫ cab | means the absolute value of the Levi-Civita tensor. The left-hand sides of Eqs. (98) and (99) define a linear transformation of the magnetic field. This transformation can therefore be represented by a matrix, Q + cd = ρ(r Rb ǫ ade r Re + r Ra ǫ bde r Re )dr = ρ|ǫ cab |r Rb ǫ ade r Re dr, or, more explicitly, The same terms appear in Q + and in the moment-of-inertia tensor Q; however, half the terms have different signs. Alignment of the coordinate to, for example, the principal axes makes the off-diagonal elements of Q + vanish. We then see that this matrix has full rank if, in such a coordinate system, ρx 2 R dr, ρy 2 R dr, and ρz 2 R dr are all distinct. This fails if the system has cylindrical symmetry. Unlike M above, Q + can also be rank-deficient in other cases, for example, when there is a 4-fold rotation symmetry axis rather than a continuous rotational symmetry. When the matrix has full rank, then

IV. KOHN-SHAM THEORY FOR LDFT
Kohn-Sham theory is formulated in a similar manner for DFT, CDFT, and LDFT, the fundamental distinction being what the Kohn-Sham system is required to reproduce. In CDFT, the Kohn-Sham system is a Slater determinant that reproduces (ρ, j m ). In LDFT, it is a Slater determinant which reproduces (ρ, p, L G ). More formally, we introduce the functionals where the minimization is over all Slater determinants Φ ∈ H 1 0 . The minimizers define the Kohn-Sham systems. A standard decomposition of the universal energy now yields where we drop the superscripts 'CDFT' and 'LDFT' and J(ρ) is the Hartree energy-that is, the classical self-repulsion of the charge distribution. A practical point related to the different constraints on the Kohn-Sham system is the calculation of magnetic properties. All properties that can be obtained as simple integrals over the current density are easy to compute in the CDFT Kohn-Sham framework. In the LDFT Kohn-Sham framework, however, only the linear momenta and magnetic moments are reproduced-the Kohn-Sham current density may differ from the current density of the interacting system. As a result, other magnetic properties are less accessible in LDFT than in CDFT.

A. Non-interacting N -and (v, A)-representability
The non-interacting N -representability problem is one of cases where the orbital and spin contributions to the magnetic moment need to be treated differently. Below, we focus on the case g s = 0 with only orbital contributions.
The Kohn-Sham theory of paramagnetic CDFT is known to suffer from N -representability problems in the case for N = 1. More precisely, a single Kohn-Sham orbital is only capable of reproducing densities (ρ, j p ) with vanishing paramagnetic vorticity ∇ × (j p /ρ) = 0 except possibly at isolated points where the curl is singular. The situation for N = 2 and N = 3 is not fully settled and for N ≥ 4 Lieb and Schrader have shown that all density pairs satisfying some regularity conditions are representable by a single Slater determinant [35].
When N = 1, LDFT shares some of the N -representability problems of paramagnetic CDFT. This follows because a one-particle wave function only has a non-vanishing canonical angular momentum when there are nodal lines [36,37]. A simple corollary is that (ρ, Λ) is not representable by a single orbital when ρ(r) > 0 and Λ = 0. Moreover, for all wave functions with a cylindrically symmetric density ρ(r) = |ψ(r)| 2 , the component of the intrinsic momentum Λ parallel to the rotation axis is quantized (see Theorem 2 in Ref. 37). The case N = 2 is much harder to analyse and it is an open question if there are any restriction on pairs (ρ, Λ) that can be represented by N = 2 and N = 3 orbitals. For N ≥ 4, Lieb and Schrader's Slater-determinant representability proof for (ρ, j p ) entails that all (ρ, Λ) are representable too. The representability problems in CDFT and LDFT are avoided if the Kohn-Sham orbitals are allowed to be fractionally occupied (see Ref. 13 for the CDFT case and Theorem 3 above for the LDFT case).
The issue of (v, A)-representability is less well characterized. However, two-dimensional counterexamples to (v, A)representability in CDFT have been constructed and these apply directly also to LDFT [38].
We remark that T s (ρ, p, L G ) and the functional F λ=0 (ρ, p, L G ) introduced in Theorem 2 coincide on the set of nondegenerate, non-interacting (v, A)-representable densities. This fact follows since the relevant Hamiltonian is in this case a pure one-electron operator, implying that the non-degenerate ground state is a Slater determinant. In general, however, they are related through T s (ρ, p, L G ) ≥ F λ=0 (ρ, p, L G ). The functional F λ=0 DM is the kinetic energy in an ensemble Kohn-Sham setting with fractional occupation numbers allowed. Clearly, T s (ρ, p, L G ) ≥ F λ=0 (ρ, p, L G ) ≥ F λ=0 DM (ρ, p, L G ) must hold, and with equalities on the set of non-degenerate, non-interacting (v, A)-representable densities.

B. Expectation-valuedness of the Kohn-Sham kinetic energy
From Theorem 6 in Ref. 39, the infimum in the definition of the CDFT functional T s in Eq. (104) is attained by some determinant. This result is premised on noninteracting N -representability, which presently has been proven generally in the case of four or more Kohn-Sham orbitals [35]. The same holds for the LDFT version of the functional defined in Eq. (105): Theorem 9. If N ≥ 4, then there exists a Slater determinant Φ ∈ H 1 0 such that Φ → ρ, p, L G and T s (ρ, p, L G ) = Φ|T |Φ . (108) Proof. Theorem 2 above with λ = 0 establishes a minimizer Φ ∈ H 1 0 with all specified properties, save for being a determinant. However, we can use the same argument as in the proof of Theorem 6 in Ref. 39 to conclude that Φ is a determinant of N orthonormal orbitals.
C. The exchange-correlation functional, gauge-shift invariance, and additive separability In Vignale and Rasolt's paramagnetic CDFT, there is a universal functional F VR (ρ, j m ), which depends on the gaugedependent paramagnetic current density. This gauge-dependence is necessary in their formalism since their universal functional includes the canonical kinetic energy. However, the exchange-correlation contribution to F VR (ρ, j m ) must be gauge invariant and should therefore only depend on the gauge-invariant "components" of (ρ, j m ). These components are the density itself and the paramagnetic vorticity ν := ∇ × (j m /ρ). In LDFT, a minimalistic version of the same problem appears: The universal functional depends on the full triple (ρ, p, L G ), but the exchange-correlation contribution should only depend on the "components" that are invariant under the restricted gauge transformation that are allowed-namely, constant shifts of the vector potential. By counting the degrees of freedom (6 in (p, L G ) and 3 in a), we conclude that these gauge-shift invariant components are contained in the pair (ρ, Λ). The fact that this pair is sufficient to uniquely determine B for a non-cylindrically symmetric ground-state density supports the conclusion that there are no additional gauge-shift invariant "components" of (ρ, p, L G ).
From these considerations, we may thus decompose the universal functional as where all the gauge-shift dependence is contained in the canonical kinetic energy term, T s (ρ, p, L G ), the arguments to F xc having been reduced to a gauge-shift independent density pair (compare Eq. (107)). Another type of constraint comes from requirement of additive separability-that is, the requirement that the energy is additive over well-separated non-interacting subsystems [40]. Additive separability is related to the stronger notions of size extensivity and size consistency [41][42][43]. This type of property is a desirable constraint on any practical approximate exchange-correlation functional. Consider a system consisting of S subsystems. Each subsystem is a ground state of some Hamiltonian H(v α , a, B) with density (ρ α , p α , L α;G ) (denoting subsystems by Greek indices). The total system is assumed to be the ground state of H(v, a, B) with v = α v α . Then additive separability amounts to the condition Let N α denote the number of electrons and R α the center of charge of subsystem α. We shall assume that the subsystems are well separated-that is, that the distances |R α − R β |, with α = β, are arbitrarily large compared with the extent of each subsystem density, so that we can write In other words, at any point r, only a single subsystem contributes to the total density. There is also a relation between the total intrinsic magnetic moment and those of the individual subsystems. To derive a detailed form of this relation, we begin by introducing notation for moments of inertia with respect to multiple reference points, The moments of inertia of the total system are defined analogously. Writing simply Q = Q R and Q α = Q α,Rα for moments of inertia with respect to the relevant center of charge, we have the relation Denoting the last term by G, we have Q = α Q α + G. In general, rank(G) = 2 when S = 1 or when the centres of charge of all subsystems are located on a line, and rank(G) = 3 otherwise. When G is invertible, with equality in the limit of infinite subsystem separations. Next, we turn to the physical magnetic moment, which by the ground-state assumption is independent of the reference point. The total physical magnetic moment can therefore be written in the form The physical magnetic moment is an extensive quantity and can thus also be written as a sum over subsystems, relying on the reference-point independence to move the reference point from R to R α . Inserting J α,Rα = Λ α + 1 2 Q α B, we find When G has full rank, we obtain with equality for infinite subsystem separation. This is another HK-like result: The total intrinsic magnetic moment and momentum-of-inertia tensor of a large total system, with three or more subsystems not located on a single line, determine the external magnetic field. This result can also be understood by introducing an effective angular velocity ω = Q −1 J R . Recall that the angular velocity vector and the angular momentum of a classical rigid body are related through J = Qω. Every point in a rigid body has the same angular velocity. When applied to a non-rigid body, ω can be interpreted as an effective, distance-weighted angular velocity, It is harder to sustain a significant effective angular velocity in a large system since it requires proportionally strong currents encircling the system. Consider, for example, a ring-shaped wire of radius R, wire cross section area A, constant electron density ρ within the wire, and current I. The angular velocity is then given by ω = I/(ρAR), which vanishes in the limit of infinite R and constant ρ, A, and I. When all subsystems are located on a single line, the matrix G is of rank 2, with zero eigenvalue in the direction vector of the line. Hence, only the components of B in the plane perpendicular to the line are uniquely determined: where G −1 is now the generalized inverse and equality holds in the limit of infinite subsystem separations. To reexpress the above relation in terms of Q instead of G, we begin by writing G = UγI 2 U T , where γ = α N α d 2 α is the two-fold degenerate non-zero eigenvalue and U is a 3 × 2 matrix containing the corresponding eigenvectors. The remaining eigenvector, with zero eigenvalue, is denotedd and is by stipulation parallel to all d α . Using the shorthand notation S = α Q α and the Woodbury matrix identity, we obtain Some further algebra yields with equality in the limit of infinite subsystem separations. Eq. (120) can now be expressed as where UU T = I −dd T is a projector onto the non-null space of G. Projection withdd T yields the remaining part of Q −1 Λ,dd which is also an intensive quantity. We also note a type of additivity for the projected intrinsic magnetic moment, Let us now consider local approximations to the exchange-correlation energy. An obvious first ansatz is the LDA-like form where only the magnitude |Λ| can enter without breaking rotational invariance. We assume that f (ρ(r), |Λ|) is negligible whenever ρ(r) is negligible. Exploiting the fact that the subsystems are separated and requiring additive separability, we now obtain This needs to hold for a large class of ground-state densities ρ α . In practice, we expect equality to hold also for the integrands, leading to We have seen above that the total |Λ| diverges to infinity in the limit of infinite subsystem separation. Consequently, the above relation can only be satisfied when f is independent of the second argument. A slightly more refined ansatz relies on the quantity β = −2Q −1 Λ instead of the intrinsic magnetic moment directly. Note that Q −1 is a non-local but tractable functional of the density. An LDA-like form is then given by Additive separability now leads to the condition where β = B or β = B ⊥ depending on whether G has full rank or not. Again, additive separability cannot be satisfied unless f is independent of its second argument. However, the prospects are better for approximate additive separability since β ≈ β α holds whenever the size or properties of the subsystems prevent large effective angular velocities ω α . A slightly more flexible LDA-like form is obtained by giving the functional access to data about the global dimensions of the system via the moment-of-inertia tensor, With this form, it becomes possible, for example, to project Λ onto the principal components of Q and to treat the components differently depending on the corresponding eigenvalue of Q.

V. THE PHYSICAL CURRENT DENSITY AND THE PHYSICAL MAGNETIC MOMENT
After the seminal work by Vignale and Rasolt, paramagnetic CDFT has become the dominant generalization of standard DFT to allow for arbitrary, external magnetic fields. It is a prima facie surprising feature that the paramagnetic current density appears as a basic variable in this formulation. A few works have attempted to construct an alternative formulation based on the physical current density [44][45][46][47][48]. However, none of the three foundations of DFT-that is, the HK theorem, the constrained-search approach, or Lieb's convex analysis framework-has been rigorously extended to physical current densities.
The existence of an adequate HK-like mapping from the physical density pair (ρ, j) to external potentials (v, A) (or an equivalence class of potentials that differ by a gauge transformation) is presently an open question [34,49]. The existence of a HK functional of the physical current density thus cannot be excluded. However, such a functional would require a different approach than in standard DFT since it would not have a straightforward variation principle [34,50,51].
A constrained-search approach suffers from the fact the inner minimization would be performed over a domain that depends on the external vector potential and the constrained-search functional would therefore not be universal, losing most of its formal and practical appeal. Consequently, we should also be careful in generalizing the common constrained-search notation 'ψ → ρ' or 'ψ → ρ, j m ' to a version with physical currents, 'ψ → ρ, j', as this leaves implicit the dependence of the wave-function domain on the vector potential.
A generalization of the Lieb formalism with the physical density as a variable would require the identification of a Legendre-Fenchel transformation pair of a concave ground-state energy functional E(v, A) and a convex universal functional F (ρ, j). A necessary condition for this is that external potentials are paired linearly with densities. However, E(v, A) is not concave in A and the interaction with the magnetic vector potential contains both linear and quadratic terms. An attempt to solve these problems using a change of variables w := v − 1 2 A 2 (in analogy with how u := v + 1 2 A 2 can be introduced in paramagnetic CDFT) leads to a direct contradiction (see Sec. IV.D in Ref. [34]). Counterexamples show that a change of variables u λ := v + λA 2 with λ < 1 2 A 2 cannot lead to a concave energy functional E(u λ , A). In summary, a satisfactory CDFT featuring the physical current density as a variable is neither presently available nor completely excluded.
Alongside the LDFT formalism presented above, we can imagine a formalism based instead on the physical magnetic moment. However, the difficulties for a physical CDFT would also be present for a physical LDFT. For example, the constrained-search domain indicated by 'ψ → ρ, J' would depend on the external potential, and therefore be nonuniversal. LDFT thus provides a simplified setting where such difficulties can be analysed and possibly circumvented. To date, no physical LDFT has been proposed or analysed in the literature. A hybrid formalism, based on the triple (ρ, L, j), has been proposed [9]. We now compare this approach with LDFT formalism.
A. Comparison of LDFT with a recent approach using (ρ, L, j) In Ref. 9 a HK-like result is introduced in a setting where only uniform magnetic fields and gauge-fixed vector potentials are allowed. After minor modifications, the result can be stated as follows: Theorem 10. Let two or more potential pairs (v k , B k ), with k = 1, . . . , K ≥ 2, be given. Let the vector potentials be gauge fixed to the form A k = 1 2 B k × r G . Suppose that ψ k is the non-degenerate ground state of H(v k , A k ) and that all ψ k give rise to the same density triple (ρ ψ k , L G;ψ k , j k ) = (ρ, L G , j). Then, for all k, l, it follows that (a) B k = B l , (b) the scalar potentials differ by at most a constant v k = v l + const, and (c) ψ k = ψ l .
Because the physical current density operator is system dependent, it is important to note that and that j m;ψ k is not independent of k. The theorem above is the point of departure for the density-functional formalism proposed in Ref. 9. In this setting, the magnetic vector potentials are restricted to a three-dimensional family, whereas the physical current densities belong to an infinite-dimensional function space, hinting at considerable redundancy of the formalism. Indeed, the allowed family of vector potentials is the subspace of Y * lin , defined in Eq. (28), corresponding to a = 0. Consequently, the results in Sec. III D can be viewed as sharpened versions of the above theorem. Firstly, the current density determines the physical magnetic moment J R . In turn, the triple (ρ, Λ, J R ) determines the pair (v, B) even without the assumption that ψ k is a ground state. An explicit expression for B is obtained from Eq. (92) above. Secondly, when the ground-state assumption is made, the physical current density and physical magnetic moment can be dropped altogether for general non-cylindrically symmetric densities since (ρ, L G ) is sufficient to determine (v, B).
For cylindrically symmetric ground-state densities, some information additional to (ρ, L G ) is needed to determine the component of B along the symmetry axis. Even in this case, however, the density pair (ρ, L G ) is sufficient to determine the ground-state wave function and therefore the exchange-correlation energy and expectation values of operators independent of the vector potential (i.e., π = p + A and r × π are not included).
In Ref. 9, a constrained-search functional that incorporates a new constraint that the wave function must reproduce a given canonical momentum is given, following a similar line of reasoning to earlier works by the same authors. It is unclear whether or not their constrained-search functional is intended to be a universal functional [34,50,51]. A density functional is said to be universal if it is independent of the external potentials (v, A). This definition is the natural generalization of the original usage by Hohenberg and Kohn [1] and common usage in textbooks and the literature (e.g., in Refs. 52 and 53, in a DFT setting where A = 0 is implicitly assumed). The term has also been used in the same manner by authors in the context of CDFT or BDFT (e.g., in Refs. 5 and 54).
Addressing this ambiguity in the use of the term 'universal', the authors of Ref. 9 have emphasized the fact that the external potentials are fixed in the Rayleigh-Ritz variation principle [9,55,56]. This answers criticism of their previous work based on inconsistency between the universality of their functional and the ordinary variation principle. Specifically, the derivations are only correct for a semi-universal functional F (ρ, L G , j, A) that is universal with respect to v and non-universal with respect to A or B.
Nonetheless, the term universal is still used to refer to their constrained-search functional (above their Eq. (45) in Ref. 9). It is possible that 'universal' in their terminology corresponds to 'semi-universal' in our terminology. However, insofar as the intent is to derive a density-functional that assigns the same intrinsic energy to a triple (ρ, L, j) irrespective of the external potential pair (v, A) their derivation is mistaken. Criticism directed at their previous work still applies, largely unaffected by the new constraint on the paramagnetic moment: The search domain of the functional defined in their Eq. (45), indicated by 'Ψ ρ,j (N, L) → ρ, j', must depend on the external vector potential and therefore be non-universal. Moreover, universality of a functional of (ρ, j) has been shown to be inconsistent with the usual Rayleigh-Ritz variation principle [34,50,51], and this holds also for a universal functional of (ρ, L, j).

A non-universal constrained-search functional of (ρ, j, A)
Let us be careful in tracking not only the overt dependence of the Hamiltonian on the external potential, but also the dependence of the various search domains on the potentials. We retain the choice of function spaces described above since they guarantee that the physical kinetic energy can split into finite canonical, paramagnetic, and diamagnetic terms and also that the magnetic moment is always finite. Because the formalism of Ref. [9] features a family of functionals parametrized by prescribed values of L, these conditions cannot be avoided.
We now note that the minimum energy for a given value of L can be written as where Vignale and Rasolt's constrained-search functional F VR (ρ, j m ) is defined as in Sec. II above. We reparametrize the outer minimization by setting j = j m + ρA, where R A = {(ρ, j)|j = j m + ρA, (ρ, j m ) ∈ X × Y }. Due to the choice of Y as a weighted L 1 space, we see that ρA ∈ Y . Consequently, the set R A is independent of A and, in fact, we have R A = X × Y . We can now write Thus, by defining the non-universal functional F PS (ρ, j, A) := F VR (ρ, j − ρA) we have Since this is simply a reformulation of Vignale and Rasolt's CDFT, it is also correct. We note, however, that the original formulation due Vignale and Rasolt yields a universal functional of (ρ, j m )-that is, the external potentials appear only in the integral pairings (j m |A) and (ρ|A 2 ). By contrast, after the reformulation, the external potential also appears in the search domain for F PS (ρ, j, A) and, via the constraint L = r × (j − ρA)dr, in the domain for the outer minimization. Different vector potentials A give rise to different search domains {ψ|ψ ∈ H 1 0 , ρ ψ = ρ, j m;ψ = j − ρA} for the functional F PS (ρ, j, A).
The functional F PS (ρ, j, A) can interpreted as a non-universal functional with potentials that remain fixed in the variation principle, rationalizing the dependence of the search domains on the external potential A. However, the advantage of such a formulation is unclear. Consider for example a light modification of the BDFT formalism due to Grayce and Harris. If universality is not required, we may also write Grayce and Harris' original work constructed a family of functionals parametrized by an arbitrary magnetic field B(r) [5]. Above we have modified this to a functional F GH;L (ρ; A) parametrized by A(r) which also incorporates the restriction on L. From the above nested minimization, it is clear that the pair (ρ, A), or alternatively the pair (ρ, B), alone determines the energy additional to the electrostatic interaction (ρ|v) with the external electric potential. In particular, the exchange-correlation energy is determined by (ρ, A). By contrast, the functional F PS (ρ, j, A) depends on the additional vector field j, yet j is formally redundant when (ρ, A) are given.

Non-standard properties of a universal functional of (ρ, j, A)
Suppose H(v k , A k ) has the ground state ψ k , with density triple (ρ k , j k , L k ) where j k = j m;k + ρ k A k . We let A N denote the set of all such (v, A)-representable ground-state density triples and A N,L is the subset with a given paramagnetic moment. To formulate a variation principle for such a density triple, we begin by rewriting the energy assigned to a state ψ l by the Hamiltonian H(v k , A k ), Alternatively, taking care to note that j l = j m;l + ρ l A l = j m;l + ρ l A k , unless A k = A l , we can also write Suppose now that there exists a universal functional F PS with the property for arbitrary ψ l → (ρ l , j l , L l ). The LDFT functional F (ρ, L), specialized to the present gauge restrictions (allowing the variable p to be dropped) is an example of such a functional. In what follows, we assume that there are also universal functionals with a non-trivial j-dependence. Identifying F PS (ρ l , j l , L l ) and taking affine combinations of the two alternative expressions now yields, where λ ∈ R is an arbitrary parameter andλ = 1 − λ. The middle expression defines a functional E λ of (v k , A k , ρ l , j l , L l ). We also have with equality when k = l. This inequality is just an alternative form of ψ l |H(v k , A k )|ψ l ≥ ψ k |H(v k , A k )|ψ k , which is just the Rayleigh-Ritz variation principle for wave functions. Identifying ψ k |H(v k , A k )|ψ k as E(v k , A k ), we now obtain the variation principle min (ρ l ,j l ,L l )∈AN (We here retain the index l to distinguish density triples based on the potentials they arise from. However, the minimization is of course performed with (ρ l , j l , L l ) ranging over all (v, A)-representable triples and not just a countable subset.) By constraining all magnetic moments to a particular value L k = L, we also obtain a variation principle for the lowest energy compatible with this constraint The case λ = 1 corresponds closely to the LDFT setting, at least if F is independent of j l . In this case, the energy functional being minimized is independent of j l . Hence, the minimization over j l is redundant, but does not lead to serious complications of the theory. The case λ = 0 corresponds to a setting akin to a physical CDFT. In this case, the appearance of the term −λ(ρ l |A k · A l ) in Eq. (141) gives the variation principle a non-standard form. The minimization over j l cannot be performed without access to the vector potential A l , which in general is different from the external potential A k that appears in the particular Hamiltonian H(v k , A k ) used to assign energies. In this regard, the constraint on the paramagnetic moment and the restriction of vector potentials to a three-dimensional space makes a crucial difference. Without the restrictions, a mapping (ρ l , j l ) → A l would be required, but there is no rigorous existence proof for this mapping and the mapping is unlikely to have a simple, practical form. With the restrictions, however, there is a simple, explicit mapping (ρ l , j l , L l ) → A l given by Eq. (92) above. A more subtle point is that the mapping (ρ l , j l , L l ) → A l makes no use of the ground-state assumption, and is not restricted to (v, A)-representable triples (ρ, j, L).
From the above discussion, we see that the non-standard variation principles Eqs. (143)-(144) delivers the correct ground-state energy and minimum energy compatible with L, respectively. However, the minimizing physical current density is not necessarily the same as that arising from the ground state. We now adapt a result from Laestadius and Benedicks [51] to the present setting: Theorem 11. Let H model a one-particle system. There exist potentials (v k , A k ) such that: (a) the minimizer (ρ l , j l , L l ) of Eq. (143) is not unique, and (b) the minimizer (ρ l , j l ) of Eq. (144) is not unique. Moreover, infinitely many of the minimizing j l are not the physical current of a ground state of H(v k , A k )-that is, j l = j k = j m;k + ρ k A k .
The proof of the above theorem follows every step in Sec.III.A of Ref. 51, which still applies and where we let the minimizing L l = ψ k |L|ψ k . Note that ψ k and the infinitely many current densities {j l } are denoted ψ 0 and {j ε } ε>0 , respectively, in Ref. 51.
Turning next to the interpretation of constrained-search functional in Ref. 9 as a universal functional, we note that it relies on a variation principle related to a functional E PS (v k , A k , ρ l , j l , L) = F (ρ l , j l , L) + (j l |A k ) + (ρ l |v k − 1 2 A 2 k ).
The relation to the functional E 1 studied above is Minimization of E PS with respect to (ρ l , j l ) corresponds to a seemingly natural variation principle. Unfortunately, this variation principle does not always deliver the correct energy. Again we adapt a result from Laestadius and Benedicks [51] to the present setting: Theorem 12. There exist potentials (v k , A k ) such that: There also exists L such that: min (ρ l ,j l )∈AN,L E PS (v k , A k , ρ l , j l , L) < E L (v k , A k ).
The essential steps of the proof of this statement is given in Sec.II of the work by Laestadius and Benedicks [51]. For the first part above, Eq. (147), we just need to add the remark given after Theorem 11 regarding L l . For the second part, Eq. (148), we set L = ψ k |L|ψ k , where ψ k is the ground state of H(v k , A k ).

VI. DISCUSSION AND CONCLUSION
A specialization of paramagnetic CDFT to uniform magnetic fields, represented by linear magnetic vector potentials, has been presented. As a result of the restriction on the vector potentials, the basic variables are reduced to the electron density, the canonical momentum, and paramagnetic moment. This theory, termed LDFT, admits a constrained-search formulation as well as a convex formulation analogous to Lieb's formulation of standard DFT. Unlike CDFT, LDFT also admits strong analogues of the HK theorem whenever the density is not rotationally symmetric: In this case, the ground-state triple (ρ, p, L G ) determines the scalar potential u and the linear magnetic vector potential A. Moreover, the ground-state pair (ρ, Λ) of density and intrinsic magnetic moment determines the scalar potential and magnetic field (u, B). For a large total system, consisting of well-separated subsystems, the magnetic part of the HK-mapping can be written in a simple explicit form involving the inverse moment-of-inertia tensor.
Many issues in CDFT find their analogue, in a simplified setting, in LDFT. For example, both the CDFT and LDFT universal functionals can be explored using Lieb optimization. The optimization over the scalar potential is known to be computationally and numerically demanding. In practice, the space of potentials X * must be replaced by a subset spanned by a finite number of basis functions, typically Gaussian-type functions. The optimization over vector potentials A ∈ Y * exacerbate these problems. By contrast, in the LDFT setting, the optimization over linear vector potentials A ∈ Y * lin requires only six additional variational parameters and is straightforward compared with the scalar potential optimization. Another example is that gauge invariance of the exchange-correlation energy leads to a CDFT functional that depends on the vorticity, whereas gauge-shift invariance leads to an LDFT functional that depends on the intrinsic magnetic moment. In the special case that the current density corresponds to rigid rotation, the vorticity is constant throughout space and related to intrinsic magnetic moment through Λ = 1 2 Qν. In more realistic cases, the two quantities have rather different properties since ν is a local, intensive quantity and Λ is a global quantity that grows superlinearly with system size. This difference manifests itself in, for example, the failure of simple LDA-like forms to satisfy additive separability. Even with a more flexible ansatz, it is likely to be challenging to construct additively separable approximations. Here Lieb optimization has a large potential to be useful in uncovering the Λ dependence of the exchange-correlation energy.
A third issue is that of the possibility of replacing the paramagnetic current density by the physical current in CDFT. The LDFT analogue is to replace the paramagnetic moment by the physical magnetic moment. In many respects this would be appealing since both vorticity and intrinsic magnetic moment have practical drawbacks. The challenges that have been pointed out for CDFT-namely, the dilemma that constrained-search formulations suffer from either nonuniversality or tension with the variation principle-carry over directly a hypothetical LDFT with physical magnetic moment. Under the restriction of uniform magnetic fields and completely gauge-fixed vector potentials, a hybrid approach was recently suggested by Pan and Sahni, where the triple (ρ, j, Λ) of density, physical current density, and paramagnetic moment serves as the basic variables. However, in light of the HK-like results proved in this work, the Pan and Sahni formulation is seen to be heavily redundant. The pair (ρ, Λ) is sufficient to determine the exchangecorrelation energy. When ρ is not rotationally symmetric, the pair (ρ, Λ) also determines the external scalar potential and magnetic field. If a HK-like result for rotationally symmetric systems is desired, the triple (ρ, J, Λ) of density, physical magnetic moment, and intrinsic magnetic momentum is sufficient.