BBUILDING KOHN-SHAM POTENTIALSFOR GROUND AND EXCITED STATES
LOUIS GARRIGUE
Abstract.
We analyze the inverse problem of Density Functional The-ory using a regularized variational method. First, we show that given k and a target density ρ , there exist potentials having k th excited mixedstates which densities are arbitrarily close to ρ . The state can be cho-sen pure in dimension d = 1 and without interactions, and we providenumerical and theoretical evidence consistently leading us to conjecturethat the same pure representability result holds for d = 2 , but thatthe set of potential-representable pure state densities is not dense for d = 3 . Finally, we present an inversion algorithm taking into accountdegeneracies. Introduction
In 1965, Kohn and Sham postulated the existence of effective one-body po-tentials which would replace the electronic interaction while keeping the sameground state density, and stated their relations to the exchange-correlationfunctionals [26]. Physical quantities of this new effective non-interacting sys-tem provide approximations of the exact ones. This led to the developementof very successful techniques enabling to predict properties of microscopicsystems in quantum chemistry and physics. The existence of such a potentialproducing a prescribed ground state density ρ is called the v -representabilityproblem, and its search is the inverse problem of Density Functional Theory.There are few works addressing the mathematical aspects of this problem,although several numerical studies were carried out. In [34], Lieb provedthat any density can be approximately represented by a ground mixed statein some external potential v , and introduced a dual variational method en-abling to find the Kohn-Sham potential. The ground state v -representabilityproblem was studied by variational methods in the cases of classical DFTat positive temperature [9] and for quantum lattices [10]. Moreover, study-ing excited states is important for DFT [13], and we do not know any workdedicated to the systematic search of corresponding inverse potentials.In this document, we address the problem of v -representability in thequantum case at zero temperature, with ground or excited states, in pureand mixed settings, both theoretically and numerically.In a first part, we present a mathematical investigation. As shown byLieb [34], the exact inverse potential of a density ρ (cid:62) maximizes thefunctional v (cid:55)→ E (0) ( v ) − (cid:82) vρ in ( L p + L ∞ )( R d , R ) , where E (0) ( v ) is the N -particle ground state energy and p is defined in (1). Nevertheless, thisfunctional is not locally coercive in this space as we will see in Section 3.2.3. Date : January 5, 2021. a r X i v : . [ m a t h - ph ] J a n L. GARRIGUE
To circumvent this ill-posedness, we regularize the problem by discretizingthe space of potentials, more precisely we restrict our attention to potentialsof the form (cid:80) i v i α i where the α i are fixed weight functions and v i ∈ R arethe parameters. The discretization amounts to integrating the problematicshort-distance degrees of freedom, and implements an ultra-violet cut-off.Our approach enables to show that for any k ∈ N , the functionals v (cid:55)→ E ( k ) ( v ) − (cid:82) vρ are coercive, where E ( k ) ( v ) denotes the k th excited stateenergy, and this implies the approximate representability of densities by k th mixed excited states, with arbitrary precision. For ( d, w ) = (1 , , where w isthe interaction, we show that we can take a pure state. Correspondingly, wedefine pseudo-discrete regularized Levy-Lieb and Lieb functionals by relaxingthe condition ρ Ψ = ρ to (cid:82) ρ Ψ α i = (cid:82) ρα i for any i , a similar approach wasapplied to optimal transport in [2, 12], where the numerical efficiency seemspromising.Numerically, computing the inverse ground state Kohn-Sham potential re-ceived significant attention, in [1, 4, 20, 24, 25, 27, 38, 39, 41–43] using the dualformulation, in [24] using the PDE-constrained optimization, and in [40] us-ing derivatives of the Levy-Lieb functional. However, degeneracy issues arecritical in the dual approach, except when ( d, k, w ) = (1 , , , and they wereno taken into account in the literature to the best of our knowledge. Indeed,the standard algorithm blocks when eigenvalues cross or when the inversepotential is degenerate. Hence in a second part, we present an algorithmwhich converges to a potential having a k th excited mixed state with thetarget density, the first part of this document proving the existence of thispotential. With w = 0 , numerical results indicate that for d = 2 , densi-ties are representable by k th excited pure states of potentials, whereas thisdepends on the targeted density for d = 3 . We also numerically remarkthat degeneracies are generic for inverse potentials, and as in the SCF pro-cedure [8] that perturbation of target densities does not lift degeneracies.Finally, we confirm the study [18], which indicates that for excited states,many potentials lead to the same density. Acknowledgement.
I warmly thank Mathieu Lewin, for having advisedme during this work, and Éric Cancès for useful discussions. This projecthas received funding from the European Research Council (ERC) under theEuropean Union’s Horizon 2020 research and innovation programme (grantagreements MDFT No 725528 and EMC2 No 810367).2.
Properties of the dual problem
Definitions.
Let d ∈ N \ { } and let Ω ⊂ R d be a (bounded or un-bounded) open set with Lipschitz boundary, representing the space in whichour quantum system lives. We do not consider spin degrees of freedom butour results can be extended in this way without complications. We define p = 1 if d = 1 , p > if d = 2 , p = d/ if d (cid:62) . (1)In all this work, we consider an even non-negative interaction potential w ∈ ( L p + L ∞ )( R d , R + ) . We take external electric potentials v ∈ ( L p + L ∞ )(Ω , R ) , UILDING INVERSE POTENTIALS 3 and consider the self-adjoint N -particle Schrödinger operator H N ( v ) := N (cid:88) i =1 − ∆ i + (cid:88) (cid:54) i For ρ (cid:62) with (cid:82) ρ = N and √ ρ ∈ H , theexact ground and excited Levy-Lieb (or pure) and Lieb (or mixed) functionals[29, 33–35], are F ( k ) ( ρ ) := sup A ⊂ H a (Ω N )dim C A = k inf Ψ ∈ A ⊥ ρ Ψ = ρ (cid:104) Ψ , H N (0) Ψ (cid:105) , (4) F ( k ) mix ( ρ ) := sup A ⊂ H a (Ω N )dim C A = k inf Γ ∈S N mix ( A ⊥ , Ω) ρ Γ = ρ Tr H N (0) Γ . L. GARRIGUE For any k , F ( k ) mix is convex and lower semi-continuous [35]. Moreover, we have F ( k ) mix (cid:54) F ( k ) , F ( k ) mix (cid:54) F ( k +1) mix and F ( k ) (cid:54) F ( k +1) . We know that in the groundstate case k = 0 , they are finite and have optimizers [34], and that theyenable to compute the energy in the sense that E (0) ( v ) = inf ρ ∈ L (Ω , R + ) √ ρ ∈ H , (cid:82) ρ = N (cid:18) F (0) ( ρ ) + (cid:90) vρ (cid:19) . As noted by Lieb in [35], for k (cid:62) we cannot recover E ( k ) ( v ) using F ( k ) ( ρ ) or F ( k ) mix ( ρ ) , or even any functional of ρ , because this would form a convexfunctional while E ( k ) is not so.2.3. The dual problem. For k ∈ N , and for v ∈ L p + L ∞ where p is as in(1), the dual functional is G ( k ) ρ ( v ) := E ( k ) ( v ) − (cid:90) Ω vρ, verifying sup v ∈ ( L p + L ∞ )(Ω , R ) G ( k ) ρ ( v ) = F ( k ) mix ( ρ ) , as showed in [35]. To prepare the approximate representability of groundand excited densities by potentials, we want to explore the maximization of G ( k ) ρ .2.4. Solution of the local dual problem. The local first order problemis to find the optimal direction(s) in which the functional G ( k ) ρ increases themost, that is solving sup u ∈ L p + L ∞ || u || Lp + L ∞ =1 + δ v G ( k ) ρ ( u ) . Here + δ v G ( k ) ρ denotes the half Gateaux differential of G ( k ) ρ , defined and de-rived in [16, Theorem 1.6] and (24). We now show that this linearized prob-lem can be “solved”, that is transformed to a simple low-dimensional problem.For a real vector space Q R ⊂ L a (Ω N , R ) formed by real wavefunctions, theset of matrices S ( Q R ) := (cid:110) Γ ∈ R dim Q R × dim Q R , Γ = Γ T (cid:111) is identified to the set of real mixed states on Q R . For k ∈ N and potentials v ∈ V ( k ) N,∂ , we define the integers m k , M k ∈ N by E ( m k − ( v ) < E ( m k ) ( v ) = · · · = E ( k ) ( v ) = · · · = E ( M k ) ( v ) < E ( M k +1) ( v ) , (5)with E ( − ( v ) := −∞ by convention. We denote by Ker R ( H N ( v ) − E ( k ) ( v )) the real vector space of real k th excited eigenfunctions. Proposition 2.1 (The local problem) . Take w (cid:62) . Let h be as in (1) , p (cid:62) h and q := p/ ( p − , and take v ∈ V ( k ) N,∂ . We have sup u ∈ L p + L ∞ || u || Lp + L ∞ =1 + δ v G ( k ) ( u ) = max Q ⊂ Ker R ( H N ( v ) − E ( k ) ( v ))dim R Q = M k − k +1 min Γ ∈S ( Q )Γ (cid:62) , Tr Γ=1 || ρ Γ − ρ || L q (Ω) , (6) UILDING INVERSE POTENTIALS 5 and the supremum is attained by u ∗ = (cid:12)(cid:12)(cid:12)(cid:12) ρ Γ ∗ − ρ || ρ Γ ∗ − ρ || L q (cid:12)(cid:12)(cid:12)(cid:12) q − sgn( ρ Γ ∗ − ρ ) , (7) where Γ ∗ is an optimizer of the right hand side of (6) . In particular, when Ker (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) =: C Ψ ( k ) ( v ) is non-degenerate, G ( k ) ρ is Fréchet differentiable at v and the problemsup u ∈ L p + L ∞ || u || Lp + L ∞ =1 d v G ( k ) ρ ( u ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ Ψ ( k ) ( v ) − ρ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) L q is solved in the unique direction (7), where ρ Γ ∗ = ρ Ψ ( k ) ( v ) .2.5. Optimality in the dual problem. Next, we analyze the optimalityconditions. For our problem of searching a potential producing a prescribeddensity, the following result shows that we have to search among the maxi-mizers of G ( k ) ρ . Theorem 2.2 (Optimality in the dual problem) . Take w (cid:62) , take a density ρ ∈ L ( R d ) , ρ (cid:62) , (cid:82) ρ = N , √ ρ ∈ H ( R d ) , and consider a binding v ∈ V ( k ) N,∂ . i ) The following assertions are equivalent • there is a k th excited mixed state Γ of v such that ρ Γ = ρ • v is a local maximizer of G ( k ) ρ • v is a global maximizer of G ( k ) ρ ii ) If v maximizes G ( k ) ρ , then it maximizes G ( (cid:96) ) ρ for all (cid:96) ∈ { m k , . . . , k } .Moreover, if k (cid:62) ( M k + m k ) / , then ρ Ψ = ρ for any normalized Ψ ∈ Ker (cid:0) H N ( v ) − E ( k ) (cid:1) . If v is a local minimizer, then k > ( M k + m k ) / . iii ) If v maximizes G ( k ) ρ and dim Ker (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) ∈ { , } , then v has a k th excited pure state Ψ such that ρ Ψ = ρ . iv ) For d = 1 and w = 0 , if v maximizes G ( k ) ρ , then there exists a purestate Ψ ∈ Ker (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) such that ρ Ψ = ρ . In iii ) , we take E ( (cid:96) ) ( v ) := −∞ for (cid:96) (cid:54) − by convention. When k = 0 and p > max(2 d/ , , the maximizer is unique by the Hohenberg-Kohn theorem,and the equivalences do not need to assume v ∈ V (0) N,∂ . Since pure states aremixed states, when we search pure states of v representing ρ , we also needto maximize G ( k ) ρ . Once the set V max ρ,k of maximizers is found, one can finallycompute inf v ∈V max ρ,k min Ψ ∈ Ker (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) (cid:82) | Ψ | =1 || ρ Ψ − ρ || L q , (8)which vanishes if and only if ρ is pure-state representable.Moreover, we think that G ( k ) ρ has no local minimum and that for maxi-mizing v ’s, k = m k .In case where w = 0 and k = 0 , the maximizing potential v =: v ks ( ρ ) iscalled the (mixed) Kohn-Sham potential [26] for ρ . We use the term “Kohn-Sham potential” for any k ∈ N . When k = 0 and ρ is a ground state density L. GARRIGUE of H N ( u ) for w = |·| − , then v ks ( ρ ) − u − ρ ∗ |·| − is called the exchange-correlation potential in the physics and quantum chemistry literature.2.6. Ill-posedness of the dual problem. Let us search for a maximizingpotential of G (0) ρ , we consider a maximizing sequence v n , if we are able toprove that v n converges weakly to some v ∈ L p + L ∞ , v would be a maximizerby weak upper semi-continuity of G (0) ρ . But we are not even able to provethat v n weakly converges locally.Moreover we now show that G ( k ) ρ is ill-posed in the sense that it is notlocally coercive in L p spaces. Let ρ ∈ L ( R d , R + ) having mass (cid:82) ρ = N be a target density which we want to represent by a potential. First of all,as a consequence of [34, Theorem 3.8] and of G (0) ρ (cid:54) G ( k ) ρ , one needs toassume that √ ρ ∈ H ( R d ) , otherwise G ( k ) ρ is not bounded from above. Nowtake p > , take ρ continuous at the origin, and take a positive potential v ∈ ( L ∩ L p )( R d ) with compact support. Consider the sequence v n ( x ) := n d v ( nx ) . Then || v n || pL p = n d ( p − (cid:82) v p → + ∞ because p > but E ( k ) ( v n ) = 0 because v n (cid:62) , and (cid:82) v n ρ → ρ (0) (cid:82) v is bounded. We remark that this lastcounterexample does not hold when d = 1 , because then p = 1 .3. Regularization We saw in Theorem 2.2 that to represent a density ρ with pure or mixedstates of a potential, we need to maximize G ( k ) ρ , but we also saw that thisfunctional is ill-posed in L p spaces. Hence we regularize it in this section tomake it coercive.3.1. Pseudo-discrete regularizations of Levy-Lieb and Lieb. We nowrelax the density constraint. Let us consider a subset I ⊂ N and a set α = ( α i ) i ∈ I of weight functions forming a partition of unity for Ω , that is (cid:80) i ∈ I α i = Ω , where α i ∈ L ∞ (Ω , R + ) . For r ∈ (cid:96) ( I, R + ) ∩ (cid:8)(cid:80) i ∈ I r i = N (cid:9) ,we introduce the regularized Levy-Lieb and Lieb functionals F α , ( k ) ( r ) := sup A ⊂ H a (Ω N )dim C A = k inf Ψ ∈ A ⊥ (cid:82) α i ρ Ψ = r i ∀ i ∈ I (cid:104) Ψ , H N (0) Ψ (cid:105) ,F α , ( k ) mix ( r ) := sup A ⊂ H a (Ω N )dim C A = k inf Γ ∈S N mix ( A ⊥ , Ω) (cid:82) α i ρ Ψ = r i ∀ i ∈ I Tr H N (0) Γ , and we define them to be + ∞ when for any A ⊂ H a (Ω N ) such that dim A = k , the minimizing sets are empty. We know that F α , (0) mix is convex. Considernow the assumption lim R → + ∞ (cid:88) i ∈ I supp α ni ∩ B c R (cid:54) = ∅ r i = 0 . (9) Theorem 3.1 (Existence of minimizers in the ground states case) . Take w (cid:62) and r ∈ (cid:96) ( I, R + ) such that (cid:80) i ∈ I r i = N , and α = ( α i ) i ∈ I a partitionof unity for Ω . Under the tightness condition (9) , F α , (0) ( r ) and F α , (0) mix ( r ) have at least one minimizer when they are finite. UILDING INVERSE POTENTIALS 7 For a given ρ ∈ L (Ω , R + ) , we define r ρ := (cid:18)(cid:90) α i ρ (cid:19) i ∈ I ∈ (cid:96) ( I, R + ) . This sequence contains the partial information on the density ρ which we aregoing to retain. Since the optimizing set of F ( k ) ( ρ ) is included in the one of F α , ( k ) ( r ρ ) , for any ρ (cid:62) with √ ρ ∈ H (Ω) and α as defined above we have F α , ( k ) ( r ρ ) (cid:54) F ( k ) ( ρ ) and F α , ( k ) mix ( r ρ ) (cid:54) F ( k ) mix ( ρ ) . In particular, F α , (0) ( r ρ ) and F α , (0) mix ( r ρ ) are finite.For k = 0 , our approximate Levy-Lieb and Lieb functionals converge tothe exact ones when the integrated weights tend to carry all the informationof the density. Theorem 3.2 (Convergence to the exact model) . Take (cid:54) w ∈ L p + L ∞ with p as in (1) . Let Ω be a connected open set. Consider a density ρ ∈ L (Ω , R + ) such that √ ρ ∈ H (Ω) and (cid:82) ρ = N . We assume that α n =( α ni ) i ∈ I n , where α ni ∈ L ∞ (Ω) , is a sequence of weights forming a partition ofunity for Ω , and such that for any f ∈ C ∞ c (Ω) , we have inf g n ∈ Span( α ni ) i ∈ In || f − g n || ( L p + L ∞ )(Ω) −→ (10) when n → + ∞ . We also assume that lim r → + ∞ sup n ∈ N (cid:88) i ∈ I supp α ni ∩ B c r (cid:54) = ∅ (cid:90) ρα ni = 0 . (11) Then lim n → + ∞ F α n , (0) ( r ρ ) = F ( ρ ) , lim n → + ∞ F α n , (0) mix ( r ρ ) = F mix ( ρ ) . Let Ψ n be a sequence of approximate minimizers for F α n , (0) ( r ρ ) , that is,such that E (Ψ n ) (cid:54) F α n , (0) ( r ρ ) + (cid:15) n where (cid:15) n → when n → + ∞ and (cid:82) α i ρ Ψ n = (cid:82) α i ρ for any i ∈ I . Then Ψ n → Ψ exact strongly in H (Ω N ) upto a subsequence, where Ψ exact is a minimizer of F (0) ( ρ ) . If Γ n is a sequenceof approximate minimizers for F α n , (0) mix ( r ρ ) , then Γ n → Γ exact strongly in thekinetic energy space S , up to a subsequence, where Γ exact is a minimizerof F mix ( ρ ) . The space S , is the set of operators A of L ( R dN ) endowed with thenorm || A || S , = Tr (cid:12)(cid:12)(cid:12) ( − ∆ + 1) A ( − ∆ + 1) (cid:12)(cid:12)(cid:12) .Our assumption (11) is used to control the decay at infinity. If all the α ni have a compact support of diameter bounded by δ independent of i and n ,then (cid:88) i ∈ I supp α ni ∩ B c r (cid:54) = ∅ (cid:90) ρα ni (cid:54) (cid:90) | x | (cid:62) r − δ ρ −→ r → + ∞ and (11) is satisfied. In [2, (3.3.4)], the authors use an inequality condition,simpler than (11). L. GARRIGUE If α ni = Ω ni is a sequence of partitions of Ω = ∪ i ∈ N Ω ni where Ω ni areconvex, and sup i ∈ N diam Ω ni → when n → + ∞ , then the assumption (10)is verified by Lemma 6.3 below. Assumption (11) is verified as well. Againby Lemma 6.3, in this case we have an explicit bound on the convergence ofdensities || ρ − ρ Ψ n || (cid:0) L ∩ L q (cid:1) (Ω) (cid:54) c d (cid:0) ||√ ρ || H + sup n ∈ N (cid:12)(cid:12)(cid:12)(cid:12) √ ρ Ψ n (cid:12)(cid:12)(cid:12)(cid:12) H (cid:1) sup i ∈ N diam Ω ni , where c d only depends on d , and q is as in (25). Note that (cid:12)(cid:12)(cid:12)(cid:12) √ ρ Ψ (cid:12)(cid:12)(cid:12)(cid:12) H and (cid:12)(cid:12)(cid:12)(cid:12) √ ρ Ψ n (cid:12)(cid:12)(cid:12)(cid:12) H are controlled by F ( ρ ) due to the Hoffmann-Ostenhof inequality.A typical choice for the α ni is given by the partition of unity finite elementmethod [3, 37],3.2. Regularization of the dual problem. Correspondingly to the previ-ous part, we change the exact model by discretizing the space of potentials.We consider a sequence of weights α = ( α i ) i ∈ I , take r ∈ (cid:96) ( I, R + ) . The dualproblem is the maximization of G ( k ) r, α ( v ) := G ( k ) ρ (cid:32)(cid:88) i ∈ I v i α i (cid:33) = E ( k ) (cid:32)(cid:88) i ∈ I v i α i (cid:33) − (cid:88) i ∈ I v i r i , over the space (cid:96) ∞ ( I, R ) of potential coefficients v = ( v i ) i ∈ I . As in the exactcase, we have E (0) (cid:32)(cid:88) i ∈ I v i α i (cid:33) = inf r ∈ (cid:96) ( I, R + ) (cid:80) i ∈ I r i = N (cid:32) F α , (0) mix ( r ) + (cid:88) i ∈ I v i r i (cid:33) , sup v ∈ (cid:96) ∞ ( I, R ) G ( k ) r, α ( v ) = F α , ( k ) mix ( r ) , as in the exact models, and by the same proofs. Again by the same proof asfor the lower semi-continuity of the exact Lieb functional [34, Theorem 3.6], G (0) r, α is weakly upper semi-continuous in the (cid:96) ∞ ( I, R ) topology. Moreover, if H N (cid:0)(cid:80) i ∈ I v i α i (cid:1) has a k th excited state Ψ v , then G ( k ) r, α ( v ) = E (Ψ v ) + (cid:88) i ∈ I v i (cid:18) − r i + (cid:90) ρ Ψ v α i (cid:19) . Gauge invariance. The gauge we are dealing with is the choice of areference for energies, corresponding to the transformation V → V + c fora constant c ∈ R . The exact dual functional V (cid:55)→ E ( k ) ( V ) − (cid:82) V ρ is gaugeinvariant, and since we want our approximate functional to be so as well, weare naturally led to take (cid:88) i ∈ I α i = 1 on Ω , (cid:88) i ∈ I r i = N. (12)The last condition is of course fulfilled for r = r ρ , which is the interestingsituation. UILDING INVERSE POTENTIALS 9 Remark 3.3. Let us explain why the previous conditions (12) are necessaryto ensure gauge invariance. Let v ∈ (cid:96) ∞ ( I, R ) be such that H N (cid:0) (cid:80) i ∈ I v i α i (cid:1) has a k th excited state, which we denote by Ψ v . Take c ∈ R , we have E ( k ) (cid:32)(cid:88) i ∈ I ( v i + c ) α i (cid:33) (cid:54) E Σ i ∈ I ( v i + c ) α i (Ψ v ) = E ( k ) (cid:32)(cid:88) i ∈ I v i α i (cid:33) + c (cid:90) (cid:88) i ∈ I α i ρ Ψ v , and hence G ( k ) r, α ( v + c ) (cid:54) G ( k ) r, α ( v ) + c (cid:32) − (cid:88) i ∈ I r i + (cid:90) Ω ρ Ψ v (cid:88) i ∈ I α i (cid:33) . To have a gauge invariant theory, we want to have (cid:90) Ω (cid:32) N − (cid:88) i ∈ I r i − (cid:88) i ∈ I α i (cid:33) ρ Ψ v = 0 , otherwise G ( k ) r, α ( v + c ) → −∞ for c → + ∞ or c → −∞ . This requirementshould not depend on v , hence we need (cid:80) i ∈ I r i = N (cid:80) i ∈ I α i a.e on Ω . Weare thus naturally led to assume (12) . Uniqueness. A Hohenberg-Kohn theorem adapted to our situationshows that the multivalued map v (cid:55)→ r ρ Ψ v , where Ψ v is a ground stateof H N (cid:0) (cid:80) v i α i (cid:1) , is essentially injective. Hence if G (0) r, α has a maximum, it isunique. Theorem 3.4 (Hohenberg-Kohn) . Let Ω ⊂ R d be an open and connectedset, with Dirichlet any boundary conditions. Let p > max(2 d/ , , and takean interaction w ∈ ( L p + L ∞ )( R d ) . Let v, u ∈ (cid:96) ∞ ( I, R ) and α = ( α i ) i ∈ I where α i ∈ L ∞ (Ω , R + ) , be such that H N (cid:0) (cid:80) i ∈ I v i α i (cid:1) and H N (cid:0) (cid:80) i ∈ I u i α i (cid:1) have at least one ground state each, which we respectively denote by Ψ v and Ψ u . If (cid:82) α i ρ Ψ v = (cid:82) α i ρ Ψ u for any i ∈ I , then v = u + c for some constant c ∈ R . The proof follows from the standard Hohenberg-Kohn theorem [15, 23] inthe form of [17, Theorem 2.1].3.2.3. Coercivity. The main goal of this section is to recover coercivity forthe discretized dual problem, in order to make it well-posed.If there is some i ∈ I such that r i = 0 , denoting by e i the i th degree offreedom of the potentials, when c → + ∞ we expect that G ( k ) r, α ( v + ce i ) → E ( k ) D (cid:0) (cid:80) j (cid:54) = i v j α j (cid:1) , where E ( k ) D (cid:0) (cid:80) j (cid:54) = i v j α j (cid:1) is finite and is the k th excited stateenergy of the system living in Ω \ supp α i with Dirichlet boundary conditions.This shows that r i > ∀ i ∈ I is a necessary condition for G ( k ) r, α to be coercive.We define c Ω := − E ( k ) (0) N (cid:54) , where E ( k ) (0) is the k th energy of N interacting particles without externalpotential. It respects E ( k ) ( c Ω Ω ) = 0 and it is non-positive because w (cid:62) .It vanishes when Ω = R d for instance. We can choose the gauge we want,so we will take potentials v such that E ( k ) (cid:0) (cid:80) i ∈ I v i α i (cid:1) = 0 for convenience.Our variational space of potentials can thus be (cid:110) v ∈ (cid:96) r ( I, R ) (cid:12)(cid:12) E ( k ) (cid:0) (cid:88) v i α i (cid:1) = 0 (cid:111) , where || v || (cid:96) r := (cid:80) i ∈ I | v i | r i .Now we can state our main result for the discretized model. Theorem 3.5 (Well-posedness of the dual problem) . Take a non-negativeinteraction w ∈ ( L p + L ∞ )( R d , R + ) where p is as in (1) . • (Coercivity) Let α be a partition of unity of Ω , with α i ∈ L ∞ ( R d , R + ) ,such that we have R > for which (supp α i ) \ ∪ j ∈ I,j (cid:54) = i supp α j contains a ball of radius R , uniformly in i ∈ I . Let r ∈ (cid:96) ( I, R + ) be suchthat (cid:80) i ∈ I r i = N and r i > for all i ∈ I . For any v ∈ (cid:96) r ( I, R ) such that E ( k ) (cid:0) (cid:80) i ∈ I v i α i (cid:1) = 0 , we have G ( k ) r, α ( v ) (cid:54) − min (cid:32) , (cid:80) v i (cid:62) c Ω r i (cid:80) v i Our bound (13) does not pass to the continuous model because then R → and c R → + ∞ . (iv) The pair ( v, Γ v ) is a saddle point of the Lagrangian L ( v, Γ) = E (Γ) + (cid:88) i ∈ I v i (cid:18) − r i + (cid:90) ρ Γ α i (cid:19) (14)and v is a Lagrange multiplier. (v) We recall that the exact existence was proved for the quantum theoryon a Z d lattice [10], and in the classical case at positive temperature [9].3.3. Building Kohn-Sham potentials. The problem of v -representabilityis, given a density ρ ∈ L (Ω , R + ) , (cid:82) ρ = N , √ ρ ∈ H , and |{ ρ = 0 } ∩ Ω | = 0 ,to find a potential v having a k th excited state Ψ v respecting ρ Ψ v = ρ . Wewill call it the inverse potential. When w = 0 and k = 0 , it is called theKohn-Sham potential [26].3.3.1. The mixed states case. In the mixed states setting (at zero temper-ature) and using the Bishop-Phelps theorem, Lieb showed in [34, Theorem3.10, Theorem 3.11, Theorem 3.14], that any such ρ can be approached toany precision in L ∩ L d/ ( d − by a v -representable ground mixed state den-sity. We can state a similar result for any k using our variational approach. Corollary 3.6 (Constructive approximate v -representability in the mixedstates setting) . Let Ω ⊂ R d be a possibly unbounded open set with Lipschitzboundary. Let ρ ∈ L (Ω , R + ) be such that √ ρ ∈ H (Ω) and |{ ρ = 0 } ∩ Ω | =0 , and k ∈ N . There exists a sequence v n ∈ L ∞ (Ω , R ) with compact supportsuch that H N ( v n ) has a mixed k th excited state Γ v n with E (Γ v n ) (cid:54) F ( k ) mix ( ρ ) and ρ Γ vn → ρ strongly in ( L ∩ L q )(Ω) , where q is as in (25) .If moreover k = 0 , Γ v n → Γ ∞ strongly in S , up to a subsequence, where Γ ∞ is a minimizer of F (0) mix ( ρ ) . Furthermore, √ ρ Γ vn → √ ρ strongly in H (Ω) ,and E (Γ v n ) → F (0) mix ( ρ ) . Although the existence part of Theorem 3.5 holds only for bounded opensets Ω , Corollary 3.6 holds when Ω is unbounded. The proof uses Theo-rem 3.5 on a sequence of growing bounded sets Ω n with well chosen weightfunctions α ni . We conjecture that if there exists a potential which exactlyproduces ρ , this sequence v n converges to this exact inverse potential, in asuitable sense.However, any ρ (cid:62) such that √ ρ ∈ H and (cid:82) ρ = N is not necessarilyexactly representable, one would need more constraints on it. For instance if ρ decreases more than exponentially, then the Kohn-Sham sequence v n wouldnot converge, it would become very large as | x | → + ∞ . It will neverthelessprobably converge locally.A consequence of Corollary 3.6 and Theorem 2.2 iv ) is the density ofnon-interacting pure v -representable densities. Corollary 3.7. Take d = 1 and w = 0 . The set (cid:26) ρ Ψ v (cid:12)(cid:12) v ∈ V ( k ) N,∂ , Ψ v ∈ Ker (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) , (cid:90) | Ψ v | = 1 (cid:27) is dense in the set of densities (cid:8) ρ ∈ L ( R , R + ) (cid:12)(cid:12) (cid:82) ρ = N (cid:9) . The pure states case. For any density ρ , F ( k ) mix ( ρ ) < F ( k ) ( ρ ) impliesthat ρ is not v -representable with pure k th excited states. To continue, wemake a conjecture. Conjecture 3.8 (Continuity of the Levy-Lieb and Lieb functionals) . Takedensities ρ, ρ n ∈ L ( R d , R + ) such that √ ρ, √ ρ n ∈ H ( R d ) . If √ ρ n → √ ρ in H ( R d ) , then F (0) ( ρ n ) → F (0) ( ρ ) and F (0) mix ( ρ n ) → F (0) mix ( ρ ) . Conjecture 3.8 would imply that the set of pure-state v -representableground densities is not dense in L ( R d , R + ) when d (cid:62) . Indeed, con-sider a density ρ such that F (0) mix ( ρ ) < F (0) ( ρ ) , the existence of such den-sities is presented in [34, Theorem 3.4 (ii)] for d = 3 but similar exampleshold for any d (cid:62) . Then by Conjecture 3.8 there exists R > such that F (0) mix ( χ ) < F (0) ( χ ) for any positive √ χ ∈ B R ( √ ρ ) , where we consideredthe ball B R ( √ ρ ) ⊂ H ( R d , R ) . Hence B R ( √ ρ ) is an open (in the set ofnon-negative square functions) set of densities which are not pure-state v -representable.However, with a different method which is not variational, it might stillbe possible to represent those densities, with excited states. As presentedin [14] for instance, the inverse potential can be seen as a Lagrange multipliercorresponding to the Euler-Lagrange equation of the Levy-Lieb functional.We give here a result for the discretized problem which only works for N = 1 . Theorem 3.9 (Pure excited v -representability, N = 1 ) . Take N = 1 , let Ω ⊂ R d be a bounded open domain with Lipschitz boundary, consider a finitepartition of unity ( α i ) i ∈ I for Ω , and r ∈ (cid:96) ( I, R + ) , (cid:80) i ∈ I r i = 1 and suchthat r i > for any i ∈ I . There exist v ∈ (cid:96) ∞ ( I, R ) and a pure one-particleground or excited state Ψ r ∈ H a (Ω) of − ∆ + (cid:80) i ∈ I v i α i such that for all i ∈ I , (cid:82) ρ Ψ r α i = r i . Applying the last result for an increasing sequence of α ’s ( α n ⊂ α n +1 ), weget the corresponding approximate representability, as we obtained Corollary3.6. For N = 1 the limit potential must be Bohm’s potential ∆ √ ρ/ √ ρ andthe state must be the ground state. We conjecture that Theorem 3.9 holdsfor any N , a sufficient condition being that minimizers Ψ of our approximateLevy-Lieb functionals have |{ Ψ = 0 }| = 0 . Conjecture 3.10. Any minimizer Ψ of F α , (0) ( r ) satisfies (cid:12)(cid:12)(cid:12)(cid:110) x ∈ R dN (cid:12)(cid:12) Ψ( x ) = 0 (cid:111)(cid:12)(cid:12)(cid:12) = 0 . This conjecture is related to a unique continuation property. In theHohenberg-Kohn theorem, one considers minimizers of the energy E v , re-specting Schrödinger’s equation, and this implies |{ Ψ = 0 }| = 0 by uniquecontinuation [17]. Here, this is a converse property in the sense that we con-sider minimizers of F α , (0) ( r ) , and the property |{ Ψ = 0 }| = 0 of minimizers,that we want to show, would imply that they respect Schrödinger’s equation(see the proof of Theorem 3.9).We remark that Conjecture 3.10 also implies that for any r ∈ (cid:96) ( I, R + ) , κ − ( r ) is a manifold, because d Ψ κ is then surjective and one can applythe preimage theorem [45, Theorem 73.C p556]. Since d Ψ κ is closed, this UILDING INVERSE POTENTIALS 13 mathematical framework can be applied to define the so-called adiabaticconnection of DFT, because one can apply the implicit functions theoremand let w decrease while the inverse potential v increases, keeping the densityfixed. 4. The dual problem when w = 0 In this section, we provide properties on the dual problem in the Kohn-Sham non-interacting case. The decomposition of the eigenfunctions intoSlater determinants will enable to get more information on the local problemand on the Euler-Langrange optimality conditions.4.1. Definitions. We first make some definitions. • We consider the exact continuous model (2) with potentials in L p + L ∞ and p as in (1), we saw in Section 2 that the discrete version of G ( k ) ρ iscoercive. We denote by D ( k ) N ( v ) := Ker R (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) (15)the real vector eigenspace of real eigenfunctions of the N -body operator,counted with multiplicity. • Let R max be the set of densities such that G (0) ρ has a maximizer. Wecan define the degeneracy of a density deg : R max −→ N \ { } ρ (cid:55)−→ dim R Ker R (cid:0) H N ( v ρ ) − E (0) ( v ρ ) (cid:1) , where v ρ is the unique maximizer of G (0) ρ . This map should have a richstructure. • In the case of w = 0 , that is when the model is an effective one-body one, we distinguish two types of degeneracies in the N -body problem.Let us denote by ( ϕ i ) i and ( E i ) i the real eigenfunctions and correspondingeigenvalues of − ∆ + v , where i (cid:55)→ E i is increasing. The eigenfunctions ofthe many-body problem (cid:80) Nj =1 ( − ∆ j + v ( x j )) are the antisymmetric products ∧ i ∈ I ϕ i , where I ⊂ I tot and | I | = N , having energies (cid:80) i ∈ I E i , hence D ( k ) N ( v ) = Span R (cid:8) ∧ i ∈ I ϕ i (cid:12)(cid:12) I ∈ I tot (cid:9) has a basis labeled by I tot . To a k th N -body energy level correspond sets I oforbitals. We say that the N -body degeneracy has a coincidental degeneracywhen several sums of energy levels “accidentally” superpose while at leastone energy level is different, that is when there are I, J ⊂ I tot such that (cid:80) i ∈ I E i = (cid:80) i ∈ J E i and there is m ∈ { , . . . , N } such that the m th (order byenergy) elements of I and J , denoted by I m and J m verify E I m (cid:54) = E J m . Weillustrate it on the left of Figure 1. When there exists j ∈ I tot \ I such that E j = E i , there is a one-body level which is partially occupied, the N -bodydegeneracy comes from a one-body degeneracy, and we say that the N -bodydegeneracy is essentially one-body , as illustrated in the middle of Figure1. We remark that coincidental and essentially one-body degeneracies cancoexist, as examplified on the right of Figure 1. Figure 1. Three kinds of N -body degeneracies: coinciden-tal, essentially one-body, and both. • Given a discrete set S = { s i } i ∈ N ⊂ R bounded below and such that s i (cid:54) s j for i < j , we define, for each j ∈ N ∪ { } , µ j ( S ) := s j . (16)For a n × n diagonalizable matrix M having real eigenvalues, we define µ j ( M ) := µ j ( D M ) where D M is the set of eigenvalues. • Consider that (Ψ tot I ) I ∈I tot is a basis of D ( k ) N ( v ) , where Ψ tot I = ∧ i ∈ I ϕ i , ϕ i being the orbitals of − ∆ + v as defined before. Let us define the “inner”orbitals I in := (cid:8) i ∈ N (cid:12)(cid:12) ∀ I ∈ I tot , i ∈ I (cid:9) , which are present in all the many-body functions of D ( k ) N ( v ) , they necessarily fill their energy levels. We definethe “inner” density ρ in := (cid:88) i ∈ N i ∈ I ∀ I ∈I tot ϕ i . We now drop these “inner” orbitals and only consider those which can changeon D ( k ) N ( v ) . Now we define the set I out := (cid:8) I \ I in (cid:12)(cid:12) I ∈ I tot (cid:9) of ( N − | I in | ) -body wavefunctions. The disjoint union of sets is denoted by ∪· . For k ∈ N ,we define I ( k ) := (cid:8) i ∈ I (cid:12)(cid:12) E v ( ϕ i ) = E k (cid:9) the subset of orbitals belonging to the k th one-body level. For I, J ∈ I out , we define the elements of the function-valued matrix of “one-body correlations” ( M ϕ ) IJ = δ IJ (cid:88) i ∈ I ϕ i + δ ∃ (cid:96),i,j ∈ N I ∪· J = I ( (cid:96) ) ∪· J ( (cid:96) ) = { i,j } I ( t ) = J ( t ) ∀ t (cid:54) = (cid:96) ϕ i ϕ j . (17)For instance when the only degeneracy is essentially one-body, and comesfrom a one-body level Span( ϕ (cid:96) + i ) (cid:54) i (cid:54) D with degeneracy D , filled with oneparticle, then M ϕ = (cid:0) ϕ (cid:96) +1 . . . ϕ (cid:96) + D (cid:1) T ⊗ (cid:0) ϕ (cid:96) +1 . . . ϕ (cid:96) + D (cid:1) . In thesame situation but when D = 3 and when there are two particles in theone-body-level, M ϕ = φ + φ φ φ φ φ φ φ φ + φ φ φ φ φ φ φ φ + φ , where φ i := ϕ (cid:96) + i .4.2. Local problem. By the next lemma, we can say that coincidentaldegeneracies “do not correlate” the many-body eigenstates in our problem. Lemma 4.1 (Local problem when w = 0 ) . Take a density ρ ∈ L (Ω) , ρ (cid:62) , (cid:82) ρ = N , √ ρ ∈ H (Ω) , consider a binding v ∈ V ( k ) N,∂ , and take w = 0 . UILDING INVERSE POTENTIALS 15 i ) We have + δ v G ( k ) ρ ( u ) = µ k − m k (cid:18)(cid:90) Ω u (cid:0) ( ρ in − ρ ) + M ϕ (cid:1)(cid:19) ii ) The matrix M ϕ = M ⊕ · · · ⊕ M T is block diagonal, where the blockscorrespond to the different essentially one-body degeneracies. iii ) In the case of only coincidental degeneracies, with Ψ , . . . , Ψ dim D ( k ) N ( v ) being an orthonormal basis of Ker (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) composed of Slaterdeterminants, then + δ v G ( k ) ρ ( u ) = µ k − m k (cid:18)(cid:90) u ( ρ Ψ i − ρ ) (cid:19) (cid:54) i (cid:54) dim D ( k ) N ( v ) . We recall that µ is defined in (16). The block-diagonalization shows thatthe degeneracies which complexify the problem are the essentially one-bodydegeneracies, not the coincidental ones. In iii ) , the normalized direction u ∗ maximizing + δ v G ( k ) ρ ( u ) is (cid:12)(cid:12)(cid:12)(cid:12) ρ Ψ i − ρ || ρ Ψ i − ρ || Lq (cid:12)(cid:12)(cid:12)(cid:12) q − sgn( ρ Ψ i − ρ ) , for the corresponding i .4.3. Optimality.Proposition 4.2 (Euler-Lagrange inequations when w = 0 ) . Take a density ρ ∈ L (Ω) , ρ (cid:62) , (cid:82) ρ = N , √ ρ ∈ H (Ω) , consider a binding v ∈ V ( k ) N,∂ , take w = 0 , and assume that v maximizes G ( k ) ρ . i ) We have, a.e in Ω , ρ in + µ k − m k ( M ϕ ) (cid:54) ρ (cid:54) ρ in + µ M k − k ( M ϕ ) . (18) ii ) If there is only one “outer” particle, that is N = 1 + (cid:82) ρ in , then k = m k and ρ in (cid:54) ρ (cid:54) ρ in + (cid:88) i ∈ G ϕ i , where G is such that ( ϕ i ) i ∈ G is a basis of the one-body degenerate level Ker R ( − ∆ + v − E (cid:96) ) producing the degeneracy of the N -body level. iii ) If Ker (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) has only coincidental degeneracies, then ithas a k th excited pure Slater state Ψ = ∧ Ni =1 ϕ i such that ρ Ψ = ρ . We remark that the larger k − m k (cid:62) , the more the Euler-Lagrange in-equations are constraining. We conjecture that m k = k in any case, similarlyas in ii ) .In particular, when k = 0 , there are only essentially one-body degeneraciesand i ) becomes ρ in + min σ ( M ϕ ) (cid:54) ρ (cid:54) ρ in + max σ ( M ϕ ) . Finally, we numerically find that when all the different partially filledone-body levels have dimensions 2, then there is a pure k th excited staterepresenting ρ at optimality. We think that the sets of mixed states densitiesand pure densities are equal in such configuration. Numerical application In this section, implement the dual problem and compute inverse poten-tials. The algorithm is presented for w = 0 , but its extension to any w isimmediate.5.1. Definition of the problem. For practical reasons, we do not use theframework of the discretized space of potentials introduced previously. In-stead we discretize the space on a finite grid ( Z ∩ [ − L, L ]) d , this correspondsto the model studied in [10]. In this case and for d = 1 for instance, thelaplacian and the potential multiplication are respectively the matrices − ∆ = − . . . − − − ... . . . , v = v ( x ) 0 . . . v ( x ) ... . . . . We implemented the algorithm in Julia, using the package DFTK [22] toefficiently compute the eigenfunctions of the one-body Schrödinger operator,this library using plane waves basis. We consider H w =0 N ( v ) := N (cid:88) i =1 ( − ∆ i + v ( x i )) . Given some density ρ belonging to the set (cid:8) ρ ∈ L (Ω) (cid:12)(cid:12) ρ (cid:62) , (cid:82) ρ = N, √ ρ ∈ H (cid:9) , (19)and called the target density, our goal is to find a potential v such that H w =0 N ( v ) has at least one k th excited state and such that there exists amixed state Γ v,k supported on Ker R (cid:0) H w =0 N ( v ) − E ( k ) ( v ) (cid:1) such that ρ Γ v,k = ρ. The existence is justified by [10] for k = 0 , and by our previous results for k (cid:62) . We recall that for k = 0 the searched potential is unique by theHohenberg-Kohn theorem [10, 23]. We would also like to know whether theset of v -representable pure state densities (cid:26) ρ Ψ ( k ) v (cid:12)(cid:12) v ∈ ( L p + L ∞ )(Ω) , Ψ ( k ) v ∈ Ker (cid:0) H w =0 N ( v ) − E ( k ) ( v ) (cid:1) , (cid:90) Ω N (cid:12)(cid:12)(cid:12) Ψ ( k ) v (cid:12)(cid:12)(cid:12) = 1 (cid:27) (20)is dense in the set of densities (19). We previously saw in Corollary 3.7 thatit is dense when d = 1 and that it must not when d = 3 due to Conjecture3.8.More explicitely, since w = 0 , we are led to study the one-body operator − ∆ + v . Following the notations defined in Section 4.1, its eigenvectorsare denoted by ϕ i and the energies by E i = (cid:82) |∇ ϕ i | + (cid:82) v | ϕ i | . The N -body k th excited states Ψ I = ∧ i ∈ I ϕ i of H w =0 N ( v ) are labeled by I ∈ I tot ,the corresponding density is ρ Ψ I = (cid:80) i ∈ I | ϕ i | and the energy is E ( k ) ( v ) = (cid:80) i ∈ I E i . UILDING INVERSE POTENTIALS 17 Algorithm. As we can see with Theorems 2.2 and 3.5, we need to max-imize G ( k ) ρ , which is a well-posed problem. Indeed, since in our simulationsthe system lives in a bounded set and since the space is discretized, G ( k ) ρ is coercive. Once we found its maximizers, knowing whether there exists apure state density equal to ρ boils down to computing (8).5.2.1. Starting point. We apply a gradient ascent algorithm on G ( k ) ρ to max-imize it. We start from Bohm’s potential v := ∆ √ ρ √ ρ , which produces exactly ρ when N = 1 and k = 0 .5.2.2. Ascent direction. We start by computing the first eigenfunctions ( ϕ i ) (cid:54) i (cid:54) M of − ∆ + v , where a (cid:62) N + k . The exact minimal needed value of a can belarger than N + k in case of degeneracies. Then we compute all the energyconfigurations (cid:80) i ∈ I E i for I ⊂ { , . . . , a } such that | I | = N , and we storethose energies in the increasing order. This gives us the N -body spectrumof H N ( v ) , E ( k ) ( v ) being the ( k + 1) th number of this list, and we deduce Ker R (cid:0) H w =0 N ( v ) − E ( k ) (cid:1) by taking the configurations having energies close to E ( k ) ( v ) as will be detailed later.A direction of stepest ascent is given by (6) and (7), where we will take p = q = 2 for practical purposes. However, we will not take a stepest ascentdirection, but one solvingsup u ∈ L || u || L =1 + δ v G ( m k ) ( u ) = min Γ ∈S (cid:0) D ( k ) N ( v ) (cid:1) Γ (cid:62) , Tr Γ=1 || ρ Γ − ρ || L . (21)The supremum is attained by u ∗ = ρ Γ ∗ − ρ || ρ Γ ∗ − ρ || L , where Γ ∗ is an optimizer of the right hand side of (21). This is justifiedbecause + δ v G ( m k ) (cid:54) + δ v G ( k ) so in this direction u ∗ , we still have (cid:54) + δ v G ( k ) ( u ∗ ) . This scheme should lead to a maximum because local maximas of G ( k ) areglobal by Theorem 2.2. Moreover, we experimentally remark that by usingthe direction given by (21) or the stepest ascent direction given by (6), having k (cid:54) = m k is very rare, and for almost converged potentials, we always have k = m k . In some situations, the min/max problem (6) is not necessarily lightto compute, it complexifies the implementation, only marginally acceleratesthe convergence, and is not convenient for implementing temperature, furtherjustifying the use of the direction (21).In case of degeneracies, which happens “most of the time”, if we choosea direction u = ( ρ Γ − ρ ) / || ρ Γ − ρ || L where Γ is a randomly choosed mixedstate of Ker R (cid:0) H w =0 N ( v ) − E ( k ) ( v ) (cid:1) , the algorithm starts to diverge. Henceoptimizing over directions is necessary. Temperature. In order to smooth the behavior of the algorithm andimprove its convergence, we introduce a “temperature” effect. We define theset of N -body Slater functions of v built on eigenfunctions F ( v ) := (cid:8) ∧ i ∈ I ϕ i (cid:12)(cid:12) | I | = N, I ⊂ N ∪ { } (cid:9) and consider the problem P ( v ) := min Γ ∈S ( F ( v ))Γ (cid:62) , Tr Γ=1 e (cid:0) E v (Γ) − E ( k )( v ) T (cid:1) | E v (Γ) − E ( k ) ( v ) | (cid:54) t (cid:90) Ω ( ρ Γ − ρ ) , (22)which is in practice computed by an ODA [5–7] algorithm. We define ρ n := ρ Γ ∗ where Γ ∗ is a minimizer of P ( v n ) and n is the n th iteration step. InAppendix 3, we provide the computations needed in the implementation ofthis part of the algorithm.There are two “temperatures” T and t . First, the cut-off | E v (Γ) − E ( k ) ( v ) | (cid:54) t considerably lowers the dimension of the optimization set F ( v ) , droppingconfigurations having energies too far from the relevant one. Then, thesmoothing factor e − (cid:0) E v (Γ) − E ( k )( v ) T (cid:1) enables to take into account many-bodystates which do not exactly have energy E ( k ) ( v ) but are close, addressingdegeneracies in a continuous way. The absence of this last factor raisesdivergence issues.Let t n , T n be the temperatures at step n . We take t = E v N + k − E v B ( N + k ) where E v i is the i th eigenenergy of the one-body operator − ∆ + v and B isa parameter. We take T n = t n /D for any n ∈ N , with D being a parameter,and progressively decrease them by choosing t n = α (cid:98) n/M (cid:99) t , where α < ,which we call the cooling factor, and M (cid:62) . We remarked that cooling“by steps” with the factor α (cid:98) n/M (cid:99) , rather than with α n/M , improves theconvergence. In practice, we find that D = B = 10 , M = 5 and α = 3 / aregood trade-offs.5.2.4. Line search optimization in the direction found. The previous proce-dure provided us an ascent direction, we now want to optimize the step inthis direction. We define v ( λ ) := v n + λ ρ n − ρ || ρ n − ρ || L , and take parameters µ > and ν > . If P ( v ( ν n )) (cid:62) P ( v (0)) , then we com-pute P ( v ( µ j ν n )) for increasing values of j ∈ N \ { } until P (cid:0) v ( µ j +1 ν n ) (cid:1) (cid:54) P ( v ( µ j ν n )) , the last value of j being denoted by j ∗ n . If P ( v ( ν n )) (cid:54) P ( v (0)) ,then we compute P ( v ( µ − j ν n )) for j ∈ N \ { } until P (cid:0) v ( µ − ( j +1) ν n ) (cid:1) (cid:54) P ( v ( µ − j ν n )) , this defines j ∗ n in this other case. Finally, the new potentialwill be v n +1 := v (cid:0) µ j ∗ n ν n (cid:1) . We “learn” the step size in the sense that ν n +1 = µ j ∗ n ν n if || ρ n +1 − ρ || L /N (cid:54) − , and ν n +1 = ν otherwise. UILDING INVERSE POTENTIALS 19 Convergence criterion. We consider that the algorithm converged when (cid:112) P ( v n ) /N (cid:54) (cid:15) and α (cid:98) n/M (cid:99) (cid:54) δ , where we take (cid:15) = 10 − and δ = 10 − . Thisensures that we found v and a mixed state Γ supported on the (cid:96) th excitedstates Ψ (cid:96)i such that || ρ Γ − ρ || L /N (cid:54) (cid:15) = 10 − , B ( N + k ) (cid:12)(cid:12) E v (Ψ (cid:96)i ) − E ( k ) ( v ) (cid:12)(cid:12) E v N + k − E v (cid:54) δ = 10 − . Pure states representability. Once we found an approximate maximiz-ing potential v n of G ( k ) ρ , if we want to know whether it produces a pure statedensity ρ , we computemin Ψ ∈ Span C F ( v n ) (cid:82) | Ψ | =1 | E vn (Ψ) − E ( k ) ( v n ) | (cid:54) t n (cid:90) (cid:0) ρ Ψ − ρ (cid:1) . (23)The above problem is computed using a particle swarm optimization algo-rithm implemented in the library Manopt.jl. The base manifold is a complex (dim D ( k ) N ( v ) − -dimensional Grassman manifold representing the optimiz-ing set of pure states.5.2.7. Remarks. • The global constant in potentials has no importance during the scheme,we only fix it in the graphs for readability purposes. • In the case where the k th level of H w =0 N ( v n ) is non-degenerate, onecan use a Newton or quasi-Newton algorithm like BFGS to accelerate theconvergence. With our notations and considering that the spectrum is purelydiscrete, for two directions u, h ∈ L p + L ∞ , d v G ( k ) ρ ( u, h ) = (cid:68) u, d v ρ ( k ) ( v ) h (cid:69) = 2 (cid:88) J ⊂ N , | J | = NI k ∪· J = { i,j } (cid:18) E ( k ) ( v ) − (cid:88) (cid:96) ∈ J E (cid:96) (cid:19) − (cid:104) ϕ i , uϕ j (cid:105) (cid:104) ϕ i , hϕ j (cid:105) , where v (cid:55)→ ρ ( k ) ( v ) is the potential-to-eigendensity map, and where I k is theconfiguration corresponding to Ker R (cid:0) H w =0 N ( v ) − E ( k ) ( v ) (cid:1) = R (cid:0) ∧ i ∈ I k ϕ i (cid:1) . Ina finite basis ( u i ) i , the Hessian ∇ v G ( k ) ρ = (cid:0) d v G ( k ) ρ ( u i , u j ) (cid:1) ij is non-degeneratewhen k = 0 . • The library DFTK is configured for periodic boundary conditions, butwe only consider densities which are very close to zero close to the boundaries,this implies potentials which are very large at the boundaries, and we recovera situation equivalent to a setting with Dirichlet boundary conditions. • As expected, for N = 1 and k = 0 , the potential v = ∆ √ ρ/ √ ρ hasa density very close to the target density ρ . For any d, N, k , the algorithmconverges significantly faster when we start from v = ∆ √ ρ/ √ ρ comparedto v = 0 .5.3. Convergence results. The convergence is theoretically justified byTheorem 3.5, and confirmed in our simulations. Up to slight adaptations ofthe parameters α , D , B , and M , the algorithm always converges both atthe levels of densities and of potentials, as expected, for any d ∈ { , , } Figure 2. Target densities and their Kohn-Sham potentialsfor k = 0 , On the first line, d = 1 and N = 3 , densities (blue)and inverse potentials (red) are in different units. On thesecond line, d = 2 , N = 2 , the target density is on the leftand its Kohn-Sham potential on the right. Figure 3. We take d = 3 , N = 4 , k = 1 , ρ is a sum of threeGaussians, and n = 50 . We plot v n and ρ first, and then ρ and log | ρ n − ρ | on the right.and any ρ, N, k . Moreover, the larger N , the faster the convergence. Weobtain arbitrary precision on (23) which scales as /n . We give a first qual-itative illustration on Figure 2, for k = 0 , related to the LDA which locallyapproximates densities by uniform electron gas partitions [19, 31, 32].As an illustration for d = 3 , we give on Figure 3 a representation of ρ , v n and log | ρ n − ρ | where N = 4 , k = 1 and where ρ is a sum of threeGaussians. UILDING INVERSE POTENTIALS 21 Figure 4. Simulations for d = 1 , N = 3 , and k = 0 on theleft, k = 1 in the middle and k = 5 on the right. Densitiesare in blue and inverse potentials are in other colors, theirdifferences come from different initial potentials v . Poten-tials are in the same units, but not densities, and n is largeenough so that || ρ n − ρ || L /N (cid:54) − .5.3.1. Uniqueness. On Figure 4, we represent inverse potentials in the case d = 1 , N = 3 , for the same density, varying k and varying v = ∆ √ ρ/ √ ρ + u i , where ( u i ) (cid:54) i (cid:54) is the same sequence for k ∈ { , , } , but the u i ’s arepairwise different. This confirms previous numerical studies [18] showingthat uniqueness of densities does not hold for k (cid:62) , in which the authorsused the kernel of d v (cid:0) u (cid:55)→ ρ Ψ ( k ) ( u ) (cid:1) to find such examples, where Ψ ( k ) ( u ) isthe non-degenerate k th eigenstate of a potential u . There must be an infinitenumber of inverse potentials, differing by “oscillations”.5.3.2. Reconstruction of potentials. We choose an initial potential v , com-pute ρ Γ v for some Γ v ∈ S (cid:0) D ( k ) N ( v ) (cid:1) with Γ v (cid:62) , Tr Γ v = 1 , and launch thealgorithm on ρ Γ v . For k = 0 , by the Hohenberg-Kohn theorem, v n shouldconverge to v , this is what we call the reconstruction of potentials.We show in Figure 5 an example of what one can obtain for d = 1 , in theground state setting k = 0 and for the third excited state k = 3 , with thesame target potential. The distance between densities is || ρ n − ρ || L (cid:39) − in the case k = 0 , and − for k = 3 .In the case d = 2 and on Figure 6 we show an example of reconstructionof a potential, where N = 5 and k = 0 . This corresponds to a precision || ρ n − ρ || L (cid:39) − for densities.On the graphs, the converged density ρ n is indistinguishable from thetarget ρ for d ∈ { , , } . Also, v n is indistinguishable from v for k = 0 in the regions where ρ ( x ) / (max ρ ) (cid:62) − , and v n is never close to v when k (cid:62) . Hence we do not plot ρ n and v n anymore. We can obtain arbitraryprecision on | v n − v | in the regions where ρ is not “small”, and as expected,we see that the precision on potentials is much lower than the precision ondensities, which is linked to the local weak-strong continuity of the map v (cid:55)→ ρ proved in [16].Moreover, when k = 0 , convergence is muh easier, and faster, because G (0) ρ is concave.5.4. Pure states, degeneracies and Levy-Lieb. We saw that our algo-rithm converges to some potential having a k th excited mixed state repre-senting ρ . Now we address the problem of finding a pure k th excited state Figure 5. Plot for d = 1 , N = 5 , in the ground state setting k = 0 on the left and in the third excited state k = 3 onthe right. We start from the same target potential v andlaunch the algorithm on the same ρ = ρ Ψ for some normalized Ψ ∈ Ker (cid:0) H w =0 N ( v ) − E ( k ) ( v ) (cid:1) , with n (cid:39) on the leftand n (cid:39) on the right. On the first line, densities andpotentials are in different units, and on the second line, weplot log | ρ n − ρ | (blue) and log | v n − v | (red) in the sameunits.representing ρ . The situation is very different as we choose d ∈ { , } or d = 3 .5.4.1. d = 1 . For one-dimensional systems, Theorem 2.2 iv ) justifies theexistence of pure states representing ρ .5.4.2. d = 2 . For two-dimensional systems, there is no non-degenerate theo-rem, and for instance the harmonic oscillator has arbitrary large degeneracyas we increase the energy level we consider. However as for d = 1 , we numer-ically remark that (22) equals (23) for any potential, and we always obtainpure-state representability, for any k .We do not have a theoretical proof of this fact, contrarily to d = 1 . Butroughly speaking, we think that this is due to the fact that the only relevantdegeneracies for our problem come from the spherical Laplacian. For d = 2 ,this operator is defined on the one-dimensional circle S and has only two-folddegeneracies, but the set of mixed states on a real vector space of dimension2 is equal to the set of pure states on this same vector space, as showed inTheorem 2.2 iii ) . Then, for instance degeneracies arising in the harmonicoscillator at the one-body level should be unessential degeneracies for ourproblem. UILDING INVERSE POTENTIALS 23 Figure 6. We initially choose d = 2 , N = 5 , k = 0 , andthe confining and non-degenerate potential v indicated thethe top left. We compute its density ρ ( v ) , displayed onthe top right, and launch the algorithm on it. We display log | v n − v | on the bottom left and log | ρ n − ρ ( v ) | on thebottom right, in their natural units, for n (cid:39) . The lengthof the squares is 4 in in space units.5.4.3. d = 3 . In the proof of [34, Theorem 3.4], Lieb identified a class of den-sities ρ such that F (0) mix ( ρ ) < F (0) ( ρ ) , using radial symmetry and degeneraciesof spherical orbitals in dimension . Since we think that those functionalsare continuous (see Conjecture 3.8), this indicates that some densities arenot (approximately) pure-state representable. Numerically, we verify it onLieb’s example, for ρ ( x ) = (2 πσ ) − / e −| x | / (2 σ ) , σ = 1 / , k = 0 , N = 2 ,min Ψ ∈ Span C F ( v n ) (cid:82) | Ψ | =1 | E vn (Ψ) − E ( k ) ( v n ) | (cid:54) t n || ρ Ψ − ρ || L (cid:38) − . , for any n ∈ N , while the distance with mixed states (22) is arbitrarily smallas n → + ∞ , confirming that the algorithm converges. Since the potentialto which we converge is unique, we conclude that there is no pure groundstate representing this ρ . Of course, we can also find infinitely many densities which are v -representablewith pure states, but a general necessary and sufficient condition for a densityto be v -representable with pure states seems out of reach.For k (cid:62) and d = 3 , we did not solve the complete problem of pure-staterepresentability (8), because we would have to test all maximizers, this set ofmaximizer must be very large and finding it should not be easy. We howeverthink that the conclusion would be the same as for k = 0 .5.5. Conclusions. • To obtain the mixed Kohn-Sahm potential, one can maximize G ( k ) ρ , inall situations. The algorithm presented here is simple and always converges. • For k (cid:62) , starting from a different potential or changing some param-eters of the algorithm leads to very different inverse potentials, emphasizingthat many potentials lead to the same density. • For d ∈ { , } , inverse potentials with pure states exist, hence thissuggests that F ( k ) mix ( ρ ) = F ( k ) ( ρ ) and that the set of v -representable densitieswith pure states (20) is dense in the space of densities. Based on the natureof the degeneracy (essentially one-body, coming from the spherical Laplaceoperator, ...), and on its number, for d = 3 and k = 0 , some densities do nothave k th pure states representing them and this indicates that (20) is notdense in the set of densities. However, Section 3.3.2, suggests that given ρ ,there exists k such that ρ is approximately a k th excited pure state density. • Degeneracies have to be taken into account in the algorithm, not onlywhen the inverse potential is non-degenerate, but all along the procedure,because intermediate degenerate potentials can block the algorithm, corre-sponding to eigenvalues crossings. The only case where we do not need totake them into account is when ( d, k ) = (1 , , because then the N -bodyground level is non-degenerate by the non-degeneracy theorem (Proposition6.2). Using our algorithm, we launched simulations on hundreds of densi-ties, choosed to be sums of gaussians with random parameters, and weresurprised by the fact that for d ∈ { , , } , “most of the time”, the inversepotential is degenerate. Perturbation of potentials generically lifts degenera-cies, but at the level of this inverse problem, perturbation of a density doesnot lift degeneracies. Hence, the fact that “generically” inverse potentials aredegenerate can seem counterintuitive. In the SCF procedure, perturbationdoes not lift degeneracies either, as shown in [8]. Finally, we remark that wealways have k = m k for v n close to optimality.6. Proofs Proofs of Proposition 2.1 and Theorem 2.2. We start by present-ing a remark. Let Q be a real vector space of real wavefunctions, then for Γ ∈ S N mix ( Q, Ω) , we have ρ Γ = ρ Re Γ . Indeed, self-adjointness Γ = Γ ∗ implies Γ( x , . . . , x N ; x , . . . , x N ) ∈ R . Incase dim Q is finite, take an orthonormal frame (Ψ i ) i ∈ I of real wavefunctionsof Q , write the hermitian matrix mat Q Γ =: (Γ ij ) (cid:54) i,j (cid:54) n , then Γ ii ∈ R and UILDING INVERSE POTENTIALS 25 we have ρ Γ = n (cid:88) i =1 Γ ii ρ Ψ i + 2 N (cid:88) (cid:54) i As explained in [16, Theorem 1.6], we have + δ v G ( k ) ρ ( u ) = min Q ⊂D ( k ) N ( v )dim Q = k − m k +1 max Ψ ∈ Span C Q (cid:82) | Ψ | =1 (cid:90) ( ρ Ψ − ρ ) u (24) = max Q ⊂D ( k ) N ( v )dim Q = M k − k min Ψ ∈ Span C Q (cid:82) | Ψ | =1 (cid:90) ( ρ Ψ − ρ ) u, where D ( k ) N ( v ) is defined in (15) and Span C Q is the complex vector spacebuilt on the real vector space Q . First, since ρ vanishes at infinity and D ( k ) N ( v ) is compact and formed by functions vanishing at infinity,sup || u || Lp + L ∞ =1 + δ v G ( k ) ρ ( u ) = sup || u || Lp =1 + δ v G ( k ) ρ ( u ) . Then, we will need the following classical lemma. Lemma 6.1. For any real linear subspaces A and B of H a ( R dN , R ) , where B is finite dimensional, and any potential v ∈ L p + L ∞ , we have inf Γ ∈S ( A )Γ (cid:62) , Tr Γ=1 E v (Γ) = inf Ψ ∈ Span C A (cid:82) | Ψ | =1 E v (Ψ) , max Γ ∈S ( B )Γ (cid:62) , Tr Γ=1 E v (Γ) = max Ψ ∈ Span C B (cid:82) | Ψ | =1 E v (Ψ) . We use the notation S ( A ) in case of infinite-dimensional vectors spaces A ,by the natural extension of the definition in finite dimension. The statementof Lemma 6.1 can be seen by taking a real mixed state Γ ∈ S N mix ( A, Ω) ,decomposing it into Γ = (cid:80) + ∞ i =1 λ i | Ψ i (cid:105) (cid:104) Ψ i | where Ψ i ∈ A is an orthonormalbasis of real wavefunctions of A , λ i ∈ R + and (cid:80) + ∞ i =1 λ i = 1 . We evaluate E v (Γ) = (cid:80) + ∞ i =1 λ i E v (Ψ i ) (cid:62) inf Ψ ∈ A, (cid:82) | Ψ | =1 E v (Ψ) , because E v (Ψ) = E v ( α Ψ) for any α ∈ C . This is similar for the max.Now we use the theorem of existence of saddle points [44, Theorem 49.Ap459] for the Lagrangian L (Γ , u ) := (cid:82) u ( ρ Γ − ρ ) , affine in its variables, tocommutesup || u || Lp (cid:54) min Γ ∈S ( A )Γ (cid:62) , Tr Γ=1 (cid:90) u ( ρ Γ − ρ ) = min Γ ∈S ( A )Γ (cid:62) , Tr Γ=1 sup || u || Lp (cid:54) (cid:90) u ( ρ Γ − ρ ) , where A is any finite-dimensional real vector space. But since the twosuprema are positive, this also holds when we optimize over {|| u || L p = 1 } .Next, for f ∈ L q , sup || g || Lp =1 (cid:90) f g = λ (cid:90) f pp − = || f || L q . Indeed, considering the Lagrangian L ( f, λ ) = (cid:82) f g − λ (cid:0)(cid:82) g p − (cid:1) , and search-ing for its extremizer, we obtain f = λp | g | p − sgn g . The condition || g || L p = 1 yields λp = || f || L q , and the extremizer is g = (cid:12)(cid:12)(cid:12)(cid:12) f || f || L q (cid:12)(cid:12)(cid:12)(cid:12) q − sgn f. Hence for any (cid:54) Γ ∈ S ( A ) such that Tr Γ = 1 ,sup || u || Lp =1 (cid:90) u ( ρ Γ − ρ ) = || ρ Γ − ρ || L q is attained by (7).Using all the previous steps, we deduce thatsup || u || Lp =1 + δ v G ( k ) ( u ) = sup || u || Lp =1 max Q ⊂D ( k ) ( v )dim Q = k − m k min Ψ ∈ Span C ( Q ⊥ ∩D ( k ) ( v ) ) (cid:82) | Ψ | =1 (cid:90) u ( ρ Ψ − ρ )= sup || u || Lp =1 max Q ⊂D ( k ) ( v )dim Q = k − m k min Γ ∈S ( Q ⊥ ∩D ( k ) ( v ) ) Γ (cid:62) , Tr Γ=1 (cid:90) u ( ρ Γ − ρ )= max Q ⊂D ( k ) ( v )dim Q = k − m k sup || u || Lp =1 min Γ ∈S ( Q ⊥ ∩D ( k ) ( v ) ) Γ (cid:62) , Tr Γ=1 (cid:90) u ( ρ Γ − ρ )= max Q ⊂D ( k ) ( v )dim Q = k − m k min Γ ∈S ( Q ⊥ ∩D ( k ) ( v ) ) Γ (cid:62) , Tr Γ=1 sup || u || Lp =1 (cid:90) u ( ρ Γ − ρ )= max Q ⊂D ( k ) ( v )dim Q = k − m k min Γ ∈S ( Q ⊥ ∩D ( k ) ( v ) ) Γ (cid:62) , Tr Γ=1 || ρ Γ − ρ || L q = max Q ⊂D ( k ) ( v )dim Q = M k − k +1 min Γ ∈S ( Q )Γ (cid:62) , Tr Γ=1 || ρ Γ − ρ || L q , where Q are real vector spaces. (cid:3) Proof of Proposition 2.2. • i ) a ) = ⇒ c ) The state Γ is a k th mixed excited state of v and hasdensity ρ Γ = ρ , hence in the energy we can restrict the optimization searchto states having density ρ , E v (Γ) = E ( k ) ( v ) = sup A ⊂ H a ( R dN )dim A = k +1 inf Λ ∈S N mix ( A ⊥ , Ω ) E v (Λ)= sup A ⊂ H a ( R dN )dim A = k +1 inf Λ ∈S N mix ( A ⊥ , Ω ) ρ Λ = ρ E v (Λ) = F ( k ) mix ( ρ ) + (cid:90) vρ = sup u G ( k ) ρ ( u ) + (cid:90) vρ. On the other hand, G ( k ) ρ ( v ) = E v (Γ) − (cid:82) vρ , thus v maximizes G ( k ) ρ . b ) = ⇒ a ) Since v is a local maximizer and V ( k ) N,∂ is open, then for u close enough to v , we have G ( k ) ρ ( u ) (cid:54) G ( k ) ρ ( v ) and u ∈ V ( k ) N,∂ . Thus sup u + δ v G ( k ) ρ ( u ) (cid:54) , and we conclude by (6) that for any Q ⊂ Ker R (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) , with dim Q = M k − k + 1 , there exists a (cid:54) Γ = Γ T ∈ S ( Q ) suchthat Tr Γ = 1 and ρ Γ = ρ . UILDING INVERSE POTENTIALS 27 Remark when k = 0 . Let v be an optimizer of G (0) ρ . We know that thereexists a mixed state Γ ρ such that ρ Γ ρ = ρ and Tr H N (0)Γ = G (0) ρ ( v ) by [34,Corollary 4.5]. This implies that Tr H N ( v )Γ ρ = E (0) ( v ) . By diagonalizing Γ ρ similarly as in the proof of Theorem 3.5, we can show that it is a groundmixed state for v . • ii ) We also know that + δ v G ( m ) ρ (cid:54) + δ v G ( n ) ρ for any m, n ∈ { m k , . . . , M k } such that m (cid:54) n . Hence for (cid:96) ∈ { m k , . . . , k } , sup u + δ v G ( (cid:96) ) ρ ( u ) (cid:54) sup u + δ v G ( k ) ρ ( u ) = 0 and v is a local maximizer of G ( (cid:96) ) ρ .We can see from (24) that + δ v G ( k ) ρ ( − u ) = − + δ v G ( M k + m k − k ) ρ ( u ) for any direction u ∈ L p + L ∞ . Thus − + δ v G ( k ) ρ ( − u ) = + δ v G ( M k + m k − k ) ρ ( u ) (cid:54) + δ v G ( k ) ρ ( u ) (cid:54) , where we used that M k + m k − k (cid:54) k and that v is a local maximizer. Wededuce that for any direction u , + δ v G ( k ) ρ ( u ) (cid:62) , so + δ v G ( k ) ρ = d v G ( k ) ρ = 0 .This yields (cid:0) d v E ( k ) (cid:1) u = (cid:82) uρ and we can conclude.If v is a local minimizer and k (cid:62) ( m k + M k ) / , by a similar reasoning wehave d v G ( k ) ρ = 0 , and by (6) there is a k th excited mixed state of v such that ρ Γ = ρ . As we saw in i ) , this implies that v is a global maximizer, but itcannot be a local minimizer at the same time. • iii ) Take real orthonormal vectors ϕ, φ we write Γ = (cid:18) a cc b (cid:19) ∈ R × in this basis. The condition Γ (cid:62) is equivalent to a (cid:62) , b (cid:62) , ab (cid:62) c , and Tr Γ = a + b = 1 . Then ρ Γ = aρ ϕ + bρ φ + 2 N c (cid:90) R d ( N − ϕφ. In the case of pure states, we have ψ = αϕ + βφ , with | α | + | β | = 1 so wecan take the parametrization α = (cos t ) e iη , β = (sin t ) e i ( η + θ ) , and ρ ψ = (cos t ) ρ ϕ + (sin t ) ρ φ + 2 (cos t ) (sin t ) (cos θ ) N (cid:90) R d ( N − ϕφ. Since (cid:8)(cid:0) (cos t ) , (sin t ) , cos t sin t cos θ (cid:1) (cid:12)(cid:12) t, θ ∈ [0 , π ] (cid:9) = (cid:8) ( a, b, c ) (cid:12)(cid:12) a, b, c ∈ R , c (cid:54) ab, a + b = 1 (cid:9) , the two spanned spaces of density are equal (cid:110) ρ ψ (cid:12)(cid:12) ψ ∈ Span C ( ϕ, φ ) , (cid:82) | ψ | = 1 (cid:111) = (cid:8) ρ Γ (cid:12)(cid:12) Γ ∈ S (Span R ( ϕ, φ )) , Γ (cid:62) , Tr Γ = 1 (cid:9) . Since we know that v has a k th excited mixed state representing ρ , it alsohas a pure one. • iv ) We will use a well-known result specific to the dimension one. Proposition 6.2 (Non-degeneracy theorem) . For d = 1 and any potential v ∈ ( L + L ∞ )( R ) , every eigenstate of − ∆ + v is non-degenerate. We recall its proof for the convenience of the reader. Proof. Let ψ, ϕ ∈ H ( R , R ) be such that ( − d / d x + v − E ) φ = 0 for φ ∈{ ψ, ϕ } , with E ∈ R . Multiplying the first equation by ϕ , the second by ψ andsubstracting, we get ψϕ (cid:48)(cid:48) − ϕψ (cid:48)(cid:48) = ( ψϕ (cid:48) − ϕψ (cid:48) ) (cid:48) . Hence ψϕ (cid:48) − ϕψ (cid:48) = c for some constant c ∈ R , but since ψ, ϕ ∈ L , then c = 0 . We have thus ψ (cid:48) /ψ = ϕ (cid:48) /ϕ on (cid:8) x ∈ R (cid:12)(cid:12) ψ (cid:54) = 0 , ϕ (cid:54) = 0 (cid:9) , which has full measure by uniquecontinuation [17]. Finally ψ = aϕ with a = ± . (cid:3) This result shows that the eigenspaces of H N ( v ) with w = 0 can only havecoincidental degeneracies, and we apply Proposition 4.2 iii ) . (cid:3) Proofs of Theorems 3.1 and 3.2. We start by proving the existenceof minimizers. Proof of Theorem 3.1. Let us denote by Ψ n a minimizing sequence for F α , (0) ( r ) .Since w (cid:62) , then (cid:90) (cid:12)(cid:12) ∇√ ρ Ψ n (cid:12)(cid:12) (cid:54) (cid:90) |∇ Ψ n | (cid:54) E (Ψ n ) → F α , (0) ( r ) , by the Hoffman-Ostenhof inequality, and Ψ n is bounded in H (Ω N ) and thereexists Ψ ∞ ∈ H a (Ω) such that Ψ n (cid:42) Ψ ∞ in H (Ω N ) and √ ρ Ψ n (cid:42) √ ρ ∞ in H (Ω) . At this step, ρ ∞ and Ψ ∞ are not related. By summing allthe constraints on the density and using that Ω (cid:80) i ∈ I α i = Ω , we have (cid:82) ρ Ψ n = N . We estimate (cid:90) B c r ∩ Ω ρ Ψ n (cid:54) (cid:90) ρ Ψ n (cid:88) supp α i ∩ B c r (cid:54) = ∅ α i = (cid:88) supp α i ∩ B c r (cid:54) = ∅ r i , and using the assumption (9) yields sup n (cid:82) B c r ∩ Ω ρ Ψ n → when r → + ∞ .This implies that √ ρ Ψ n converges strongly in L (Ω) and weakly in H (Ω) ,up to extraction of a subsequence. The tightness of ρ Ψ n also implies that Ψ n → Ψ ∞ strongly in L (Ω N ) , that the limit of ρ Ψ n is ρ Ψ ∞ , and eventuallythat (cid:82) ρ Ψ ∞ α i = r i . By lower semi-continuity of the energy functional since w (cid:62) , E (Ψ ∞ ) (cid:54) F α , (0) ( r ) , hence Ψ ∞ is a minimizer. By equivalence of thequadratic form E with the one of H , Ψ → Ψ ∞ strongly in H (Ω) .In the mixed state case, let us denote by Γ n a minimizing sequence. Weuse the compactness of the Fock space of particle number less than N , S N mix (cid:0) F (cid:54) N (cid:1) for the geometric convergence [30, Lemma 2.2], we thus have Γ n (cid:42) g Γ ∞ for some Γ ∞ ∈ S N mix (cid:0) F (cid:54) N (cid:1) . As before, the tightness of ρ Γ n implies Γ n → Γ ∞ strongly in trace-class by [30, Lemma 2.3], hence Γ ∞ is an N -particle density matrix. (cid:3) One does not need the weight functions α i to have the diameters of theirsupports converging to zero to get that our regularized functionals convergeto the exact one. Nevertheless, when this is the case, we can deduce boundson the rate of convergence of the densities of minimizers to the target den-sity. More precisely, the following result quantifies the distance between two UILDING INVERSE POTENTIALS 29 densities satisfying (cid:82) ρα i = (cid:82) χα i for all i ∈ I . We consider exponents q ∈ [1 , + ∞ ) if d = 1 , q ∈ [1 , if d = 2 and q = d/ ( d − if d (cid:62) . (25) Lemma 6.3 (Bounds on approximate densities) . Let q be as in (25) . Let Ω be an open set, α i = β k i i where k i ∈ R + and β i are non-negative concavefunctions on supp α i , with sup i ∈ N diam supp α i < + ∞ , Ω (cid:80) + ∞ i =1 α i = Ω .Let ρ, χ ∈ L (Ω , R + ) such that √ ρ, √ χ ∈ H (Ω) . If (cid:82) ρα i = (cid:82) χα i for any i ∈ I , then || ρ − χ || ( L ∩ L q )(Ω) (cid:54) c d (cid:0) ||√ ρ || H + ||√ χ || H (cid:1) sup i ∈ N diam supp α i , where c d only depends on d .Proof of Lemma 6.3. Take p (cid:62) . We use the weighted Poincaré-Wirtingerinequality from [11, Theorem 1.1] with c f,i := (cid:0)(cid:82) α i (cid:1) − (cid:82) f α i . We obtain (cid:90) | f − c f,i | p α i (cid:54) c p (diam supp α i ) p (cid:90) |∇ f | p α i , for f ∈ { ρ, χ } . Thus, since c ρ,i = c χ,i by assumption, (cid:90) | ρ − χ | p α i = (cid:90) | ρ − c ρ,i − ( χ − c ρ,i ) | p α i (cid:54) c p (cid:18)(cid:90) | ρ − c ρ,i | p α i + (cid:90) | χ − c ρ,i | p α i (cid:19) (cid:54) c p (cid:18) sup i ∈ N diam supp α i (cid:19) p (cid:90) (cid:0) |∇ ρ | p + |∇ χ | p (cid:1) α i . Summing over i and raising to the power /p yields || ρ − χ || L p (cid:54) c p (cid:18) sup i ∈ N diam supp α i (cid:19) (cid:0) ||∇ ρ || L p + ||∇ χ || L p (cid:1) . (26)We take p = q , so by writing ∇ f = 2 √ f ∇√ f and by the Sobolev inequality,we have ||∇ f || L p (cid:54) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∇ (cid:112) f (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) L (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:112) f (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) L p − p (cid:54) c d (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∇ (cid:112) f (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) L . Applying (26) concludes the proof. (cid:3) We prove now the convergence of our regularized functionals to the exactones. Proof of Theorem 3.2. Let us denote by Ψ n a sequence of approximate min-imizers for F α n , (0) ( r ρ ) . Since w (cid:62) , then (cid:90) (cid:12)(cid:12) ∇√ ρ Ψ n (cid:12)(cid:12) (cid:54) (cid:90) |∇ Ψ n | (cid:54) E (Ψ n ) (cid:54) F α n , (0) ( r ρ ) + (cid:15) n (cid:54) F (0) ( ρ ) + (cid:15) n , where (cid:15) n → . Hence (Ψ n ) n ∈ N is bounded in H (Ω N ) and there exists Ψ ∞ ∈ ∧ N H (Ω) such that Ψ n (cid:42) Ψ ∞ weakly in H (Ω N ) . By summing allthe constraints on the density and using that Ω (cid:80) i ∈ I α i = Ω , we have (cid:82) ρ Ψ n = N . We estimate (cid:90) B c r ∩ Ω ρ Ψ n (cid:54) (cid:90) ρ Ψ n (cid:88) supp α ni ∩ B c r (cid:54) = ∅ α ni = (cid:90) ρ (cid:88) supp α ni ∩ B c r (cid:54) = ∅ α ni , and using the assumption (11) yields sup n (cid:82) B c r ∩ Ω ρ Ψ n → when r → + ∞ .This implies that √ ρ Ψ n converges strongly in L (Ω) and weakly in H (Ω) ,up to extraction of a subsequence. The tightness of ρ Ψ n also implies that Ψ n → Ψ ∞ strongly in L (Ω N ) and that the limit of ρ Ψ n is ρ Ψ ∞ .Let f ∈ C ∞ c (Ω) , by assumption (10), there exists a sequence of functions f n ∈ Span { α ni , i ∈ I n } such that || f − f n || L p + L ∞ → when n → + ∞ . Wealso have (cid:82) f n ρ Ψ n = (cid:82) f n ρ because f n ∈ Span { α ni , i ∈ I n } . By using (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) f ( ρ Ψ n − ρ ) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) ( f − f n )( ρ Ψ n − ρ ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:54) c d (cid:18) ||√ ρ || H + sup n ∈ N (cid:12)(cid:12)(cid:12)(cid:12) √ ρ Ψ n (cid:12)(cid:12)(cid:12)(cid:12) H (cid:19) || f − f n || L p + L ∞ , we deduce that (cid:82) f ρ Ψ n → (cid:82) f ρ . This is a convergence of ρ Ψ n to ρ in thesense of distributions and by uniqueness of the limit, we have then ρ = ρ Ψ ∞ .We deduce that Ψ ∞ belongs to the minimizing set of F (0) ( ρ ) , consequently F (0) ( ρ ) (cid:54) E (Ψ ∞ ) . By also using lower semi-continuity of the energy func-tional since w (cid:62) , we have F (0) ( ρ ) (cid:54) E (Ψ ∞ ) (cid:54) lim inf E (Ψ n ) (cid:54) F (0) ( ρ ) . We have thus equality and we conclude that Ψ ∞ is a minimizer of F ( ρ ) . Letus consider the quadratic form q (Ψ) := (cid:104) Ψ , H N (0)Ψ (cid:105) . The convergence onthe Levy-Lieb functionals F α n , (0) ( r ρ ) → F (0) ( ρ ) gives q (Ψ n ) → q (Ψ ∞ ) , andsince w (cid:62) is ( − ∆) -bounded as a quadratic form, the associated norm of q is equivalent to the H norm, and hence Ψ n → Ψ ∞ in H .In the mixed states case, we follow a similar adaptation as for provingTheorem 3.1. As in the pure states case, the norm of E is equivalent to thenorm of S , , hence Γ n → Γ ∞ strongly in S , . (cid:3) The dual problem: proof of Theorem 3.5. In this section, we proveTheorem 3.5 on the coercivity of the dual functional G ( k ) r, α . In the proofs wewill use the notation E ( k ) dis ( v ) := E ( k ) (cid:32)(cid:88) i ∈ I v i α i (cid:33) , V ( v ) := (cid:88) i ∈ I v i α i . We recall that c Ω := − E ( k ) (0) /N is the constant such that the energy E ( k ) (cid:0) (cid:80) i ∈ I c Ω α i (cid:1) = E ( k ) ( c Ω ) = 0 vanishes. We present a fact about the signof the potential. Lemma 6.4. Let v ∈ (cid:96) ∞ ( I, R ) be such that E (0) dis ( v ) = 0 and v (cid:54) = c Ω . If H N (cid:0) V ( v ) (cid:1) has a ground state, then there exists i ∈ I such that v i < c Ω . If Ω is bounded, there exist i, j ∈ I such that v i < c Ω < v j .Proof of Lemma 6.4. UILDING INVERSE POTENTIALS 31 • Let Ψ be a ground state of H N ( V ( v )) . We have E (0) dis ( v ) = E V ( v ) (Ψ) = E c Ω (Ψ) + (cid:90) ( V ( v ) − c Ω ) ρ Ψ (cid:62) E (0) (cid:0) c Ω (cid:1) + (cid:90) ( V ( v ) − c Ω ) ρ Ψ = (cid:88) i ∈ I ( v i − c Ω ) (cid:90) ρ Ψ α i . Thus (cid:88) v i >c Ω | v i − c Ω | (cid:90) ρ Ψ v α i (cid:54) (cid:88) v i Proof of Theorem 3.5. • We first prove (13). We assumed that there are points y i ∈ R d suchthat for any i ∈ I , B R ( y i ) ⊂ (supp α i ) \ ∪ j ∈ I,j (cid:54) = i supp α j . We write X = ( x , . . . , x N ) and Y i = ( y i , . . . , y i ) . Take normalized Φ , . . . , Φ k ∈∧ N H ( B R ) with disjoint supports. Take some non-empty Q ⊂ I and for j ∈ { , . . . , l } , form Ψ j,Q := 1 (cid:113)(cid:80) i ∈ Q r i (cid:88) i ∈ Q √ r i Φ j ( X − Y i ) . This satisfies (cid:82) R dN | Ψ j,Q | = 1 , T (Ψ j,Q ) = T (Φ j ) , W (Ψ j,Q ) = W (Φ j ) and ρ Ψ j,Q ( x ) = (cid:88) i ∈ Q r i − (cid:88) i ∈ Q r i ρ Φ j ( x − y i ) . We use the expression E ( k ) ( V ) = inf dim A = k +1 max Ψ ∈ A (cid:82) | Ψ | =1 E V (Ψ) and choose the frame A := (cid:0) Ψ ,Q , . . . , Ψ k,Q (cid:1) so that G ( k ) r, α ( v ) (cid:54) − (cid:88) i ∈ I v i r i + max λ j ∈ C (cid:80) kj =0 | λ j | =1 E V ( v ) k (cid:88) j =0 λ j Ψ j,Q . For any i ∈ I , the only non-vanishing element of α in B R ( y i ) is α i , so α i = 1 on B R ( y i ) and (cid:90) Ω α i ρ Ψ j,Q = N r i δ i ∈ Q (cid:80) (cid:96) ∈ Q r (cid:96) , (cid:90) Ω V ( v ) ρ (cid:80) kj =0 λ j Ψ j,Q = k (cid:88) j =0 | λ j | (cid:90) Ω V ( v ) ρ Ψ j,Q = N (cid:80) i ∈ Q v i r i (cid:80) (cid:96) ∈ Q r (cid:96) . We see that the external potential energy of the trial state does not dependon the λ j ’s. Defining c R := max λ j ∈ C (cid:80) kj =0 | λ j | =1 E k (cid:88) j =0 λ j Ψ j,Q , we deduce that G ( k ) r, α (cid:0) v (cid:1) (cid:54) c R + N (cid:80) i ∈ Q r i (cid:88) i ∈ Q v i r i − (cid:88) i ∈ I v i r i (29) = c R + (cid:80) i ∈ I \ Q r i (cid:80) i ∈ Q r i (cid:88) i ∈ Q v i r i − (cid:88) i ∈ I \ Q v i r i . Since G is gauge invariant, for any µ ∈ R and any non-empty Q ⊂ I , wehave G ( k ) r, α ( v ) = G ( k ) r, α ( v − µ ) (cid:54) c R + (cid:80) I \ Q r i (cid:80) i ∈ Q r i (cid:88) i ∈ Q ( v i − µ ) r i − (cid:88) i ∈ I \ Q ( v i − µ ) r i . (30)We define the two sets I ± v := (cid:8) i ∈ I (cid:12)(cid:12) ± v i > c Ω (cid:9) . In the case I − v (cid:54) = ∅ , wetake Q = I − v and µ = c Ω yielding G ( k ) r, α ( v ) − c R (cid:54) (cid:80) v i (cid:62) c Ω r i (cid:80) v i In the case v = c Ω , we have G ( k ) r, α ( v ) = − c Ω N = E ( k ) (0) (cid:54) c R by using thesame trial state as before, hence the bound also holds.If v (cid:62) c Ω , then inf v = c Ω , because otherwise inf v > c Ω and E ( k ) dis ( v ) (cid:62) E ( k ) dis (inf v ) > E ( k ) dis ( c Ω ) = 0 . We take a sequence ( v ϕ ( n ) ) n ∈ N where ϕ ( n ) ∈ I is such that v ϕ ( n ) → c Ω when n → + ∞ . If I is finite, then there is (cid:96) ∈ I such that v (cid:96) = 0 and we take ϕ ( n ) = (cid:96) for any n ∈ N . We choose Q = { ϕ ( n ) } with only one element, and(29) yields G ( k ) r, α ( v ) (cid:54) c R + N v ϕ ( n ) − (cid:88) i ∈ I v i r i . By gauge invariance, we also have G ( k ) r, α ( v ) = G ( k ) r, α ( v − c Ω ) (cid:54) − || v − c Ω || (cid:96) r + N (cid:0) v ϕ ( n ) − c Ω (cid:1) + c R , where we used v (cid:62) c Ω . Taking the limit n → + ∞ yields G ( k ) r, α ( v ) (cid:54) − || v − c Ω || (cid:96) r + c R . In case k = 0 and I is infinite, we can still prove that there is a maximizer.Let v n be a maximizing sequence. By coercivity, (cid:80) i ∈ I | v ni | r i is boundedhence | v ni | (cid:54) c/r i uniformly in i, n . There exists ( v ∞ i ) i ∈ I ∈ (cid:96) ∞ ( I, R ) suchthat v ni → v ∞ i for all i ∈ I , up to a subsequence and finally v ∞ ∈ (cid:96) r ( I, R ) byFatou’s lemma. We conclude by using weak upper semi-continuity of G (0) r, α .For k (cid:62) , G ( k ) r, α is not upper-semi continuous or concave but when I is finite,since it is coercive and lives in a finite-dimensional space, it has a maximum. • Now assume that Ω is bounded, so that every potential is binding,moreover I is necessarily finite. We first give the beginning of a prove whichwould not use Theorem 2.2 to see why it only works for k = 0 .We would first need that F α , ( k ) mix ( r ) has an optimizer (Γ r , A r ) , this wasproved by Lieb [34] for k = 0 and we are not able to prove it for k (cid:62) . Since (cid:82) ρ Γ r α i = r i , then E V ( v ∞ ) (Γ r ) = Tr H N (0)Γ r + (cid:88) i ∈ I v ∞ i r i = G ( k ) r, α ( v ∞ ) + (cid:88) i ∈ I v ∞ i r i = E ( k ) dis ( v ∞ ) . We diagonalize Γ r =: (cid:80) (cid:96) ∈ N λ (cid:96) | ϕ (cid:96) (cid:105) (cid:104) ϕ (cid:96) | , where (cid:80) (cid:96) ∈ N λ (cid:96) = 1 and ϕ (cid:96) ∈ A ⊥ r .By linearity of the energy functional, we have (cid:88) (cid:96) ∈ N λ (cid:96) E V ( v ∞ ) ( ϕ (cid:96) ) = E ( k ) dis ( v ∞ ) . We need again k = 0 , because then E V ( v ∞ ) ( ϕ (cid:96) ) (cid:62) E (0) dis ( v ∞ ) and thus E V ( v ∞ ) ( ϕ (cid:96) ) = E (0) dis ( v ∞ ) for any (cid:96) ∈ N , and finally ϕ (cid:96) ∈ Ker (cid:16) H N ( V ( v ∞ )) − E (0) dis ( v ∞ ) (cid:17) . Nevertheless, by an adaptation of Theorem 2.2 to the discretized case, weknow that since v maximizes G ( k ) r, α , then it has a k th mixed excited state Γ with density respecting (cid:82) ρ Γ α i = r i . (cid:3) Remark 6.5. When Ω = R d , the situation is different. Assume also E ( k ) dis ( v ) =0 . We apply (30) with Q = { v i < (cid:15) } for some (cid:15) > , which is not empty since inf v = 0 . This yields G ( k ) r, α ( v ) − c R (cid:54) − min (cid:32) , (cid:80) v i (cid:62) (cid:15) r i (cid:80) v i <(cid:15) r i (cid:33) || v − (cid:15) || (cid:96) r, α and we conclude by letting (cid:15) → , G ( k ) r, α ( v ) − c R (cid:54) − min (cid:32) , (cid:80) v i > r i (cid:80) v i (cid:54) r i (cid:33) || v || (cid:96) r (cid:54) − min (cid:18) , (cid:80) v i > r i N (cid:19) || v || (cid:96) r . The problem is that we are not able to find a strictly positive lower bound for (cid:80) v i > r i , which would provide coercivity. Remark 6.6. A natural norm on potentials is the gauge invariant quotientnorm || v || (cid:96) r / ∼ = inf µ ∈ R || v − µ || (cid:96) r , where v ∼ u when v − u is constant. If in (30) we take Q = (cid:8) i ∈ I (cid:12)(cid:12) v i = inf v (cid:9) , µ = inf v we obtain G ( k ) r, α ( v ) − c R (cid:54) − (cid:88) v i > inf v ( v i − inf v ) r i = − || v − inf v || (cid:96) r (cid:54) − || v || (cid:96) r / ∼ . Hence G ( k ) r, α is coercive in the (cid:96) r / ∼ norm, but this is not a convenient normbecause by definition we do not control the constant. Building Kohn-Sham potentials: proofs of Corollary 3.6 andTheorem 3.9. In this section, we show how Theorem 3.5 yields approximate v -representability when Ω is unbounded.6.4.1. The mixed states case.Proof of Corollary 3.6. Consider copies of the cube C n := [ − /n, /n ) d , cen-tered on the grid points of ( Z /n ) d . Take Ω n to be the union of all thosecubes ( C ni ) i ∈ I n which are included in Ω ∩ B n , where B r is the ball of radius n . They form an increasing sequence Ω n ⊂ Ω n +1 . We choose α ni := C ni . Weapply Theorem 3.5 to r ( n ) := c n r ρ Ω n with c n := N/ (cid:82) Ω n ρ . The condition |{ ρ = 0 }| = 0 ensures that r ( n ) i > for any i ∈ I n , n ∈ N . We denote by v n ∈ (cid:96) ∞ ( I, R ) the maximizer of G ( k ) r ( n ) , α n , we have (cid:82) ρ Γ n α ni = c n (cid:82) ρα ni andwe apply Theorem 3.2 for each n .Each Γ n lives in S N mix (Ω n ) and each H N ( (cid:80) v ni α ni ) is an operator of L a (Ω Nn ) ,but by “digging pits” close to where the density is localized, we can createa potential V n = − λ n Ω n + (cid:80) i v ni α ni with λ n large enough so that we keepthe same properties for systems living in an unbounded domain Ω .For k = 0 , we take Γ n a minimizer of F α n , (0) ( r ( n ) ) and apply Theorem 3.2. (cid:3) UILDING INVERSE POTENTIALS 35 The pure states case.Proof of Theorem 3.9. We assume N (cid:62) and only restrict to N = 1 at theend of the argument. We define the map κ : H a (Ω N ) ∩ (cid:110)(cid:82) | Ψ | = 1 (cid:111) −→ ( R + ) | I | Ψ (cid:55)−→ (cid:0)(cid:82) ρ Ψ α i (cid:1) i ∈ I . It is C ∞ and since I is finite, Ran d Ψ κ is closed for any Ψ ∈ L a (Ω N ) becausethe image lives in a finite-dimensional space. Now d Ψ G has closed rangefor any Ψ ∈ ∧ N H (Ω) because its target space is finite-dimensional. Wecompute, for any ϕ ∈ H a (Ω N ) and any v ∈ R | I | , (cid:88) i ∈ I v i ((d Ψ κ ) ϕ ) i = 2 N Re (cid:88) i ∈ I v i (cid:90) Ω N (Ψ ϕ )( x , . . . , x N ) α i ( x )d x · · · d x N = 2 Re (cid:90) Ω N N (cid:88) j =1 (cid:32)(cid:88) i ∈ I v i α i (cid:33) ( x j )(Ψ ϕ )( x , . . . , x N )d x · · · d x N = 2 Re (cid:10) ϕ, (cid:0) (cid:80) i V ( v ) i (cid:1) Ψ (cid:11) Furthermore, h : Ψ (cid:55)→ (cid:104) Ψ , H N (0)Ψ (cid:105) is C ∞ with differential (d Ψ h ) ϕ = 2 Re (cid:104) ϕ, H N (0)Ψ (cid:105) . Let Ψ ∈ H a (Ω N ) with unit norm be a minimizer of F α , (0) ( r ) , then (cid:82) ρ Ψ α i = r i for any i ∈ I . We apply [44, Prop 43.19 p291], ensuring the exis-tence of Lagrange multipliers ( v i ) i ∈ I ∈ (cid:0) (cid:96) ( I, R ) (cid:1) ∗ = (cid:96) ∞ ( I, R ) such that H N ( V ( v ))Ψ = 0 weakly, or such that (cid:0) (cid:80) i V ( v ) i (cid:1) Ψ = 0 weakly and v (cid:54) = 0 .We prove by contradiction that the second case is impossible, so let usassume that (cid:0) (cid:80) i V ( v ) i (cid:1) Ψ = 0 . This is where we need N = 1 , which implies V ( v )Ψ = 0 and v i (cid:82) ρ Ψ α i = v i r i , but since r i > , we conclude that v i =0 , which is a contradiction. If Conjecture 3.10 holds, then (cid:0) (cid:80) i V ( v ) i (cid:1) Ψ = 0 implies (cid:0) (cid:80) i V ( v ) i (cid:1) = 0 a.e and v = 0 , and this avoids the second case forall N .We hence know that H N ( V ( v ))Ψ = 0 (32)with V ( v ) ∈ L p (Ω) , p as in (1), so Ψ is in the domain of H N ( V ( v )) andconsequently is an eigenvalue. Since Ω is bounded, Ψ is in the discretespectrum. (cid:3) Proofs of Section 4: Lemma 4.1 and Proposition 4.2. Proof of Lemma 4.1. i ) We define D N := dim Ker R (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) . We choose a basis (Ψ I ) I ∈I tot of Ker R (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) composed of antisymmetric productsof one-body real orbitals. We compute (cid:90) R d ( N − (Ψ I Ψ J ) ( x, x , . . . , x N )d x · · · d x N = δ I,J N − (cid:32)(cid:88) i ∈ I ϕ i ( x ) (cid:33) + δ I ∪· J = { i,j } N − ϕ i ( x ) ϕ j ( x ) . Figure 7. Vanishing crossing terms between two configurations.The second term does not vanish if there is exactly one of the degenerateone-body energy levels of I and J differing by exactly one particle, the otherparticles of the outer levels should have the same distribution. Indeed, asillustrated in Figure 7, if the difference belongs to two different levels (left),the energies are different so those terms actually do not appear, and if theenergies are equal but differences belong to two different levels, then thereare more than one difference. We deduce that (cid:90) R d ( N − (Ψ I Ψ J ) ( x, x , . . . , x N )d x · · · d x N (33) = δ I,J N − (cid:32)(cid:88) i ∈ I ϕ i ( x ) (cid:33) + N − δ ∃ (cid:96),i,jI ∪· J = I ( (cid:96) ) ∪· J ( (cid:96) ) = { i,j } I ( t ) = J ( t ) ∀ t (cid:54) = (cid:96) ϕ i ( x ) ϕ j ( x ) . Finally we can rewrite + δ v E ( k ) ( u ) = µ k − m k (cid:0) P (cid:0) (cid:80) i u i (cid:1) P (cid:1) , where P is the projector onto Ker R (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) . Using (33), thematrix elements of P (cid:0) (cid:80) i u i (cid:1) P are (cid:10) Ψ I , (cid:0) (cid:80) i u i (cid:1) Ψ J (cid:11) = (cid:90) ( δ IJ ρ in + M ϕ,IJ ) u. Defining the elementwise action of the integral on matrices (cid:0)(cid:82) u M ϕ (cid:1) IJ := (cid:82) u M ϕ,IJ , we have + δ v G ρ ( u ) = (cid:90) ( ρ in − ρ ) u + µ k − m k (cid:18)(cid:90) u M ϕ (cid:19) .ii ) As we see in (33) and Figure 7, (cid:82) R d ( N − Ψ I Ψ J = 0 when I and J havemore than one particle difference. iii ) In this case, M ϕ is diagonal by ii ) . (cid:3) Proof of Proposition 4.2. i ) The optimality condition sup u + δ v G ( k ) ρ ( u ) (cid:54) is equivalent to µ k − m k (cid:18)(cid:90) Ω u M ϕ (cid:19) (cid:54) (cid:90) Ω ( ρ − ρ in ) u (cid:54) µ M k − k (cid:18)(cid:90) Ω u M ϕ (cid:19) , (34)for all u ∈ L p + L ∞ . These are the Euler-Lagrange inequalities, replacingequalities because of the degeneracies. Applying this condition (34) to thesequences u n := n d B /n ( x ) and − u n , and taking n → + ∞ , we obtain (18).However, (18) does not imply (34). More generally, for a function-valuedreal symmetric matrix S , and for k, K ∈ N , µ k ( S ( x )) (cid:54) (cid:54) µ K ( S ( x )) a.e. is local and it does not imply µ k (cid:0)(cid:82) uS (cid:1) (cid:54) (cid:54) µ K (cid:0)(cid:82) uS (cid:1) for all u , UILDING INVERSE POTENTIALS 37 which is global. A counterexample is the D × D diagonal matrix S ii = − if x ∈ [( i − /D [ and S ii = 1 otherwise, indeed with u = 1 we obtain min σ (cid:0)(cid:82) uS (cid:1) = ( D − /D > . ii ) We define N in := (cid:82) ρ in the number of “inner” particles, and N out := N − N in the number of outer particles. If N out = 1 , we are in the situa-tion where the only degeneracy comes from a one-body degeneracy at theone-body eigenspace Span( φ i ) (cid:54) i (cid:54) D N . At almost every x ∈ R d , M ϕ = (cid:0) φ . . . φ D N (cid:1) T (cid:0) φ . . . φ D N (cid:1) , all the eigenvalue are except one, which is (cid:80) D N i =1 φ i . Indeed, the columns are proportional hence the rank is one, andthe last eigenvalue is equal to the trace. Applying (18) and integrating yields N out (cid:54) (cid:90) µ M k − k ( M ϕ ) = (cid:26) (cid:82) (cid:80) D N i =1 φ i = D N if k = m k otherwise.We deduce that m k = k . iii ) We apply Lemma 4.1 iii ) . (cid:3) Appendix 1: two remarks on density functionals In this appendix, we make two remarks on excited states functionals.First, we notice that the inner problem in the Levy-Lieb functional (4) isfinite and has an optimizer. Lemma 6.7. Take k ∈ N with k (cid:62) , A ⊂ L a ( R dN ) , dim C A = k and ρ ∈ L ( R d , R + ) such that (cid:82) ρ = N and √ ρ ∈ H ( R d ) . There exists Ψ ∈ H a ( R dN , C ) such that Ψ ⊥ A and ρ Ψ = ρ , and the infimum inf Ψ ∈ A ⊥ ρ Ψ = ρ (cid:104) Ψ , H N (0)Ψ (cid:105) is finite and attained.Proof. We take a frame Span ( φ , . . . , φ k ) = A , and then we consider theHarriman-Lieb [21, 34] orbitals ϕ , . . . , ϕ N ∈ H ( R d , C ) , which are such that ρ ∧ Nj =1 ϕ j = ρ . We then use the orthonormalization procedure [28, Corol-lary 1.3] with the functions φ , . . . φ k , ϕ , . . . , ϕ N so that there exists func-tions g , . . . , g k , f , . . . , f N ∈ C ∞ c ( R , R ) such that φ e ig ( x ) , . . . φ k e ig k ( x ) , ϕ e if ( x ) , . . . , ϕ N e if N ( x ) is an orthonormal familly. Eventually, by defining Ψ := ∧ Nj =1 ϕ j e if j ( x ) ∈ H ( R dN , C ) , we have ρ Ψ = ρ and Ψ ∈ A . We concludethat the set (cid:8) Ψ ∈ H ∩ A ⊥ (cid:12)(cid:12) ρ Ψ = ρ (cid:9) is not empty.The minimum is attained by an adaptation of the proof of the case k = 0 which is [34, Theorem 3.3]. (cid:3) Then, we present a remark about the other possibility of defining theLevy-Lieb functional for excited states, which is ˜ F ( k ) ( ρ ) := inf A ⊂ H a ( R dN )dim C A = k +1 ∃ Ψ ∈ A,ρ Ψ = ρ max Ψ ∈ Aρ Ψ = ρ (cid:104) Ψ , H N (0)Ψ (cid:105) . The advantage of this one is that we can directly prove that ˜ F ( k ) ( ρ ) is finite.However, we now show why it seems to be not of much use. The first Levy-Lieb functional provides an upper bound to the energy E ( k ) ( v ) (cid:54) inf ρ (cid:62) , (cid:82) ρ = N √ ρ ∈ H (cid:18) F ( k ) ( ρ ) + (cid:90) vρ (cid:19) . However, there are potentials v ∈ L p + L ∞ such that E ( k ) ( v ) < sup ρ (cid:62) , (cid:82) ρ = N √ ρ ∈ H (cid:18) ˜ F ( k ) ( ρ ) + (cid:90) vρ (cid:19) and hence ˜ F ( k ) ( ρ ) does not provide a lower bound to the energy. Indeed,otherwise we would have ˜ F ( k ) ( ρ ) (cid:54) E ( k ) ( v ) − (cid:82) vρ for all v ∈ L p + L ∞ and all ρ (cid:62) such that (cid:82) ρ = N , √ ρ ∈ H , hence ˜ F ( k ) ( ρ ) (cid:54) inf v (cid:0) E ( k ) ( v ) − (cid:82) vρ (cid:1) .We assumed that (cid:82) B r ρ > without loss of generality. With v n = n B r forinstance, then E ( k ) ( v n ) − (cid:82) v n ρ ∼ − n (cid:82) B r ρ → −∞ when n → + ∞ , and thiswould imply that ˜ F ( k ) ( ρ ) = −∞ for all ρ . Appendix 2: maximizing + δ v G ( k ) ρ ( u ) in the two-folddegenerate case Here we show that in the case dim Ker (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) = 2 = p = q ,we can reduce the -dimensional optimization problem (6) of maximizing + δ v G ( k ) ρ to a 1-dimensional problem. The interaction w is general. Thiswould enable to further accelerate the ODA algorithm.Take Ψ and Φ real such that they form a real orthonormal basis of thedegenerate level Ker R (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) , we have M k = m k + 1 and let usdefine (cid:15) := 1 if k = m k and (cid:15) := − if k = M k . Then + δ v G ( m k ) ρ ( u ) = min a,b ∈ C | a | + | b | =1 (cid:90) u ( ρ a Ψ+ b Φ − ρ ) (35)and + δ v G ( M k ) ρ ( u ) has the same formula but with a maximization. With theparametrization a = (cos α ) e iη , b = (sin α ) e i ( η + β ) , we have ρ (cos α )Ψ+ e iβ (sin α )Φ = (cos α ) ρ Ψ + (sin α ) ρ Φ + 2 sin α cos α cos βρ Ψ , Φ = ρ Ψ + ρ Φ ρ Ψ − ρ Φ α + sin 2 α cos βρ Ψ , Φ and (cid:90) uρ (cos α )Ψ+ e iβ (sin α )Φ = (cid:90) u ρ Ψ + ρ Φ A cos 2 α + B sin 2 α cos β where A u := 12 (cid:90) u ( ρ Ψ − ρ Φ ) , B u := (cid:10) Ψ , (cid:0) (cid:80) i u i (cid:1) Φ (cid:11) = (cid:90) uρ Ψ , Φ . This yields + δ v G ( m k ) ρ ( u ) = (cid:90) u (cid:18) ρ Ψ + ρ Φ − ρ (cid:19) + min α,β ∈ [0 , π ] ( A u cos 2 α + B u cos β sin 2 α ) . UILDING INVERSE POTENTIALS 39 Optimizing over α yields the optimal value α u ∈ π N +arctan ( B u (cos β ) /A u ) and using the classical formulas for cos arctan and sin arctan , we get A u cos 2 α u + B u cos β sin 2 α u = ± (cid:112) A u + (cos β ) B u . Finally optimizing over β gives β u = 0 , so + δ v G ( k ) ρ ( u ) = (cid:90) (cid:18) ρ Ψ + ρ Φ − ρ (cid:19) u − (cid:15) (cid:112) A u + B u , (36)and for B (cid:54) = 0 , cos 2 α u = − (cid:15)A u (cid:112) A u + B u , sin 2 α u = − (cid:15)B u (cid:112) A u + B u In order to compute the supremum over directions u , we consider theLagrangian L ( u, λ ) := + δ v G ( k ) ρ ( u ) − λ (cid:0)(cid:82) u − (cid:1) , which is C ∞ when u is nota constant, and the Euler-Lagrange equation is λu ∗ = ρ Ψ + ρ Φ − ρ − (cid:15) ( ρ Ψ − ρ Φ ) A u ∗ + 2 B u ∗ ρ Ψ , Φ (cid:112) A u ∗ + B u ∗ = ρ Ψ + ρ Φ − ρ + (cos 2 α u ∗ ) ρ Ψ − ρ Φ α u ∗ ) ρ Ψ , Φ , hence the optimal direction belongs to the directions u θ := c θ (cid:18) ρ Ψ + ρ Φ − ρ + (cos 2 θ ) ρ Ψ − ρ Φ θ ) ρ Ψ , Φ (cid:19) , where θ ∈ [0 , π ] . The prefactor c θ (cid:62) is chosed such that (cid:82) u = 1 . Wecould then reintroduce this expression back into (36) and reoptimize over θ to find an equation that θ has to verify, to maximize + δ v G ( k ) ρ ( u θ ) , but thisrelation is quite involved. Denoting by θ ∗ the optimizing angle, which is thesame for k = m k and k = M k , we reduced the problem to a circle searchsup || u || L =1 + δ v G ( k ) ρ ( u ) = max θ ∈ [0 , π ] + δ v G ( k ) ρ ( u θ ) = + δ v G ( k ) ρ ( u θ ∗ )= (cid:90) (cid:18) ρ Ψ + ρ Φ − ρ (cid:19) u θ ∗ − (cid:15) (cid:113) A u θ ∗ + B u θ ∗ . Appendix 3: computations for the implementation Here we provide complementary computations for the description of thealgorithm of Section 5.6.6. Optimal direction over mixed states. We now provide the detailsof the computations linked to the problem (22). We express the cost func-tion in terms of the parameters Γ IJ and provide its gradient, needed in theimplementation of ODA.Once again we take the notations of Section 4.1, where we approximate Ker (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) by (cid:110) Ψ ∈ H a (Ω N ) (cid:12)(cid:12) (cid:12)(cid:12)(cid:12) E v (Ψ) − E ( k ) ( v ) (cid:12)(cid:12)(cid:12) (cid:54) t (cid:111) , having dimension |I tot | = |I out | , and we define Q out := Span R (cid:0) ∧ i ∈ I Ψ I (cid:1) I ∈I out .Mixed states are decomposed into Γ = (cid:80) I,J ∈I tot Γ IJ | Ψ I (cid:105) (cid:104) Ψ J | , where (Γ IJ ) IJ is a real positive matrix with unit trace, and ρ Γ = N (cid:88) I,J ∈I tot Γ IJ (cid:90) Ω N − Ψ I Ψ J = ρ in + (cid:88) I,J ∈I out Γ IJ ρ out IJ , where ρ out IJ := N out (cid:82) R d ( N out − Ψ I Ψ J . The problem (22) can be reformulatedby P ( v ) = min Γ ∈S ( Q out )Γ (cid:62) , Tr Γ=1 e ( e in + E v (Γ) − E ( k )( v ) ) T (cid:90) Ω ( ρ Γ + ρ in − ρ ) , where e in := (cid:80) i ∈ I in E i . For any function g : M → R , where M is a smoothmanifold modeled on a Hilbert space H with scalar product (cid:104)(cid:104)· , ·(cid:105)(cid:105) , we recallthat we can define the gradient ∇ x g ∈ T x M (cid:39) H by Riesz’ theorem, via (d x g ) a = (cid:104)(cid:104)∇ x g, a (cid:105)(cid:105) . The energy of a mixed state is E v (Γ) = (cid:88) I,J ∈I tot Γ IJ (cid:42) Ψ J , N (cid:88) i =1 ( − ∆ i + v ( x i )) Ψ I (cid:43) = e in + (cid:88) I,J ∈I out Γ IJ e vIJ , where for I, J ∈ I out , e vIJ := (cid:42) Ψ J , N out (cid:88) i =1 ( − ∆ i + v ( x i )) Ψ I (cid:43) = δ IJ (cid:88) i ∈ I E i + δ I ∪· J = { i,j } (cid:90) ϕ j ( − ∆ + v ) ϕ i Hence ( ∇ Γ E v ) IJ = e vIJ which does not depend on Γ . The function we opti-mize is f (Γ) = e (cid:18) e in − E ( k )( v )+ (cid:80) K,L ∈I out Γ KLevKL (cid:19) T (cid:18) (cid:90) ( ρ in − ρ ) + 2 (cid:88) I,J ∈I out Γ IJ (cid:90) ( ρ in − ρ ) ρ KL + (cid:88) I,J ∈I out K,L ∈I out Γ IJ Γ KL (cid:90) ρ IJ ρ KL (cid:19) , and has gradient ( ∇ Γ f ) IJ = e (cid:18) e in − E ( k )( v )+ (cid:80) K,L ∈I out Γ KLevKL (cid:19) T (cid:32) (cid:90) ( ρ in − ρ ) ρ out IJ + 2 (cid:88) K,L ∈I out Γ KL (cid:90) ρ out IJ ρ out KL + T − e vIJ (cid:16) e in − E ( k ) ( v ) + (cid:88) K,L ∈I out Γ KL e vKL (cid:17) × (cid:18) (cid:90) ( ρ in − ρ ) + 2 (cid:88) K,L ∈I out Γ KL (cid:90) ( ρ in − ρ ) ρ out KL + (cid:88) K,L ∈I out X,Y ∈I out Γ KL Γ XY (cid:90) ρ out KL ρ out XY (cid:19)(cid:33) . UILDING INVERSE POTENTIALS 41 Hence, to launch the ODA, one has to compute e in − E ( k ) ( v ) , (cid:90) ( ρ in − ρ ) , (cid:90) ( ρ in − ρ ) ρ out KL , e vKL , (cid:90) ρ out KL ρ out XY . Optimization over pure states. We take the notations of Section 4.1,and here give an extra computation concerning (23). The familly ( ∧ i ∈ I in ϕ i ∧ Ψ I ) I ∈I out is a basis of Ker R (cid:0) H N ( v ) − E ( k ) ( v ) (cid:1) . We can represent the eigen-functions of the k th N -body level by complex columns vectors C , Ψ = (cid:88) I ∈I tot C I Ψ I = (cid:94) j ∈ I in ϕ j (cid:88) I ∈I out C I (cid:94) i ∈ I ϕ i . The coefficients verify (cid:80) I ∈I out | C I | = 1 and the set of such C ’s forms aGrassmann manifold. Then the density is ρ Ψ = ρ in + (cid:88) I ∈I out | C I | (cid:88) i ∈ I ϕ i + 2 (cid:88) I,J ∈I out ∃ (cid:96),i,jI ∪· J = { i,j } (cid:0) Re C I C J (cid:1) ϕ i ϕ j , and is invariant under U (cid:0) |I out | (cid:1) transformations of C . References [1] G. Accorto, P. Brandolini, F. Marino, A. Porro, A. Scalesi, G. Colò,X. Roca-Maza, and E. Vigezzi , First step in the nuclear inverse Kohn-Shamproblem: From densities to potentials , Phys. Rev. C, 101 (2020), p. 024315.[2] A. Alfonsi, R. Coyaud, V. Ehrlacher, and D. Lombardi , Approximation of op-timal transport problems with marginal moments constraints , Math. Comput, (2020).[3] I. Babuška and J. M. Melenk , The partition of unity method , Int. J. Numer.Meth. Eng, 40 (1997), pp. 727–758.[4] T. J. Callow, N. N. Lathiotakis, and N. I. Gidopoulos , Density-inversionmethod for the Kohn-Sham potential: Role of the screening density , J. Chem. Phys,152 (2020), p. 164114.[5] É. Cancès , SCF algorithms for HF electronic calculations , in Mathematical modelsand methods for ab initio quantum chemistry, vol. 74 of Lecture Notes in Chem,Springer, Berlin, 2000, ch. 2, pp. 17–43.[6] É. Cancès and C. Le Bris , Can we outperform the DIIS approach for electronicstructure calculations? , Int. J. Quantum Chem, 79 (2000), pp. 82–90.[7] , On the convergence of SCF algorithms for the Hartree-Fock equations , M2ANMath. Model. Numer. Anal., 34 (2000), pp. 749–774.[8] E. Cancès and N. Mourad , A mathematical perspective on density functional per-turbation theory , Nonlinearity, 27 (2014), p. 1999.[9] J. Chayes, L. Chayes, and E. H. Lieb , The inverse problem in classical statisticalmechanics , Comm. Math. Phys, 93 (1984), pp. 57–121.[10] J. Chayes, L. Chayes, and M. B. Ruskai , Density functional approach to quantumlattice systems , J. Stat. Phys, 38 (1985), pp. 497–518.[11] S.-K. Chua and R. L. Wheeden , Estimates of best constants for weighted poincaréinequalities on convex domains , Proc. London Math. Soc., 93 (2006), pp. 197–226.[12] R. Coyaud , in preparation , PhD thesis, 2021.[13] N. Ferré, M. Filatov, M. Huix-Rotllant, and C. Adamo , Density-functionalmethods for excited states , Springer, 2016.[14] K. F. Freed and M. Levy , Direct first principles algorithm for the universal elec-tron density functional , J. Chem. Phys, 77 (1982), pp. 396–398.[15] L. Garrigue , Unique continuation for many-body Schrödinger operators and theHohenberg-Kohn theorem , Math. Phys. Anal. Geom, 21 (2018), p. 27. [16] L. Garrigue , Some properties of the potential-to-ground state map in quantum me-chanics , arXiv e-prints, (2020), p. arXiv:2012.04054.[17] L. Garrigue , Unique continuation for many-body Schrödinger operators and theHohenberg-Kohn theorem. II. The Pauli Hamiltonian , Doc. Math, (2020).[18] R. Gaudoin and K. Burke , Lack of Hohenberg-Kohn theorem for excited states ,Phys. Rev. Lett, 93 (2004), p. 173001.[19] G. Giuliani and G. Vignale , Quantum Theory of the Electron Liquid , CambridgeUniversity Press, 2005.[20] A. Gonis and M. Däne , On the v-representability of ensemble densities of electronsystems , J. Phys. Chem. Solids, 116 (2018), pp. 100–112.[21] J. E. Harriman , Orthonormal orbitals for the representation of an arbitrary density ,Phys. Rev. A, 24 (1981), pp. 680–682.[22] M. F. Herbst, A. Levitt, G. Kemlin, S. Sirajdine, E. Berquist, L. Ponet,and Tzsuzsi , Juliamolsim/dftk.jl: v0.2.3 , Dec. 2020.[23] P. Hohenberg and W. Kohn , Inhomogeneous electron gas , Phys. Rev, 136 (1964),pp. B864–B871.[24] D. S. Jensen and A. Wasserman , Numerical methods for the inverse problem ofdensity functional theory , Int. J. Quantum Chem, 118 (2018), p. e25425.[25] B. Kanungo, P. M. Zimmerman, and V. Gavini , Exact exchange-correlationpotentials from ground-state electron densities , Nat. Commun, 10 (2019), pp. 1–9.[26] W. Kohn and L. J. Sham , Self-consistent equations including exchange and corre-lation effects , Phys. Rev. (2), 140 (1965), pp. A1133–A1138.[27] A. Kumar, R. Singh, and M. K. Harbola , Universal nature of different methodsof obtaining the exact Kohn-Sham exchange-correlation potential for a given density ,J. Phys. B, 52 (2019), p. 075007.[28] O. Lazarev and E. H. Lieb , A smooth, complex generalization of the Hobby-Ricetheorem , Indiana Univ. Math. J., 62 (2013), pp. 1133–1141.[29] M. Levy , Universal variational functionals of electron densities, first-order densitymatrices, and natural spin-orbitals and solution of the v -representability problem ,Proc. Natl. Acad. Sci. USA, 76 (1979), pp. 6062–6065.[30] M. Lewin , Geometric methods for nonlinear many-body quantum systems , J. Funct.Anal., 260 (2011), pp. 3535–3595.[31] M. Lewin, E. H. Lieb, and R. Seiringer , Statistical mechanics of the uniform elec-tron gas , J. Éc. polytech. Math, 5 (2018), pp. 79–116.[32] , The local density approximation in density functional theory , Pure Appl. Anal,2 (2019), pp. 35–73.[33] , Universal functionals in density functional theory , arXiv preprintarXiv:1912.10424, (2019).[34] E. H. Lieb , Density functionals for Coulomb systems , Int. J. Quantum Chem, 24(1983), pp. 243–277.[35] , Density functional methods in physics , NATO ASI Series B, 123 (1985).[36] E. H. Lieb and M. Loss , Analysis , vol. 14 of Graduate Studies in Mathematics,American Mathematical Society, Providence, RI, 2nd ed., 2001.[37] J. M. Melenk and I. Babuška , The partition of unity finite element method: basictheory and applications , in Research Report/Seminar für Angewandte Mathematik,vol. 1996, Eidgenössische Technische Hochschule, Seminar für Angewandte Mathe-matik, 1996.[38] J. R. Moreno, G. Carleo, and A. Georges , Deep learning the Hohenberg-Kohnmaps of density functional theory , Phys. Rev. Lett, 125 (2020), p. 076402.[39] T. Naito, D. Ohashi, and H. Liang , Improvement of functionals in density func-tional theory by the inverse Kohn-Sham method and density functional perturbationtheory , J. Phys. B, 52 (2019), p. 245003.[40] M. Penz, A. Laestadius, E. I. Tellgren, and M. Ruggenthaler , Guaranteedconvergence of a regularized Kohn-Sham iteration in finite dimensions , Phys. Rev.Lett, 123 (2019), p. 037401.[41] D. Schnieders and J. Neugebauer , Accurate embedding through potential recon-struction: A comparison of different strategies , J. Chem. Phys, 149 (2018), p. 054103. UILDING INVERSE POTENTIALS 43 [42] L. O. Wagner, T. E. Baker, E. Stoudenmire, K. Burke, and S. R. White , Kohn-Sham calculations with the exact functional , Phys. Rev. B, 90 (2014), p. 045109.[43] Q. Wu and W. Yang , A direct optimization method for calculating density func-tionals and exchange–correlation potentials from electron densities , J. Chem. Phys,118 (2003), pp. 2498–2509.[44] E. Zeidler , Nonlinear functional analysis and its applications. III : Variationalmethods and optimization , Springer Science & Business Media, 2013.[45] , Nonlinear functional analysis and its applications. IV : Applications to math-ematical physics , Springer Science & Business Media, 2013. CERMICS, École des ponts ParisTech, 6 and 8 av. Pascal, 77455 Marne-la-Vallée, France Email address ::