Density-potential inversion from Moreau-Yosida regularization

For a quantum-mechanical many-electron system, given a density, the Zhao-Morrison-Parr method allows to compute the effective potential that yields precisely that density. In this work, we demonstrate how this and similar inversion procedures mathematically relate to the Moreau-Yosida regularization of density functionals on Banach spaces. It is shown that these inversion procedures can in fact be understood as a limit process as the regularization parameter approaches zero. This sheds new insight on the role of Moreau-Yosida regularization in density-functional theory and allows to systematically improve density-potential inversion. Our results apply to the Kohn-Sham setting with fractional occupation that determines an effective one-body potential that in turn reproduces an interacting density.


I. INTRODUCTION
In cases where the one-body ground-state density ρ gs of a quantum-mechanical many-electron system is known, the Zhao-Morrison-Parr (ZMP) method [1; 2], among others [3][4][5], allows to determine the effective potential of a reference system that reproduces exactly that density. Such inverseproblem methods can either be used to directly match experimental data [6], or to gain valuable insight into approximations in density-functional theory (DFT) [7][8][9][10]. The resulting density-potential map has been recently studied from a mathematical point of view in great detail [11; 12]. In particular, these methods allow to determine the effective potential, v KS , of a Kohn-Sham (KS) system [13] of non-interacting particles and thus also the exchange-correlation potential, v xc . Specialized methods for inversion in KS systems have also been devised [14; 15] and the time-dependent case recently received attention as well [16]. The ZMP method itself has been implemented and used on numerous occasions, mainly in the context of KS-DFT [17][18][19][20][21][22][23].
In this work, we focus on methods for non-relativistic Nelectron quantum systems. Here, v will be a one-body potential and V the corresponding lifted N -body operator, i.e., V = N j=1 v(r j ). Now take ρ gs as a ground-state density of some Hamiltonian H and H ref 0 a reference Hamiltonian. A motivating problem for our discussion is as follows: How does the (effective) potential v needs to be chosen such that H ref 0 + V has exactly the given ground-state density ρ gs ? For H modeling an interacting system in an external potential v ext and taking H ref i.e., the KS potential that contains the external potential of H, the Hartree potential v H = ρ gs * | · | −1 , and the exchange-correlation potential. This effectively maps the interacting problem to a non-interacting one. Contrary to that, in the case H = H ref 0 + V , our result addresses the direct density-potential inversion. a) Electronic mail: andre.laestadius@oslomet.no The ZMP method was originally derived using a penalty parameter λ with the density constraint appearing in terms of a Coulomb integral. For a fixed value λ > 0 this method defines v λ (r) = λ ρ λ (r ) − ρ gs (r ) where ρ gs is a ground-state density of the Hamiltonian H. The pair (ρ λ , v λ ) is determined self-consistently and the procedure is repeated for increasing values of λ. Finally, careful numerical extrapolation λ → ∞ yields a potential v λ → v that has the ground-state density ρ gs , as required [2]. In the following, we will demonstrate how the ZMP and related methods can be justified in a mathematical precise manner by relying on the Moreau-Yosida regularization of a (not necessarily universal) density functional on a suitably chosen function space. The possibility of such a connection between the ZMP method and Moreau-Yosida regularization was already mentioned in Ref. [24], but has not yet been made explicit.

A. Outline
In Section II A, we briefly review the ZMP method for the reader's convenience. Next, in Section II B, we discuss the necessary convex-analysis concepts and Moreau-Yosida regularization with applications to DFT in mind. We then formulate our rather general setting in Section II C, to which we will apply the density-potential inversion procedure. Our main results are discussed in Section III, which contains the derivation of the ZMP method with an additional correction term (Theorem 3). This is achieved by first formulating an abstract density-potential inversion procedure (Procedure 1 and also Procedure 2) by exploiting a property of the derivative of the Moreau-Yosida regularized functional (Theorem 2). We then specialize the abstract procedure to obtain the ZMP method both on bounded (Theorem 3), and on unbounded domains (Theorem 4). In Section III D, we describe some of the numerical experiments. Section IV summarizes the established results and gives an outlook to possible future work. The proofs of our results may be found in Section V.

A. The Zhao-Morrison-Parr method
The main idea of the ZMP method [1; 2] is to rephrase the pointwise constraint ρ(r) = ρ gs (r) to the equivalent D(ρ − ρ gs ) = 0, where D(ρ) denotes the Hartree energy of ρ (see Eq. (12)), and add it to the one-body energy with a penalty parameter λ > 0. We then arrive at the minimization problem where we have set ρ λ (r) = N j=1 |φ λ j (r)| 2 and v ext is the external potential. Requiring that the derivative of the above functional with respect to the orbitals φ λ j vanishes, we obtain the Kohn-Sham-type equations of the ZMP method [2]. These equations are then solved for λ → ∞. We may summarize the method as follows.

Remark 1.
(i) We have included v ext (r)ρ λ (r) dr and D(ρ λ ) in the minimization problem such that v λ targets v xc only, as can be seen from Eq. (2). The Hartree potential in Eq.
(2) is sometimes substituted for a Fermi-Amaldi term, as for example in Ref. [2].
(ii) The choice of D(·) as a penalty term in the minimization problem is rather ad hoc at this point, and Ref. [1] suggested to use the L 2 -norm instead. One can imagine that a wide variety of functionals should work and we will later show how the corresponding different ZMP schemes arise from different choices for the basic function spaces.
(iii) Note that λ is here just a parameter introduced to penalize for densities differing too much from ρ gs . In our interpretation below, the regularization parameter ε replaces the penalty parameter λ.

B. Moreau-Yosida regularization in DFT
Differentiability of the exact universal (interacting or noninteracting) density functionals cannot be guaranteed, they are in fact everywhere discontinuous in the standard formulation [25]. A way to achieve differentiability of convex functionals, introduced in the context of DFT in Ref. [24] (see Refs. [26; 27] for recent reviews), is by Moreau-Yosida regularization. This program was extended by two of the present authors in Refs. [28; 29] to different DFT settings, and was also used to show convergence of a modified KS iteration on finite lattices [30; 31]. In this section, we briefly review the basic properties of Moreau-Yosida regularization and their consequences for DFT.
We start in the very general setting of a Banach space X, containing all densities under consideration, that is reflexive and where both X and its dual X * are strictly convex. A normed space is strictly convex if the function ρ → ρ 2 is strictly convex. All Lebesgue spaces L p , 1 < p < ∞, and in general every Hilbert space verify these assumptions, they are even uniformly convex. We note that the non-reflexive spaces L 1 and L ∞ are excluded.
In the following, an important role will be played by the duality mapping J that maps elements from X to X * in a canonical way. It is defined as The duality mapping is always homogeneous and in the given setting it is also single-valued and bijective [32, Prop. 1.117], so J −1 : X * → X is well-defined. The relevance of convex analysis to DFT was recognized very early [33], and so a series of definitions from this very useful field are in order. We refer the reader to various textbooks [32; 34; 35] for more information. As mentioned above, exact density functionals are (typically) not differentiable, however, a more general concept of differentiation is applicable. The subdifferential ∂f : X → P(X * ) (mapping to the power set of X * ) of a convex functional f : X → R ∪ {+∞} at ρ ∈ X is defined as the set Intuitively, the convex set ∂f (ρ) collects all the supporting hyperplanes of the graph of f at the point (ρ, f (ρ)). If f happens to be continuously differentiable at a point ρ, then ∂f (ρ) is the singleton set {f (ρ)}. If ∂f (ρ) = ∅, then we say that f is subdifferentiable at ρ. Analogously, the notion of the superdifferential ∂f can be introduced for concave functions f via ∂f = −∂(−f ). Similarly, whenever ∂f (ρ) = ∅, then we say that f is superdifferentiable at ρ.
The Legendre transform f * : The functionals f * and g * are always convex (even if f and g are not). The famous biconjugation theorem says that whenever f is proper, lower semicontinuous and convex, then f * is proper, (weak- * ) lower semicontinuous and convex, and furthermore (f * ) * = f . This result is very important for the convex formulation of DFT.
The Moreau-Yosida regularization f ε : X → R with a fixed ε > 0 of a convex, lower semicontinuous functional f : X → R∪{+∞} is given by the lower envelope of the parabola ρ → 1 2ε ρ 2 X tracing along the graph (ρ, f (ρ)). In a formula this means Since the parabola is strictly convex and X is reflexive, the above infimum is attained at a unique point. Consequently, the following definition of the proximal mapping Π ε f : X → X makes sense, Both the regularization of a functional and the proximal mapping are exemplified in Fig. 1. We now collect all the relevant properties of the Moreau-Yosida regularization, the proofs of which may be found in Refs. [24; 28; 32; 35].

lower semicontinuous functional. Then the following properties hold true for the Moreau-Yosida
(v) The proximal mapping Π ε f is always single-valued and obeys the limit Π ε f (ρ) → ρ as ε → 0.
Let us briefly discuss the consequences of these properties on a concrete example furnished by the convex and (weakly) lower semicontinuous [33], but the everywhere discontinuous [25] The Lieb functional may also be defined through the concave ground-state energy E : X * → R via Legendre transformation as which is sometimes referred to as the Lieb variational principle [27]. As usual, the dual pairing between potential and density means v, ρ = v(r)ρ(r) dr. Now, consider the Moreau-Yosida regularization F ε of F . The regularized ground-state energy E ε (more specifically, the energy corresponding to the regularized functional) given by −E ε (−·) = (F ε ) * can then be expressed using property ((viii)) above simply as The one-parameter family of convex and (Gâteaux-) differentiable functionals {F ε } ε>0 converges pointwise and increasingly to the original, unregularized F from below as ε → 0 by property ((iii)) above. The regularized ground-state energy E ε (v) also converges to the true ground-state energy E(v) according to Eq. (6), and since J is just the derivative of the concave functional − 1 2 v 2 X * . This also means that the degeneracy of the ground state, manifested through multiple ground-state densities ∂E(v), is retained exactly for the regularized case. These facts reflect the "lossless" character of the Moreau-Yosidaregularization, i.e., no information is lost in the regularization process. Note, however, that the Moreau-Yosida regularization of F ε does not make E ε differentiable. Nevertheless, it makes E ε superdifferentiable everywhere, while in the unregularized case E is superdifferentiable only at potentials v ∈ X * that support a ground state (a.k.a. binding potentials). In such a case, ∂E(v) is the set of densities coming from ensemble ground states associated to v.
By the biconjugation theorem, the well-known duality relation between F and E [33], sometimes called the Hohenberg-Kohn variational principle, also holds in the regularized setting, According to Theorem 1 ((i)), F ε is finite everywhere (i.e., even at densities which are not representable by a wavefunction), so the minimizer ρ in Eq. (8) might be unphysical: ρ can be negative, might not integrate to N , or |∇ √ ρ| 2 dr might diverge (see [33]). This is in contrast to minimizers of Eq. (7), which are of course always physical. The way out from this apparent difficulty is the following. We know from general convex-analysis considerations [32,Prop. 2.33] that ρ is an optimizer in Eq. (8) if and only if −v ∈ ∂F ε (ρ), in other words (Theorem 1 ((iv))), if and only if Combining this with Theorem 1 ((vi)), this is equivalent to (recall that we use the notation Inverting this equation, we find that a minimizer ρ of Eq. (8) is always of the form We call ρ ε the proximal density of ρ. It is important to note that ρ ε will always be physical: ρ ε ≥ 0, ρ ε dr = N and |∇ √ ρ ε | 2 dr < ∞. This is because ρ ε is indeed a minimizer of Eq. (7) by Theorem 1 ((vii)). Again, we see the "lossless" character of Moreau-Yosida regularization: from a possibly unphysical minimizer ρ (a.k.a. a quasidensity) of Eq. (8), we can always reconstruct a unique physical density We can conclude that solving the original minimization problem Eq. (7) is completely equivalent to solving a regularized problem Eq. (8). One simply needs to combine Eq. (6) and Eq. (9) to obtain the ground-state energy E(v) and the respective ground-state density. Similar considerations hold true for the (mixed state) kineticenergy functional which is also convex and lower semicontinuous [33], but discontinuous everywhere. One obtains the regularized functional T ε with equivalent properties as before. However, it is important to note at this point that the functional used in conventional KS-DFT is not convex [33], so that the Moreau-Yosida regularization (or the convex-analysis approach) is not directly applicable in this case. Therefore, Moreau-Yosida regularized KS methods will always involve fractional occupation numbers since one considers density matrices in Eq. (10). We note in passing, that in the case ρ is (non-interacting) v-representable with a non-degenerate ground state, then T (ρ) = T S (ρ).
We are now ready to recognize that the ZMP minimization problem initially posed in Section II A really has the form of a Moreau-Yosida regularization, where the regularization In order to show that this is not just a accidental similarity and to be able to really formulate the ZMP method rigorously in terms of a Moreau-Yosida regularization, we have to find the suitable setting with a convex functional.

C. Model setting
In KS-DFT the universal density functional (introduced above in Eq. (4)) where D(ρ) = D(ρ, ρ), D(ρ, σ) = 1 2 are the Hartree energy and Coulomb inner product. The crucial object in Eq. (11) is the exchange-correlation energy E xc that accounts for the non-classical contributions and that actually gets defined through Eq. (11). E xc has been referred to as "nature's glue" [36].
To set the stage, we will relate ground-state densities and their corresponding potentials through convex energy functionals. First, we introduce the density functional F(ρ), given below as a constrained-search functional over density matrices, F CS (ρ), and an additional term, a convex (and lower semicontinous) functional G(ρ) > −∞, i.e.,

Remark 2.
(i) The functional G will typically be non-linear and differentiable.

(iii) Relevant for the forthcoming discussion is the choice
, so that the potential from the inversion procedure is just the exchange-correlation potential. This will be further explained in Section III A below.
(iv) Another possibility is to follow Ref. [4] and set H ref 0 = 0 and for G(ρ) use, e.g., one of the functionals where α > 0. These functionals are all convex (and lower semicontinuous).
The density functional F can then be used to compute the (reference) ground-state energy E F (v) of a system described by F (H ref 0 and G) with an additional potential v, Apart from a difference in sign, E F (v) is the convex conjugate of F(ρ), and as such, a concave function. We can also revert back from E F (v) to F(ρ) by duality as just like in Eq. (5). The condition for a density being a minimizer in Eq. (14) can be equivalently rephrased as (see Section II B) − v ∈ ∂F(ρ).
This equation already constitutes a density-potential map and will serve as the starting point for the development of a more practical method. The inversion of this relation, mapping from potentials to densities, entails solving the corresponding KStype (SCF-) equation with potential v and determining the ground-state density (or densities, in case of degeneracy). By the reciprocity theorem [32,Prop. 2.33], this amounts to the superdifferential of the convex conjugate functional E F , which is just the condition for v being a maximizer in Eq. (15). In order to establish a relation to the notion of vrepresentability, we consider the case when G is differentiable, such that For a potential −v ∈ ∂F the difference to ∂F CS just appears as a fixed shift +G and the cardinality of the set ∂F CS is the same for ∂F. This allows to discern three qualitatively different cases when it comes to solutions of Eq. is v-representable and if we are in a setting where the Hohenberg-Kohn theorem holds [37].
(iii) Finally, the subdifferential can also contain multiple elements in cases where the density is non-uniquely vrepresentable, a situation recently discovered in finite lattice system with ground-state degeneracy [38]. Now, returning to the general case, the situation would be much simpler if the functional F(ρ) could be assumed to be functionally (i.e. Gâteaux-) differentiable, since then the subdifferential always contains exactly one element (the functional derivative). As we discussed above, this is typically not the case, and in particular not the case for exact functionals. To remedy the situation, we recall that the desired differentiability do hold for the regularized functionals. It will be shown below that by appropriately combining these ideas, Moreau-Yosida regularization serves as a basis for the rigorous reformulation of ZMP and related methods.

A. Abstract density-potential inversion
The main concrete example for the forthcoming discussion is the (fractional occupation number) Kohn-Sham scheme, which we choose to model with the functional that was already mentioned in Remark 2 ((iii)). The functional F KS is convex and lower semicontinuous for every v ext ∈ X * . In the KS setting, ρ ∈ X is noninteracting v-representable if and only if ∂F KS (ρ) = ∅. Note that we have chosen the functional F KS to be such that only the exchange-correlation energy is left out as the remaining, unknown contribution. Then, comparing with Eq. (17) above, the set of KS potentials reads . This is because ( v ext , · ) ≡ v ext and D (ρ) is just the Hartree potential. Using this language, the density-potential inversion consists of finding a representative of −∂F KS (ρ gs ), where ρ gs is a ground-state of the interacting problem. The following fundamental result says that it is always possible to find this potential as the limit of the functional derivative −(F ε KS ) (ρ) of the regularized functional F ε KS (ρ) as ε → 0. We state the result in its full generality, since we want to apply it in different situations, i.e., not just for F KS .

Theorem 2.
Suppose that X is a strictly convex Banach space and X * a uniformly convex one. Let F : X → R ∪ {+∞} be a convex, lower semicontinuous functional. Let ρ gs ∈ X be such that the subdifferential is nonempty, ∂F(ρ gs ) = ∅. Setting ρ ε := Π ε F (ρ gs ) as before, we have that the (strong) limit of the sequence in X * is the unique element v ∈ −∂F(ρ gs ) ⊂ X * with minimal norm.
The proof can be found in an even more general form in Section V. Recall that (F ε ) (ρ gs ) always exists due to to regularization. To implement the procedure suggested by the theorem, one needs to determine the proximal density ρ ε and v ε := 1 ε J(ρ ε − ρ gs ) ∈ −∂F(ρ ε ), and let ε → 0. Of course, in general, neither step is trivial.
Here, the relation v ε ∈ −∂F(ρ ε ) is equivalent to ρ ε ∈ ∂E F (v ε ) by the reciprocity relation, where E F (v) = inf ρ [F(ρ) + v, ρ ]. This step avoids the necessity to calculate the proximal mapping. For instance, if F = F KS , then solving ρ ε ∈ ∂E FKS (v ε ) for ρ ε amounts to solving the (fractional occupation number) Kohn-Sham equations with v ε in place of the exchange-correlation potential. We return to this case in more detail below.
All this suggests the following abstract self-consistent scheme. Suppose that besides the ground-state density ρ gs ∈ X, the initial v ε 0 and ρ ε 0 are given. In a general step we compute ρ ε i ∈ ∂E F (v ε i ) and then update the potential via When sufficiently converged, this is repeated for smaller values of ε to facilitate the limit ε → 0 in the form of an extrapolation. We summarize the above considerations in a form of an algorithm. We use mixing to facilitate convergence, like in the optimal-damping algorithm for electronic-structure calculations [39].

Procedure 1 (abstract density-potential inversion algorithm).
Fix a strictly decreasing sequence ε k → 0 and a mixing parameter 0 < µ < 1. For a given v-representable ρ gs ∈ X, choose corresponding initial values v k 0 . In each i-iteration determine the ground-state density ρ k i ∈ ∂E F (v k i ) and calculate the succeeding potential as Switch to the next k-iteration when v k i is sufficiently converged to some v k and finally let v k → v as k → ∞.
Here, the initial value v k 0 for the effective potential may be chosen to be the zero potential or the result from the previous k-iteration.
It must be noted that in general the convergence of the given procedure is not assured (see Section III D for a numerical test and some remarks about the choices on µ and the ε-sequence). Especially the occurrence of degeneracies-when ∂E F (v) includes multiple elements to choose from-will be problematic. Yet the potentials that produce such degeneracies are expected to be rare in potential space and thus they are unlikely to be hit by an iteration step v k i . Procedure 1 can be simplified by letting the mixing parameter µ go to zero, while keeping α := µε −1 k constant, thereby dropping the extrapolation steps in k. This gives an alternative procedure. Procedure 2 (simplified abstract density-potential inversion algorithm). Fix a parameter 0 < α < 1. For a given vrepresentable ρ gs ∈ X choose an initial value for the potential v 0 . In each iteration step determine the ground-state density ρ i ∈ ∂E F (v i ) and calculate the succeeding potential as

Stop when v i is sufficiently converged.
Clearly, the choice of function space X determines the duality mapping J that is crucial for concrete realizations and we now explore a few possibilities that will include two variants of the ZMP method as special cases.

B. Interpretation of the Zhao-Morrison-Parr method
We already presented the original ZMP approach to densitypotential inversion in Section II A. It is now fairly straightforward to interpret the method in a rigorous way as the densitypotential inversion following Eq. (18). The two procedures developed for abstract density-potential inversion would then already serve as possible refinements. Additionally, the discussed cases also add correction terms to the usual ZMP method.
Furthermore, by setting λ = (4πε) −1 , the first term in the right-hand side of Eq. (20) corresponds exactly to the ZMP method Eq. (1), while the second term is a correction that depends on the shape of Ω.
The setting Ω = R 3 is the one usually presented in molecular DFT. In order to formulate the ZMP method on a unbounded domain, we may proceed as follows. Let X = H −1 (R 3 ) and which gives J(ρ) = v. The integral kernel 1 4π |r| −1 e −γ|r| appearing here is known as the Yukawa potential. We can therefore conclude Theorem 4 (ZMP method on unbounded domains). Let X = H −1 (R 3 ) and X * = H 1 (R 3 ). Then the potential update step of Eq.
With this we can proceed like in Theorem 3.

C. Density-potential inversion on L p spaces
Before presenting our numerical results in Section III D, we wish to discuss the Hilbert-space setting X = L 2 (Ω), X * = L 2 (Ω), Ω ⊆ R d . This choice implies J = J −1 = id. Irrespective of Ω being bounded or not, X does not include the space L 1 ∩ L 3 , so the setting is unsuited for dealing with densities in the continuum.
This changes in a discrete setting that represents a finite lattice, since for X = R M , M the number of vertices, all L p norms, p ≥ 1, are equivalent. Using the identity map in Procedure 1, we get the naïve update scheme which tells us to choose a positive (repulsive) potential where the current density is larger than the target density and a negative (attractive) one if it is smaller. Similarly, Procedure 2 is just This update scheme already appeared in the reverseengineering procedure in Ref. [42] on quantum rings that is itself inspired by the density-potential inversion method of Ref. [3]. The same method was used in Ref. [43], but with an adaptive parameter α i . The discrete L 2 Hilbert-space setting will be further discussed with a numerical experiment in Section III D where we compare the different methods.
We also briefly comment on the choice X = L 3 (Ω), X * = L 3/2 (Ω), Ω ⊆ R d . This example does not consist of Hilbert spaces and it is chosen especially in order to include L 1 ∩L 3 in X [33]. It was previously featured in Ref. [28], where the regularized Kohn-Sham iteration was generalized to such density spaces. An example of how the duality mapping looks like for the case L p ([0, 1]) can be found in Ref. [44,Prop. 3.14], but it is not a simple expression even in just one dimension (d = 1).

D. Numerical illustration
In this section, we consider a numerical example based on a single system that serves as a basic proof-of-concept of the density-potential inversion method. The system under consideration for our numerical exploration is a quantum ring with M = 50 lattice sites, next-neighbor hopping τ = −1, and a single spinless particle. The space choice is X = X * = R M with the standard L 2 -norm. We fix a periodic target density ρ gs and solve for the corresponding effective potential by three different methods. The first is just applying a standard non-linear optimization algorithm (BFGS, in its scipy implementation) to Eq. (5), while the other two are Procedures 1 and 2. The parameter choice is µ = 0.05, α = 0.5, and the relatively short ε-sequence (1, 0.7, 0.4, 0.1). (Note that these parameters cannot be considered universal, a different setting requires a recalibration, so a more careful choice including system size and other parameters seems desirable.) Then all three potentials are again put into the Schrödinger equation and we solve for the ground-state density that is then compared to the target density. Since the potential from the BFGS method turns out to be the most accurate one, the other two potentials are compared to this one. We summarize the results in Fig. 2. For the chosen tolerance, convergence was achieved in 136 iterations in Procedure 2 while in every ε-iteration between 43 and 70 steps were required to converge the inner i-iteration (with a total of 211 iterations) in Procedure 1. This means that within this very specific example, Procedure 1, our generalization of the ZMP method, has no benefit compared to the simpler Procedure 2, neither regarding iteration speed nor accuracy. Yet, this finding was only established for the given example and other setups might show a more pronounced difference between the two procedures built upon Moreau-Yosida regularization. Still, both procedures have comparable overall convergence speed, while the reference optimization algorithm BFGS has a considerably longer runtime (by a factor of 10 in the chosen example).
Note that here the L 2 -norm was mainly chosen for practicability, since it yields the trivial duality mapping J = id. A different, non-uniform choice of norm could for example focus on specific areas of interest and help to suppress problems with low-density regions.

IV. SUMMARY AND FUTURE WORK
Starting from the original ZMP method [1; 2], we have formulated a general framework aimed at understanding the density-potential mapping using the Moreau-Yosida regularization scheme applied to DFT. For the convenience of the reader, we reviewed the Moreau-Yosida regularization in detail and summarized its consequences for DFT (Section II B). Next, we formulated our rather general setting that allows to handle a wide variety of density functionals (Section II C). Our main results were stated in Section III, which also contain our reformulation of the ZMP method. We first proposed an abstract density-potential inversion procedure (Procedure 1) by exploiting a property of the derivative of the Moreau-Yosida-regularized density functional (Theorem 2). Further, Procedure 2 yields an even simpler density-potential inversion scheme that was already considered in the literature before [42; 43]. We then specialized the abstract procedure to obtain the ZMP method both on bounded domains (Theorem 3) and on unbounded domains (Theorem 4). In the bounded case, we found that a correction term to the original ZMP method, Eq. (1), is necessary, which accounts for the finite size of the domain. In the unbounded case, the integration kernel of the original ZMP method gets changed to a Yukawa potential. This means that a refinement of the current ZMP inversion algorithms seems entirely possible, e.g., by considering such additional terms or by switching to different, more appropriate spaces for density and potential. Section III D finally describes a small proof-of-concept of the different inversion schemes within a numerical experiment.
We stress that the focus of this work is not the immediate development of a new and more efficient practical densitypotential inversion method, but to understand the ZMP method and related techniques within a generalized setting based on the Moreau-Yosida regularization. In conclusion, our findings provide a framework for devising and analyzing density- Comparison of different density-potential-inversion methods implemented on a quantum ring with one particle and M = 50 lattice sites. A target density was given, then different inversion schemes (Procedure 1, Procedure 2 and the standard non-linear optimization algorithm BFGS) were applied to get the potential and to calculate again the ground-state density from there. The left panels show the respective functions, the right panels show the difference to the target density (for densities) and to the BFGS result (for potentials) in a logarithmic scale. potential inversion procedures. It would be also interesting to see how the proposed modifications of the ZMP method perform in practice. Possible modifications are the inclusion of the additional terms derived in Theorems 3 and 4 or the choice of an altogether different density space X, which influences the method through the corresponding duality mapping.

V. PROOFS
We first restate Theorem 2 in a slightly more general way. The theorem itself is then an immediate consequence if one just remembers that uniform convexity of a space implies it is reflexive and that if a space is reflexive, also its dual space has this property.

Lemma 1.
Suppose that X, X * are both strictly convex and reflexive. Let f : X → R ∪ {+∞} be a convex, lower semicontinuous functional and x ∈ X such that the subdifferential ∂f (x) is nonempty. Then there is a unique ξ ∈ ∂f (x) ⊆ X * with minimal norm that is the weak limit of for ε 0. Moreover, if X * is uniformly convex then strong convergence holds.
Proof. Since f is convex and lower semicontinuous, its subdifferential is a maximal monotone operator from X to X * [32,Th. 2.43]. Then the first equality and weak convergence follow directly from a result for maximal monotone operators [32, Prop. 1.146 (iv)]. We will just add a note about the notation in the reference for the reader's orientation. In Ref. [32], the parameter ε is replaced by λ, the duality mapping J is called f , and the subdifferential ∂f is the operator A, while the gradients (f ε ) are A λ . The connection with regularization is made in Ref. [32,Sec. 2.2.3]. That the density is in the domain of the subdifferential means that the subdifferential is not empty at this point. That there exists a unique element with minimal norm in ∂f follows from a basic theorem of functional analysis [45,Cor. 5.1.19].
We finally show the set membership in Eq. (21). Taking the definition of the proximal mapping as the minimizer in Eq. (3) and translating that into the condition that the subdifferential includes zero at the minimum, we have Here we first used that the subdifferential is additive for convex functions and then that ∂ 1 2 · 2 X = J [32, Ex. 2.32].
Proof of Theorem 3. As usual X * = H 1 0 (Ω) is the Sobolev space with norm v X * = ( v 2 L 2 + ∇v 2 L 2 ) 1/2 (here ∇ is the weak gradient operator) and zero trace on ∂Ω. The Poincaré inequality [40, 6.30] gives an equivalent norm v X * = ∇v L 2 that we will henceforth choose. The dual space X = H −1 (Ω) consists of distributions as explained in Ref. [40, 3.12]. Now since with the inverse duality mapping and a potential v ∈ X * it must hold v, J −1 (v) = v 2 X * = ∇v 2 L 2 = v, −∆v , we just get J −1 = −∆ D , the Dirichlet-Laplace operator on Ω. (This also corresponds to the Riesz-Fréchet map since we are in the Hilbert space setting.) To obtain v = J(ρ) we must thus solve −∆v = ρ, v| ∂Ω = 0, i.e., Poisson's equation with homogeneous Dirichlet boundary conditions. If ρ is continuous and we assume a twice continuouslydifferentiable solution v to exist, then this solution is given by an integral involving the Green function for the domain Ω This Green function can be expressed as G(r, r ) = Φ(r−r )− φ r Ω (r ) with Φ(r) = 1/(4π|r|) the fundamental solution of Laplace's equation and φ r Ω (r ) the corrector function solving ∆ φ r Ω (r ) = 0 on Ω with boundary condition φ r Ω (r ) = Φ(r − r) on ∂Ω. Inserting this into Eq. (18) gives the stated formula. The convergence to the exchange-correlation potential is then an application of Theorem 2.
Proof of Theorem 4. The only difference between the previous proof is a change in the duality map (because the function space has been changed). The usual norm of X * = H 1 (R 3 ) is v X * = ( v 2 L 2 + ∇v 2 L 2 ) 1/2 but for any γ > 0 one has v X * = (γ 2 v 2 L 2 + ∇v 2 L 2 ) 1/2 as an equivalent norm. Then the condition for the duality map is v, J −1 (v) = γ 2 v 2 L 2 + ∇v 2 L 2 = v, γ 2 v + v, −∆v and thus J −1 = −∆ + γ 2 . Consequently, application of J means solving the screened Poisson equation.

DATA AVAILABILITY STATEMENT
The numerical example of Section III D was implemented in a small Python script that is part of the public domain and can be found at https://github.com/ERC-REGAL/REGAL/.