January 01, 1970
The inverse Kohn–Sham (KS) problem seeks a local effective potential whose noninteracting ground state reproduces a prescribed electron density. Although many inversion formulations and schemes have been developed, they are often formulated in disparate languages, including reduced variational optimization, penalty regularization, response-based iteration, and PDE-constrained optimization. In this work, we develop a unified framework for inverse KS theory in two steps. First, we identify the fixed-density noninteracting constrained search embedded in exact density functional theory as the natural variational anchor of inverse KS inversion. In this setting, the KS potential appears as the variational dual object associated with density reproduction, reducing to the familiar multiplier picture in regular regimes. Second, building on this anchor, we classify major inversion formulations according to how the KS state equations and density-reproduction condition are treated as objectives, constraints, penalties, or feasibility relations within the optimization architecture. Within this framework, the Wu–Yang formulation appears as a reduced exact-multiplier realization, the Zhao–Morrison–Parr construction as a quadratic-penalty relaxation, and PDE-constrained approaches as explicit state-constraint formulations at the orbital level. Rather than comparing inversion formulations primarily at the algorithmic level, the present work uses variational and convex-analytic structure as its organizing language. This viewpoint clarifies the roles of additive-constant ambiguity, asymptotic normalization, nonsmooth variational structure, and weak-gap instability across inversion formulations. The aim is not a complete existence theory of inverse KS representability, but a unified mathematical framework for understanding how major inversion formulations are connected and why their numerical failure modes are closely related.
Density functional theory (DFT) is one of the central frameworks for the electronic structure of interacting many-body systems in chemistry and condensed-matter physics. Its formal basis is provided by the Hohenberg–Kohn theorems [1], which establish, under appropriate assumptions, a one-to-one correspondence between the ground-state density and the external scalar potential, up to an additive constant. Within the Kohn–Sham (KS) construction [2], this correspondence motivates the introduction of an auxiliary noninteracting system whose ground-state density matches that of the interacting problem, while exchange and correlation are absorbed into an effective one-body potential.
A persistent question in exact and practical DFT is the inverse Kohn–Sham problem: given a prescribed target density \(\rho_{\mathrm{tar}}(\mathbf{r})\), determine an effective local potential \(v_\mathrm{s}(\mathbf{r})\) whose noninteracting ground state reproduces that density [3], [4]. In practice, target densities may be obtained from highly accurate electronic-structure calculations, and inversion then serves both as a conceptual probe of the density–potential correspondence and as a practical route to reference effective and exchange-correlation potentials [3], [5]. It also provides a useful setting for diagnosing the limitations of approximate functionals and for generating benchmark data for data-driven electronic-structure modeling [5]–[7].
Over the past three decades, several influential inversion procedures have been developed, including the Zhao–Morrison–Parr (ZMP) method [8] and the Wu–Yang variational formulation [9]. A broader class of inversion and stabilization strategies has since been developed, including variational reconstructions, basis-set-stabilized formulations, and response-based schemes [10]–[14]. Closely related wavefunction- and reduced-density-matrix-based reconstructions of KS potentials have also been introduced [15]–[17]. These methods are often introduced in rather different mathematical languages, including reduced optimization over potentials, quadratic-penalty formulations, response-based iteration, orbital-level reconstruction, and explicit state-equation approaches [3], [4]. A main goal of the present work is to place these formulations within a unified variational framework. In this sense, the aim is not primarily an algorithmic comparison, but a structural reorganization in which variational and convex-analytic ideas serve as the main language for relating the principal inversion formulations.
The perspective adopted here has two connected levels. At the first level, we argue that the inverse KS problem should be understood not merely as the inversion of a nonlinear density-to-potential map, but as the fixed-density noninteracting constrained-search problem already embedded in exact DFT. More precisely, inverse KS theory is anchored here in the density-constrained inner noninteracting variational problem that appears within the Levy–Lieb and Lieb formulations of exact DFT [18], [19]. In this setting, the KS potential appears naturally as the multiplier associated with density reproduction, so that inverse KS theory is placed directly within the variational architecture of exact DFT rather than outside it. This interpretation is particularly transparent when exact DFT is organized into parallel interacting and noninteracting ensemble variational hierarchies: inverse KS theory then appears as the problem of recovering the supporting potential of \(T_s[\rho]\) at a prescribed density.[20]
At the second level, building on this variational anchor, we develop a broader structural classification of inverse KS formulations. The key observation is that inverse KS inversion is governed by two defining relations: the KS state equations and the density-reproduction condition. Major inversion formulations may then be classified according to how these two relations are assigned the roles of objective, constraint, penalty, or feasibility relation within the optimization architecture. In this sense, Wu–Yang appears as a reduced exact-multiplier formulation [9], ZMP as a quadratic-penalty relaxation [8], [21], and PDE-constrained approaches as explicit state-constraint formulations at the orbital level [4], [22], [23]. The same structural perspective also naturally suggests augmented-Lagrangian enrichments, which interpolate between exact multiplier enforcement and pure penalty regularization, as well as enlarged all-at-once formulations in which state consistency and density reproduction are treated within a broader coupled optimization architecture. It also accommodates limiting feasibility-type formulations in which both defining relations are treated primarily as constraints.
This viewpoint does not remove the familiar analytical difficulties of inversion. Issues including noninteracting \(v\)-representability, nonuniqueness up to constants, nonsmooth variational structure, weak-gap instability, and discretization-induced ill-conditioning remain central in the inversion literature [3], [24]–[27]. The point of the present framework is instead to place these phenomena in a common mathematical setting and to identify how they enter different inversion formulations.
The aim is not a complete theorem-driven treatment of inverse KS representability in full generality. Rather, we seek to provide a unified mathematical framework within which the major inversion formulations, their regularity assumptions, and their characteristic numerical failure modes may be understood in common terms.
A common heuristic description of inverse Kohn–Sham theory treats the problem as the inversion of a nonlinear forward map \[\mathcal{M}: v \mapsto \rho[v].\] This notation is useful at an intuitive level, but it should be understood only schematically, in a regular regime where the ground-state density associated with \(v\) is well defined. Its domain and range are not meant to coincide with the full ambient potential and density spaces, and the potential is understood only up to the usual additive-constant ambiguity. More importantly, this picture suppresses the variational structure of exact density functional theory. In particular, it obscures the fact that the effective Kohn–Sham potential arises naturally from the fixed-density noninteracting constrained search, rather than merely as the output of an abstract inverse map [18], [19].
In this section, we reformulate the inverse Kohn–Sham problem in a function-space variational setting. The purpose is not to replace the familiar physical interpretation, but to place it within the constrained-search architecture of exact DFT. In this formulation, the inverse problem appears as the fixed-density noninteracting constrained-search problem, whose value later becomes the noninteracting kinetic functional \(T_\mathrm{s}[\rho]\), and its connection to multiplier theory, regularity, and density–potential duality becomes explicit.
Throughout this work, we use \(\langle \cdot,\cdot\rangle\) for single-particle inner products and natural dual pairings, while Dirac bra–ket notation \(\langle \cdot|\hat{A}|\cdot\rangle\) is reserved for many-body operator expectation values. The precise meaning is determined by context. We write \(\rho\) for a generic density, \(\rho_{\mathrm{tar}}\) for a prescribed target density, and \(v\) for a generic multiplier or dual potential variable. Once this dual variable is identified with the effective Kohn–Sham potential, we denote it by \(v_\mathrm{s}\).
Let \(\Omega \subseteq \mathbb{R}^3\) denote the spatial domain. Depending on the physical setting, one may take \(\Omega=\mathbb{R}^3\) for isolated systems or a periodic torus for extended systems. We write \(H^1(\Omega)\) for the usual Sobolev space \[H^1(\Omega)=\left\{\phi\in L^2(\Omega)\mid \nabla \phi \in [L^2(\Omega)]^3 \right\},\] which provides the natural single-particle energy space for KS orbitals.
For the density variable, we adopt the standard Lieb setting [19]. The natural Banach space is \[X=L^1(\Omega)\cap L^3(\Omega),\] equipped with the norm \[\|\rho\|_X=\|\rho\|_{L^1}+\|\rho\|_{L^3}.\] Its dual may be identified with \[X^*=L^\infty(\Omega)+L^{3/2}(\Omega),\] which serves as a natural ambient space for external and effective scalar potentials [19]. In the present work, this dual space is used not only for external potentials in the sense of Lieb’s formulation, but also as the natural ambient space for the dual variables associated with the noninteracting kinetic functional. In particular, the KS potential \(v_\mathrm{s}\) is most naturally viewed as an element of \(X^*\) at the level of variational duality, even though additional assumptions are needed to identify such a dual element with a regular local Schrödinger potential.
A physically admissible electron density must satisfy more than integrability. We therefore consider densities in the set \[\mathcal{I}_N := \left\{ \rho\in X \;\middle|\; \rho\ge 0,\; \int_\Omega \rho = N,\; \nabla \sqrt{\rho}\in [L^2(\Omega)]^3 \right\}.\] Here and below, inequalities such as \(\rho\ge 0\) are understood almost everywhere. The condition \(\nabla \sqrt{\rho}\in L^2\) is the familiar von Weizsäcker-type regularity requirement. It is central in mathematical DFT because it controls finiteness of kinetic energy and is closely tied to the existence of finite-energy states reproducing the given density [19].
At the many-body level, one works over an admissible antisymmetric \(N\)-electron wavefunction class with finite kinetic energy and associated density in \(\mathcal{I}_N\). For present purposes, the structural point is that the density variable lives naturally in a Banach space, whereas the state variables live in Sobolev-type energy spaces. This mixed Banach–Sobolev structure is one source of the analytical subtlety of the inverse KS problem [19], [26].
At the most basic level, the inverse KS problem asks whether a prescribed target density can be represented as the noninteracting ground-state density of some effective local potential. In a regular regime, this may be viewed schematically as a map \[\rho \mapsto v_\mathrm{s}[\rho],\] defined on an appropriate representable class and only modulo the usual additive-constant ambiguity.
Operationally, the correspondence is defined through the single-particle eigenvalue problem \[\left( -\frac{1}{2}\nabla^2+v_\mathrm{s}[\rho](\mathbf{r}) \right)\phi_i(\mathbf{r}) = \varepsilon_i \phi_i(\mathbf{r}), \qquad i=1,\dots,N, \label{eq:ks95state}\tag{1}\] together with the density reconstruction condition \[\rho(\mathbf{r})=\sum_{i=1}^N |\phi_i(\mathbf{r})|^2. \label{eq:density95reproduction}\tag{2}\] Thus the inverse problem may be viewed operationally as the search for a local potential whose ground-state orbitals reproduce a prescribed density.
The mathematical status of this correspondence is subtler than the schematic notation suggests. In standard ground-state DFT, the cleanest Hohenberg–Kohn uniqueness statements are obtained in nondegenerate settings, while degeneracy and related nonuniqueness phenomena require more careful treatment of the density–potential relation and its domain of validity [1], [24]–[26]. Inverse KS theory inherits these subtleties, and practical inversions are further complicated by ill-posedness and instability in finite discretizations [3], [4].
Accordingly, the discussion below is understood primarily in the regular zero-temperature, fixed-\(N\), local-potential regime in which a sufficiently well-behaved inverse KS correspondence may be meaningfully discussed. Our aim is not to settle the full existence theory of the map \(\rho\mapsto v_\mathrm{s}[\rho]\), but to explain why, in the regime where such a map is meaningful, the effective KS potential arises naturally from the multiplier structure of the noninteracting constrained search.
Let \(\rho_{\mathrm{tar}}\in \mathcal{I}_N\) be a prescribed target density. The corresponding fixed-density noninteracting constrained search is \[T_\mathrm{s}[\rho_{\mathrm{tar}}] = \inf_{\Phi \mapsto \rho_{\mathrm{tar}}} \langle \Phi | \hat{T} | \Phi \rangle, \label{eq:Ts95rhotar}\tag{3}\] where the infimum is taken over admissible noninteracting \(N\)-electron states \(\Phi\) with density \(\rho_\Phi=\rho_{\mathrm{tar}}\). Thus inverse KS theory is not, at core, an unconstrained search for a potential. It is the study of the constrained minimization problem that defines the noninteracting kinetic energy at a prescribed density [18], [19].
To analyze 3 , one introduces the density constraint explicitly and studies the associated Lagrangian. Writing the constraint abstractly as \(\rho_\Phi=\rho_{\mathrm{tar}}\), together with the usual orthonormality constraints on the occupied orbitals underlying \(\Phi\), one is led formally to \[\begin{align} \mathcal{L}(\Phi,v_\mathrm{s},\Lambda) :={}& \langle \Phi|\hat{T}|\Phi\rangle + \langle v_\mathrm{s},\rho_\Phi-\rho_{\mathrm{tar}}\rangle \\ &\quad - \operatorname{Tr}\!\big[\Lambda(\langle \phi_i,\phi_j\rangle-\delta_{ij})\big]. \end{align}\] where \(v_\mathrm{s}\) is the multiplier associated with the density constraint and \(\Lambda\) collects the orbital orthonormality multipliers.
Stationarity with respect to the orbitals yields the familiar Euler–Lagrange structure \[\left( -\frac{1}{2}\Delta + v_\mathrm{s}(\mathbf{r}) \right)\phi_i(\mathbf{r}) = \sum_j \Lambda_{ij}\phi_j(\mathbf{r}),\] which, after diagonalization of the Hermitian matrix \(\Lambda\), becomes \[\left( -\frac{1}{2}\Delta + v_\mathrm{s}(\mathbf{r}) \right)\phi_i(\mathbf{r}) = \varepsilon_i \phi_i(\mathbf{r}).\] Thus the inverse problem may be read in a precise variational sense: one fixes the density and recovers the multiplier field associated with constrained minimization of the noninteracting kinetic energy.
This interpretation is conceptually important. It shows that the effective potential of inverse KS theory is not an ad hoc reconstruction variable introduced from outside the variational structure of DFT. It is the multiplier naturally attached to the fixed-density constrained search 3 . In this way, the inverse problem already appears as an intrinsic variational object before one passes to any particular inversion algorithm.
At the same time, this multiplier picture also clarifies why inverse KS theory is analytically delicate. The density constraint is posed in a Banach-space setting appropriate to admissible densities, whereas the state variable \(\Phi\) lives in a Sobolev-type energy space. The resulting Banach–Sobolev coupling is precisely what makes the constrained search mathematically nontrivial and what underlies representability, regularity, and stability issues in practice [19], [26].
The multiplier origin of \(v_\mathrm{s}\) acquires its full meaning only when placed in the broader setting of the Levy–Lieb constrained-search formulation [18], [19]. In that framework, the exact ground-state energy is written as \[E_0 = \inf_{\rho\in \mathcal{I}_N} \left\{ F_{\mathrm{LL}}[\rho] + \langle v_{\mathrm{ext}},\rho\rangle \right\},\] where \[F_{\mathrm{LL}}[\rho] = \inf_{\Psi \mapsto \rho} \langle \Psi | \hat{T}+\hat{W} | \Psi \rangle.\] Thus the interacting many-body problem already has the form of an outer variational problem over densities together with an inner constrained search over states compatible with those densities.
The fixed-density minimization problem introduced above is not a separate construction specific to inversion. Rather, it is precisely the fixed-density instance of the noninteracting constrained search defining the KS kinetic functional, \[T_\mathrm{s}[\rho] = \inf_{\Phi\mapsto \rho} \langle \Phi | \hat{T} | \Phi \rangle.\] For inverse KS theory, one studies this constrained search at the prescribed density \(\rho=\rho_{\mathrm{tar}}\) and analyzes the associated multiplier structure.
The KS construction preserves the Levy–Lieb variational logic while reorganizing the universal functional into noninteracting, Hartree, and exchange–correlation parts. Accordingly, one may write the forward KS problem as \[E_0 = \inf_{\rho\in \mathcal{I}_N} \left\{ \inf_{\Phi\mapsto \rho} \langle \Phi|\hat{T}|\Phi\rangle + E_\mathrm{H}[\rho] + E_{\mathrm{xc}}[\rho] + \langle v_{\mathrm{ext}},\rho\rangle \right\}.\] Thus the forward problem already contains an outer–inner variational structure: the outer level selects the energetically optimal density, while the inner level identifies the noninteracting state of minimal kinetic energy compatible with that density.
In a sufficiently regular regime, combining outer stationarity with the inner multiplier structure yields the standard KS relation \[v_\mathrm{s}[\rho](\mathbf{r}) = v_{\mathrm{ext}}(\mathbf{r}) + v_\mathrm{H}[\rho](\mathbf{r}) + v_{\mathrm{xc}}[\rho](\mathbf{r}),\] modulo the usual additive constant ambiguity. In this way, the effective potential appearing as a multiplier in the inner constrained search is identified with the usual KS potential of the forward theory.
From this perspective, the inverse KS problem is conceptually very natural. In the forward theory, the density is determined by the outer minimization and the inner constrained search is carried implicitly inside \(T_\mathrm{s}[\rho]\). In the inverse problem, by contrast, the density is prescribed in advance, and one studies directly the fixed-density constrained search defining \(T_\mathrm{s}[\rho]\) at \(\rho=\rho_{\mathrm{tar}}\). The inverse problem is therefore not an artificial optimization overlay imposed on a nonvariational theory. It is the direct study of the multiplier structure of the inner noninteracting constrained search already present in KS theory itself.
For the purposes of the present work, the constrained-search problem 3 serves as the common prototype underlying the inversion formulations discussed later. We use \(\Phi\) for the noninteracting many-electron state when the constrained-search structure is emphasized, and \(\{\phi_i\}_{i=1}^N\) when the orbital-level KS equations are written explicitly.
The variational picture developed above is conceptually clean, but it relies on nontrivial assumptions. In particular, it presupposes that the constrained search is attained in a sufficiently regular class, that the target density lies in a regime where a meaningful noninteracting variational description is available, and that the associated variational objects are regular enough for the multiplier interpretation to make classical sense [19], [26], [27]. These issues are closely tied to the long-recognized subtleties of representability, differentiability, and density–potential structure in DFT [24], [25].
This subsection is included not as a separate mathematical digression, but because the regularity of \(T_\mathrm{s}\) is precisely what determines how far the clean multiplier picture developed above can be pushed. It therefore provides the conceptual bridge from the regular constrained-search setting of Sec. II to the taxonomy of exact-multiplier, penalty, and explicit state-constraint formulations developed in Sec. III.
From the viewpoint of convex analysis, the noninteracting kinetic-energy functional \(T_\mathrm{s}[\rho]\) is naturally treated as a convex, lower semicontinuous functional on an appropriate density space [19], [28], [29]. In general, however, convexity does not imply differentiability. Even when \(T_\mathrm{s}\) is well defined on a physically relevant domain, it need not admit a unique Fréchet derivative at every density. This distinction is fundamental for inversion: if \(T_\mathrm{s}\) fails to be differentiable at a target density, then the naive identification of the inverse KS potential with a unique classical derivative \(\delta T_\mathrm{s}/\delta \rho\) is no longer justified [25], [26]. This should be distinguished from formal discussions of functional derivatives of noninteracting kinetic-energy density functionals in smoother or more regular settings [30], where the central issue is not the possible failure of differentiability at the target density itself.
A more appropriate generalized replacement is provided, when available, by the subdifferential. More precisely, if one works in a convex-analytic setting in which a suitable extension of \(T_\mathrm{s}\) is treated as a proper lower semicontinuous convex functional on an appropriate density space, then for densities \(\rho_0\) belonging to its effective domain one may consider the subdifferential \[\partial T_\mathrm{s}[\rho_0] = \left\{ \xi\in X^* \;\middle|\; T_\mathrm{s}[\rho]\ge T_\mathrm{s}[\rho_0]+\langle \xi,\rho-\rho_0\rangle \quad \forall \rho\in X \right\}.\]
When \(T_\mathrm{s}\) is Fréchet differentiable at \(\rho_0\), the subdifferential reduces to a singleton containing the usual derivative. At points of nondifferentiability, however, \(\partial T_\mathrm{s}[\rho_0]\) may contain multiple subgradients, reflecting loss of uniqueness in the corresponding dual variable [28], [29]. This does not imply that every admissible target density lies in the effective domain of a suitably defined noninteracting kinetic functional, nor that the subdifferential is nonempty at every density where \(T_\mathrm{s}\) is finite.
In the present setting, the relevance of this construction is direct. The regular multiplier picture suggests that the effective KS potential should be associated, up to sign convention and additive constants, with the variational object dual to the density variable in the constrained search for \(T_\mathrm{s}\). In the differentiable case, this reduces to the familiar local effective potential. More generally, however, whenever the subdifferential is nonempty, the relevant dual variational object is an element of \(\partial T_\mathrm{s}[\rho_0]\) rather than a classical derivative. Accordingly, the subdifferential framework should be viewed not as separate from the multiplier picture, but as its natural convex-analytic extension beyond the fully regular differentiable regime [21], [26], [27].
This is closely tied to representability. Levy’s constrained-search formulation and Lieb’s convex reformulation largely remove the need to assume interacting ground-state \(v\)-representability at the foundational level by working instead with admissible density classes and a variational definition of the universal functional [18], [19]. By contrast, the inverse KS problem reintroduces a distinct noninteracting representability question: whether a prescribed target density can be realized as the ground-state density of some effective local one-body potential. It is this noninteracting \(v\)-representability problem, rather than the original interacting one, that becomes central for the existence and regularity of the inverse KS map [26], [31].
These observations also help explain several familiar numerical pathologies in practical inversion methods, including basis-set artifacts, oscillatory potentials, and instability in weak-gap regimes [3], [11], [13]. In practice, difficulties often arise when the target density is obtained from incomplete basis representations, strongly correlated states, near-degenerate systems, or limiting regimes in which the noninteracting description becomes fragile. In such cases, inversion procedures that implicitly assume a smooth and single-valued density-to-potential map may become unstable, produce strongly oscillatory effective potentials, or fail to converge robustly.
For this reason, differentiability of \(T_\mathrm{s}\) should not be viewed as generic, but as a structural assumption. In favorable regimes, the constrained search gives rise to a well-behaved local multiplier picture, and inverse KS theory may be formulated in terms of an effective potential that is locally unique modulo constants. Outside such regimes, one must instead expect multivalued subgradient structure, loss of smooth duality, or the need for regularized surrogate formulations. The transition from derivatives to subdifferentials therefore marks the point at which inverse KS theory becomes genuinely nonsmooth.
In this sense, representability, nondifferentiability, and the subdifferential structure of \(T_\mathrm{s}\) are not peripheral technicalities. They mark the regularity-level distinctions underlying the three principal realizations discussed in Sec. III: exact multipliers, quadratic penalties, and explicit state constraints.
As discussed above, the inverse KS problem determines the effective potential only modulo an additive constant. The natural mathematical object is therefore not a single potential \(v_\mathrm{s}\), but an equivalence class \[[v_\mathrm{s}]\in X^*/\mathbb{R}.\] This leaves a residual gauge question: once inverse KS theory determines the potential only up to an additive constant, how should the remaining constant be fixed? For finite Coulomb systems, the exact asymptotic structure of the many-body density and of the KS solution selects a distinguished representative from the class \([v_\mathrm{s}]\) and thereby removes the residual gauge freedom [32]–[34].
Within a sufficiently regular noninteracting \(v\)-representable class, one expects the inverse correspondence to be injective modulo constants, in the same Hohenberg–Kohn sense that two distinct potentials differing by more than a constant should not generate the same admissible ground-state density [1], [24], [26]. By contrast, one should not expect global surjectivity on naive ambient spaces such as \(X\) and \(X^*\). Not every admissible density need arise from a sufficiently regular local noninteracting potential, and not every dual-space potential belongs to the image of a well-behaved inverse KS construction [26].
The issue of local regularity is subtler still. Even within representable classes, the dependence of \(v_\mathrm{s}[\rho]\) on \(\rho\) should be viewed as at most locally regular in favorable regimes. Near degeneracies, gap closings, or boundaries of noninteracting representability, differentiability may fail and the inverse correspondence may cease to be single-valued or stable in the classical sense [3], [26], [27]. It is precisely for this reason that the quotient-valued character of the potential and the need for asymptotic normalization become structurally important rather than merely conventional.
For a fixed external potential, the exchange-correlation potential may be viewed pointwise as the derived quantity \[v_{\mathrm{xc}}[\rho](\mathbf{r}) = v_\mathrm{s}[\rho](\mathbf{r})-v_{\mathrm{ext}}(\mathbf{r})-v_\mathrm{H}[\rho](\mathbf{r}).\] Accordingly, at fixed \(v_{\mathrm{ext}}\), whatever regularity or nonuniqueness affects the inverse map \(\rho\mapsto v_\mathrm{s}[\rho]\) is inherited, modulo the explicit Hartree contribution, by the exchange-correlation map \(\rho\mapsto v_{\mathrm{xc}}[\rho]\).
For finite Coulomb systems, the residual additive ambiguity may be removed by imposing the exact asymptotic structure of the KS solution. Consider a finite \(N\)-electron Coulomb system with total nuclear charge \(Z\), external potential satisfying \(v_{\mathrm{ext}}(r)\sim -Z/r\), and exact ground-state density \(\rho(r)\). Assume that the system is bound with first ionization energy \(I>0\) and that the exact density is noninteracting \(v\)-representable by a KS potential \(v_\mathrm{s}[\rho]\), unique up to an additive constant. In the asymptotic region, away from nodal directions, the exact many-body density has the form \[\rho(r)\sim C\,r^{2\beta}e^{-2\kappa r}, \qquad \kappa=\sqrt{2I}, \qquad \beta=\frac{Z-N+1}{\kappa}-1\] [32], [34].
On the other hand, if the KS potential has the asymptotic form \(v_\mathrm{s}[\rho](r)\sim -\gamma/r\), then the associated highest occupied orbital obeys \[\varepsilon_{\mathrm{HOMO}}=-\frac{\kappa^2}{2}, \qquad \beta=\frac{\gamma}{\kappa}-1.\] Using the ionization-potential theorem \(\varepsilon_{\mathrm{HOMO}}=-I\) [34], [35] and comparing the two expressions for \(\beta\), one obtains \[v_\mathrm{s}[\rho](r)\sim -\frac{Z-N+1}{r}, \qquad r\to\infty.\] Since \(v_{\mathrm{ext}}(r)\sim -Z/r\) and \(v_\mathrm{H}[\rho](r)\sim N/r\), it follows that \[v_{\mathrm{xc}}[\rho](r)=-\frac{1}{r}+\mathcal{O}(r^{-2}), \qquad r\to\infty.\]
This asymptotic analysis has a direct gauge-theoretic interpretation. In a regular noninteracting \(v\)-representable regime, the inverse KS problem determines only the equivalence class \([v_\mathrm{s}]\in X^*/\mathbb{R}\). For finite Coulomb systems, however, the exact large-\(r\) asymptotic condition selects a distinguished representative from that class. In this sense, asymptotics does not merely provide additional physical information; it resolves the residual gauge indeterminacy of the inverse map and promotes a quotient-valued density-to-potential correspondence to a normalized representative [32], [33].
The preceding discussion may be reorganized in direct analogy with Lieb’s convex formulation of DFT [19]. In the interacting case, one passes from the constrained-search functional of Levy and Lieb to the convex-analytic density functional naturally paired with the ground-state energy through Fenchel duality. More precisely, if \[E[v] = \inf_{\rho\in X} \left\{ F[\rho]+\langle v,\rho\rangle \right\},\] then the corresponding interacting density functional may be written formally as \[F[\rho] = \sup_{v\in X^*} \left\{ E[v]-\langle v,\rho\rangle \right\}.\] Here \(F\) denotes the interacting Lieb functional in the convex-analytic sense, namely the lower-semicontinuous convex functional paired with the ground-state energy through Fenchel duality [19]. It should be distinguished from the Levy–Lieb constrained-search functional \(F_{\mathrm{LL}}\) introduced earlier: while \(F_{\mathrm{LL}}\) is defined directly by constrained minimization over interacting wavefunctions, the Lieb functional \(F\) is the lower-semicontinuous convex envelope, equivalently the biconjugate, associated with the dual ground-state energy formulation [19], [28], [36].
On the noninteracting side, one may introduce the noninteracting ground-state energy \[E_\mathrm{s}[v] = \inf_{\Phi} \left\{ \langle \Phi|\hat{T}|\Phi\rangle + \langle v,\rho_\Phi\rangle \right\}, \label{eq:Es95def}\tag{4}\] where \(v\) lies in a suitable potential space dual to the density space. Formally, this yields the density-side variational principle \[E_\mathrm{s}[v] = \inf_{\rho\in X} \left\{ T_\mathrm{s}[\rho]+\langle v,\rho\rangle \right\}, \label{eq:Es95dual}\tag{5}\] with \(T_\mathrm{s}\) playing the role of the relevant noninteracting density functional. Conversely, in the same convex-analytic spirit, \[T_\mathrm{s}[\rho] = \sup_{v\in X^*} \left\{ E_\mathrm{s}[v]-\langle v,\rho\rangle \right\}, \label{eq:Ts95dual}\tag{6}\] up to the usual sign convention for the conjugate pairing [19], [28], [29]. In this sense, the noninteracting kinetic functional admits a Lieb-type dual interpretation closely parallel to that of the interacting functional.
This dual perspective clarifies the mathematical status of the density–potential mapping. In a sufficiently regular, nondegenerate regime, one may use the schematic notation \(\mathcal{M}:v\mapsto \rho[v]\) and its formal inverse as convenient shorthand. More generally, however, the relevant density–potential relations should be understood through subdifferentials on the density side and superdifferentials, or equivalently subdifferentials of the negated energy functionals, on the potential side [26], [28], [29]. When the relevant functional is differentiable at a given point, the associated differential object reduces to a singleton, and the corresponding density–potential relation is locally single-valued there. At points of nondifferentiability, by contrast, the relevant differential object is naturally set-valued, and the corresponding relation must be understood in a multivalued sense. A distinct failure mode occurs when the density or potential under consideration does not belong to the effective domain of the corresponding functional at all. In that case, one does not obtain a multivalued relation, but rather no density–potential relation is defined there within the given variational framework. Thus the schematic notation \(\mathcal{M}\) and its formal inverse are best regarded as shorthand valid only in sufficiently regular, nondegenerate, and representable regimes.
A further structural issue is that the interacting and noninteracting dual frameworks need not have identical effective domains and ranges. A density that is admissible on the interacting side need not automatically belong to the density class compatible with the noninteracting dual structure. From this viewpoint, the noninteracting \(v\)-representability problem may be understood as a compatibility problem between the density class accessible in exact DFT and the density class admitting realization within the noninteracting density–potential duality associated with \(T_\mathrm{s}\) [19], [26]. This should not, however, be read as a claim that the full KS decomposition inherits the same global convex structure as the Lieb functional framework. The relevant convex-analytic structure attaches to the noninteracting component \(T_\mathrm{s}\), while the Hartree term is explicit and the exchange-correlation contribution is not in general expected to define a convex functional. The point is therefore not that full KS theory reduces to a globally clean convex program, but rather that its noninteracting component inherits a Lieb-type dual structure that provides the appropriate variational setting for the inverse problem.
The purpose of this section is to show that several major inversion formulations may be understood as distinct variational realizations of the same inverse Kohn–Sham problem. The variational anchor developed in Sec. II is the density-constrained noninteracting constrained search associated with \(T_\mathrm{s}[\rho]\). That anchor provides the natural starting point for interpreting inverse KS theory within exact DFT. At the same time, the final part of Sec. II showed that this noninteracting constrained-search structure may also be reorganized in a Lieb-type density–potential dual perspective, with \(T_\mathrm{s}\) and \(E_\mathrm{s}[v]\) playing the roles of dual noninteracting density and potential functionals. The broader classification developed in the present section is not limited to direct optimization reformulations of \(T_\mathrm{s}\) itself. Rather, it rests on the observation that inverse KS inversion is governed by two defining relations: the KS state equations and the density-reproduction condition. What distinguishes the principal formulations is how these relations are organized within the optimization architecture.
From this viewpoint, three broad structural classes can already be distinguished. In exact-multiplier formulations, the density condition is enforced sharply and the noninteracting state variables are reduced out. In quadratic-penalty formulations, exact density reproduction is relaxed and the mismatch is penalized. In explicit state-constraint formulations, the KS equations are retained as visible constraints. These distinctions form the basis of the structural classification developed below.
For notational convenience, when a symmetric positive kernel \(K\) is used, we define \[\|f\|_K^2 := \iint_{\Omega\times\Omega} f(\mathbf{r})\,K(\mathbf{r},\mathbf{r}')\,f(\mathbf{r}')\, d\mathbf{r}\,d\mathbf{r}'.\] In orbital-based formulations, we also write \[D[\{\phi_i\}] := \sum_{i=1}^N |\phi_i|^2-\rho_{\mathrm{tar}}, \qquad R_i[\phi_i,v,\varepsilon_i] := \left(-\frac{1}{2}\nabla^2+v-\varepsilon_i\right)\phi_i.\] Unless otherwise specified, norms on state and density residuals are understood schematically and need not be the same.
Operationally, the Wu–Yang method may be understood as a formulation in which the density-reproduction condition 2 is enforced exactly, while the KS state equations 1 are treated implicitly through their dependence on the trial effective potential. Variationally, this is the exact-multiplier treatment of the fixed-density constrained search introduced in Sec. II, followed by elimination of the noninteracting state variables [3], [4], [9], [10]. In the language of the Lieb-type dual perspective developed at the end of Sec. II, Wu–Yang may be viewed as the reduced-potential realization of the noninteracting density–potential dual structure associated with \(T_\mathrm{s}\).
Introducing a multiplier \(v(\mathbf{r})\) for the condition \(\rho_\Phi=\rho_{\mathrm{tar}}\) leads formally to the constrained Lagrangian \[\mathcal{L}(\Phi,v) = \langle \Phi|\hat{T}|\Phi\rangle + \langle v,\rho_\Phi-\rho_{\mathrm{tar}}\rangle,\] supplemented by the orbital orthonormality constraints. This already has the structure of a saddle problem, \[\inf_{\Phi}\sup_v \mathcal{L}(\Phi,v).\]
The characteristic Wu–Yang reduction is obtained by eliminating the state variables. Fix a multiplier field \(v\in X^*\) and minimize the Lagrangian over noninteracting states subject to orthonormality. Equivalently, solve the KS state equations for the current trial potential and absorb the state variables into a reduced potential functional. This yields \[\begin{align} W[v] :={}& \inf_{\Phi} \left\{ \langle \Phi|\hat{T}|\Phi\rangle + \langle v,\rho_\Phi-\rho_{\mathrm{tar}}\rangle \right\} \\ ={}& E_\mathrm{s}[v]-\langle v,\rho_{\mathrm{tar}}\rangle , \end{align}\] where \[E_\mathrm{s}[v] = \inf_\Phi \left\{ \langle \Phi|\hat{T}|\Phi\rangle + \langle v,\rho_\Phi\rangle \right\}\] denotes the noninteracting ground-state energy in the potential \(v\).
This reduced functional is exactly what one expects from the noninteracting Lieb-type dual perspective of Sec. II. There, \(E_\mathrm{s}[v]\) was identified as the potential-side object dual to the density functional \(T_\mathrm{s}[\rho]\). From this viewpoint, the Wu–Yang functional is simply the reduced exact-multiplier functional obtained by pairing the noninteracting dual energy with the prescribed target density. In this sense, the Wu–Yang formulation is not merely an ad hoc reduction to potential space, but the most direct reduced-potential expression of the density-centered \(T_\mathrm{s}\) framework developed in Sec. II.
At this point it is useful to distinguish between the reduced-potential viewpoint and the stronger language of exact duality. The reduced functional \(W[v]\) is obtained directly from exact multiplier enforcement followed by elimination of the state variables, and this already suffices to explain the operational content of the Wu–Yang algorithm. However, the reduced problem is not automatically identical to the original constrained problem. In general, one only has the weak-duality relation \[\sup_v \inf_\Phi \mathcal{L}(\Phi,v) \;\le\; \inf_\Phi \sup_v \mathcal{L}(\Phi,v).\] Thus the reduced maximization over \(v\) should not be identified automatically with the exact density-constrained inverse problem without additional assumptions ensuring an appropriate saddle structure, dual attainment, or strong-duality property [19], [28], [29].
Assuming sufficient regularity of the noninteracting ground-state energy with respect to the potential, one has \[\frac{\delta E_\mathrm{s}[v]}{\delta v(\mathbf{r})}=\rho_v(\mathbf{r}),\] where \(\rho_v\) is the noninteracting ground-state density associated with \(v\). Hence \[\frac{\delta W[v]}{\delta v(\mathbf{r})} = \rho_v(\mathbf{r})-\rho_{\mathrm{tar}}(\mathbf{r}).\] A stationary point of \(W[v]\) therefore corresponds formally to exact satisfaction of the density-reproduction condition.
Under nondegeneracy and suitable spectral-gap assumptions, the second variation of \(W[v]\) is governed by the static noninteracting density-response operator, \[\frac{\delta^2 W}{\delta v(\mathbf{r})\,\delta v(\mathbf{r}')} = \chi_\mathrm{s}(\mathbf{r},\mathbf{r}'),\] where, in a standard orbital representation, \[\chi_\mathrm{s}(\mathbf{r},\mathbf{r}') = 2 \sum_{i\in \mathrm{occ}} \sum_{a\in \mathrm{unocc}} \frac{ \phi_i^*(\mathbf{r})\phi_a(\mathbf{r})\phi_a^*(\mathbf{r}')\phi_i(\mathbf{r}') }{ \varepsilon_i-\varepsilon_a }.\] This makes the local optimization structure transparent: Wu–Yang may be interpreted as a Newton or quasi-Newton method in potential space, with gradient given by density mismatch and Hessian given by the KS response operator. It also explains why the method becomes fragile in small-gap or nearly degenerate regimes, where the same response operator becomes poorly conditioned [3], [4], [13], [14].
Structural conclusion. The Wu–Yang formulation is the reduced exact-multiplier realization of the common density-constrained inverse Kohn–Sham problem.
A natural alternative to exact multiplier enforcement is to relax the density-reproduction condition through a quadratic penalty term. Operationally, the KS state equations remain in place as the forward state relation associated with the current effective potential, but exact density matching is no longer imposed as a hard constraint. Instead, density mismatch is assigned a quadratic cost in the objective. From the present viewpoint, this is the essential variational structure of the Zhao–Morrison–Parr method [8]. In other words, ZMP is best understood as a quadratic-penalty relaxation of the density-constrained inverse KS problem [3], [8], [21].
At the schematic level, one replaces the exact constrained problem by a penalized functional of the form \[\mathcal{J}_\lambda[\Phi] = \langle \Phi|\hat{T}|\Phi\rangle + \frac{\lambda}{2}\|\rho_\Phi-\rho_{\mathrm{tar}}\|^2,\] or, more generally, by a kernel-weighted quadratic penalty \[\mathcal{J}_\lambda[\Phi] = \langle \Phi|\hat{T}|\Phi\rangle + \frac{\lambda}{2}\|\rho_\Phi-\rho_{\mathrm{tar}}\|_K^2.\] Here \(K\) is a symmetric positive kernel and \(\lambda>0\) is the penalty parameter. In the original ZMP construction, the Coulomb kernel is the canonical choice [8], but the structural point is more general: the hard density constraint has been replaced by a quadratic penalty that drives the solution toward the target density as \(\lambda\) increases.
Thus ZMP does not solve the hard-constrained inverse problem exactly at finite \(\lambda\); rather, it solves a penalized surrogate problem. Accordingly, the potential obtained at finite \(\lambda\) should be viewed as a penalty-surrogate potential, not by itself as the exact inverse KS potential. For feasible states satisfying \(\rho_\Phi=\rho_{\mathrm{tar}}\), the penalty term vanishes, so the penalized objective agrees with the exact constrained objective on the constraint manifold. Away from that manifold, however, the minimizer may trade density mismatch against kinetic optimality. ZMP should therefore be understood as a family of approximate inverse problems that approach the exact density-constrained problem only in the large-penalty limit.
The corresponding effective potential is obtained by taking the first variation of the penalty term with respect to the density. One is then led to the density-dependent correction \[v_\lambda(\mathbf{r}) = \lambda\,(K*(\rho_\Phi-\rho_{\mathrm{tar}}))(\mathbf{r}).\] Thus the density mismatch is fed back directly into the one-body potential through the kernel \(K\). This is the characteristic operational content of ZMP: rather than introducing an independent multiplier variable and solving a reduced optimization problem over potentials, one iteratively corrects the effective potential by penalizing the current density residual [6], [8].
The practical attraction of the penalty viewpoint is robustness. Even when the target density lies near the boundary of regular noninteracting \(v\)-representability, the penalized problem remains meaningful as an approximate inversion task. One no longer asks for exact feasibility at every stage, but for the best compromise between noninteracting kinetic optimality and density matching at finite \(\lambda\). This frequently makes the optimization more forgiving than an exact reduced-multiplier or reduced-dual formulation [3], [13], [14].
At the same time, the weakness of quadratic-penalty methods is classical. If one wants to recover the exact density constraint, the penalty parameter must be taken large. As \(\lambda\) increases, however, the geometry of the optimization landscape becomes increasingly anisotropic. At a schematic second-order level, one may write \[\nabla^2 \mathcal{J}_\lambda \approx \nabla^2 T_\mathrm{s} + \lambda (\nabla C)^*(\nabla C), \qquad C[\Phi]=\rho_\Phi-\rho_{\mathrm{tar}},\] which makes clear that density-changing directions acquire stiffness proportional to \(\lambda\), while tangential directions do not [28], [29].
In practical implementations, one often increases the penalty parameter progressively. The resulting family of penalized problems therefore admits a continuation, or homotopy-like, interpretation: one follows a path from relatively soft density matching toward increasingly accurate enforcement of the target density as \(\lambda\) grows. More abstractly, the ZMP construction lies close in spirit to Moreau–Yosida-type regularization [21]. The hard density constraint is replaced by a smoother family of penalized problems in which exact feasibility is traded for improved local regularity of the objective landscape. The precise extent to which a given ZMP formulation coincides with a classical Moreau envelope depends on the function-space setting and metric structure, but this convex-analytic interpretation is best understood as a posterior explanation of the quadratic-penalty structure rather than as the primary definition of the method.
Like Wu–Yang, this formulation remains closely tied to the density-centered constrained-search picture, but replaces exact multiplier enforcement by quadratic penalization.
Structural conclusion. The ZMP formulation is the quadratic-penalty relaxation of the same density-constrained inverse Kohn–Sham problem.
A third route is to enlarge the constrained system rather than to eliminate the state variables through a reduced potential optimization or to relax the density condition through a quadratic penalty alone. In this approach, the KS state equations 1 are retained explicitly as constraints, while the density-reproduction condition 2 is enforced through a tracking-type objective or, in some formulations, through additional explicit constraints. From the present viewpoint, this is the essential structure of the PDE-constrained paradigm: the inverse KS problem is formulated as an explicit state-constraint optimization problem rather than as a reduced multiplier problem or a pure penalty relaxation [3], [4], [22], [23].
At the schematic level, one considers an optimization problem of the form \[\min_{\{\phi_i\},\,v}\; \mathcal{J}[\{\phi_i\}]\] subject to the KS state equations 1 and the usual orbital orthonormality conditions. Here \(\{\phi_i\}_{i=1}^N\) are the state variables, \(v\) is the control variable, and \(\mathcal{J}\) is typically a density-mismatch functional such as \[\mathcal{J}[\{\phi_i\}] = \frac{1}{2} \|\rho-\rho_{\mathrm{tar}}\|_K^2, \qquad \rho=\sum_{i=1}^N |\phi_i|^2.\] The key structural feature is that the KS equations remain explicitly visible as part of the constrained variational system.
Unlike Wu–Yang, this approach does not eliminate the state variables in favor of a potential-only optimization. Unlike ZMP, it does not replace the hard state relation by a finite-penalty surrogate. Instead, it keeps the relevant state equations visible throughout the optimization.
To derive first-order optimality conditions, one introduces adjoint variables \(\{p_i\}_{i=1}^N\) and forms a control Lagrangian \[\begin{align} \mathcal{L}(\{\phi_i\},v,\{\varepsilon_i\},\{p_i\},\Lambda) :={}& \mathcal{J}[\{\phi_i\}] \\ &+ \sum_{i=1}^N \left\langle p_i,\, R_i[\phi_i,v,\varepsilon_i] \right\rangle \\ &+ \sum_{i,j=1}^N \Lambda_{ij}\bigl(\langle \phi_i,\phi_j\rangle-\delta_{ij}\bigr). \end{align}\] where \(\Lambda\) enforces orbital orthonormality. The resulting stationarity conditions couple the state equations, adjoint equations, normalization conditions, and control gradient relation into a single optimality system [22], [23]. For notational economy, we reuse the symbol \(\mathcal{L}\) for the relevant Lagrangian in each formulation.
The adjoint equations arise from variation with respect to the orbitals and have the characteristic form of projected Sternheimer or linear-response equations. Schematically, \[\left( -\frac{1}{2}\nabla^2+v(\mathbf{r})-\varepsilon_k \right) p_k(\mathbf{r}) = -\frac{\delta \mathcal{J}}{\delta \phi_k}(\mathbf{r}) - \sum_{j=1}^N \Lambda_{kj}\phi_j(\mathbf{r}),\] to be understood in the appropriate weak or projected sense. The essential computational point is that, once the adjoint system is solved, the gradient of the objective with respect to the control variable may be obtained without explicitly differentiating the orbitals with respect to the potential [4], [22], [23].
The resulting control gradient is obtained by variation with respect to the potential. Formally, it has the schematic structure \[\frac{\delta \mathcal{L}}{\delta v(\mathbf{r})} \sim \sum_{k=1}^N \phi_k(\mathbf{r})p_k(\mathbf{r}),\] up to the precise pairing convention and weak-form normalization being used. Thus the inversion may be attacked by gradient-based optimization in the control variable, with derivative information supplied by state–adjoint pairs.
At the same time, the weakness of the explicit state-constraint paradigm is equally clear. The projected adjoint equations involve shifted KS operators whose inverses are controlled by spectral gaps. In well-gapped systems, the reduced operator behaves stably and the adjoint equations provide a reliable route to accurate gradients. In weak-gap, nearly degenerate, or metallic regimes, however, the same projected operators become poorly conditioned or nearly singular. The resulting loss of robustness is therefore not accidental; it reflects the spectral structure of the underlying linearized inverse problem itself [3], [4], [13], [14].
Thus the PDE-constrained formulation is most naturally understood at the broader structural level, where the KS state equations are retained explicitly rather than absorbed into a \(T_s\)-centered reduced formulation.
Structural conclusion. PDE-constrained inversion is the explicit state-constraint realization of the same inverse Kohn–Sham problem at the orbital level.
Response-based and stabilized inversion schemes can often be interpreted as algorithmic realizations, regularized variants, or iterative solvers associated with the same underlying structural classes identified here, rather than as fundamentally separate variational or optimization-theoretic classes.
Having identified the three principal realizations in Sec. III, we now abstract the optimization pattern behind them. The classification is governed by the two defining inverse-KS relations: the KS state equations and the density-reproduction condition. Different formulations arise according to whether these relations are treated as objectives, constraints, penalties, or feasibility relations.
From this viewpoint, the formulations of Sec. III are not merely a list of representative inversion schemes. They are specific realizations of a more general taxonomy of inverse KS optimization problems. This taxonomy is defined at the level of optimization structure rather than at the level of existing algorithmic implementations or standard inversion workflows. One may formulate the problem with one relation entering as an optimization objective and the other as an explicit constraint; one may promote both relations to residuals or tracking terms inside a joint objective; or one may treat both as coupled constraints, leading primarily to a feasibility or equation-solving viewpoint rather than to a standard optimization problem. [4], [22], [23]
In this sense, Sec. III identifies principal structural realizations, whereas the present section makes explicit the broader optimization-theoretic taxonomy into which those realizations fit.
The first and most important class consists of formulations in which one of the two defining inverse-KS relations is assigned the role of an optimization objective, while the other is retained as an explicit constraint. This class already contains the main paradigms discussed in Sec. III.
One possibility is to retain the density relation as the central constraint and to absorb the state relation into the optimization architecture. This is the viewpoint most closely connected to the fixed-density constrained search of Sec. II. The Wu–Yang formulation belongs to this branch: the density condition remains the essential target, while the KS state equations are treated implicitly through reduction to a potential-space problem [9]. The ZMP construction belongs to the same branch, except that the density condition is no longer enforced sharply and is instead softened into a quadratic penalty [8], [21].
The complementary possibility is to retain the KS state equations as explicit constraints and to move density reproduction into the objective. This is the natural setting of PDE-constrained inversion. There the state equations remain visible throughout the formulation, while density mismatch is measured by a tracking functional over the coupled state–control variables [4], [22], [23].
Within this same class, augmented-Lagrangian ideas arise naturally as enriched versions of one-objective/one-constraint formulations. Their role is not to define a separate inversion class, but to strengthen formulations in which one relation is retained as a hard constraint by combining multiplier information with moderate regularization [28], [29].
Consider the density-centered constrained-search viewpoint of Sec. II, \[\min_{\Phi}\; \langle \Phi|\hat{T}|\Phi\rangle \qquad \text{subject to } \rho_\Phi=\rho_{\mathrm{tar}}.\] An augmented-Lagrangian enrichment of this formulation takes the form \[\mathcal{L}[\Phi,\mu] = \langle \Phi|\hat{T}|\Phi\rangle + \langle \mu,\rho_\Phi-\rho_{\mathrm{tar}}\rangle + \frac{\lambda}{2}\|\rho_\Phi-\rho_{\mathrm{tar}}\|^2,\] where \(\mu\) is a density multiplier and \(\lambda>0\) is a finite penalty parameter. If the quadratic term is removed, one recovers the exact-multiplier viewpoint underlying Wu–Yang; if the multiplier is frozen or omitted, one recovers a ZMP-type penalty structure. In this sense, the augmented-Lagrangian construction is a natural extension within the same one-objective/one-constraint class.
Relative to a pure exact-multiplier formulation, the quadratic term provides additional local stabilization and may improve robustness of the constrained search. Relative to a pure penalty formulation, the multiplier absorbs part of the burden of feasibility, so that exact density reproduction need not be approached solely through a numerically dangerous large-penalty limit. At the level of optimization structure, this viewpoint yields a natural primal–dual coupled formulation; algorithmically, it suggests iterations of Uzawa type, for example, \[\Phi^{(k+1)} \approx \arg\min_{\Phi}\, \mathcal{L}\big[\Phi,\mu^{(k)}\big],\] followed by \[\mu^{(k+1)}(\mathbf{r}) = \mu^{(k)}(\mathbf{r}) + \lambda\bigl(\rho_{\Phi^{(k+1)}}(\mathbf{r})-\rho_{\mathrm{tar}}(\mathbf{r})\bigr).\]
A second broad class arises when neither of the two defining inverse-KS relations is singled out as the unique retained outer constraint. Instead, both relations are promoted to residuals, penalties, or tracking terms inside a joint objective. This leads to enlarged all-at-once formulations in which state consistency and density reproduction are treated simultaneously at the objective level.
Writing the occupied orbitals as \(\{\phi_i\}_{i=1}^N\), one natural example is \[\min_{\{\phi_i\},\,v,\,\{\varepsilon_i\}} \; w_1\sum_{i=1}^N \|R_i[\phi_i,v,\varepsilon_i]\|^2 + w_2\|D[\{\phi_i\}]\|^2,\] together with the usual orbital orthonormality constraints and, where desired, additional regularity penalties on the potential or gauge-fixing terms. Here the first term measures violation of the KS state equations, while the second measures violation of density reproduction. In such a formulation, neither relation is treated solely as a hidden solver nor solely as a hard constraint. Instead, both enter symmetrically as coupled optimization targets.
This enlarged formulation unifies structures that otherwise appear separately in reduced exact-multiplier formulations, quadratic-penalty relaxations, and explicit state-constraint methods. At the same time, it introduces nontrivial design choices: one must decide how the two residual blocks should be balanced, whether through fixed weights \(w_1,w_2\), continuation in one or both weights, or adaptive merit strategies that respond to the evolving scale of the state and density residuals.
At the algorithmic level, this opens a substantially richer design space than the classical nested inversion picture. One could imagine alternating minimization in \((\{\phi_i\},v,\{\varepsilon_i\})\), block-coordinate descent, sequential quadratic programming, Newton–Krylov methods on the coupled residual system, or hybrid schemes in which adjoint-based derivative information is combined with multiplier updates and continuation in regularization parameters [22], [23]. Likewise, modern parameterized representations of orbitals and potentials, such as finite-dimensional basis expansions, symmetry-adapted ansätze, and neural-network or PINN-style parameterizations, may permit genuinely joint optimization strategies that are difficult to formulate naturally within rigidly asymmetric inversion formulations and schemes [37].
At this level of abstraction, the issue is not whether a particular member of this class has already become standard in the inversion literature, but that assigning objective-type roles to both defining inverse-KS relations defines a distinct optimization class in its own right.
Structural conclusion. Enlarged all-at-once residual formulations constitute the two-objective class in the optimization-theoretic taxonomy of inverse Kohn–Sham inversion.
A third limiting class is obtained when both defining inverse-KS relations are treated as constraints rather than as objectives. In this case, the inverse problem is viewed primarily as a coupled nonlinear system to be solved rather than as a conventional optimization problem with a distinguished cost functional. One then seeks orbitals, a potential, and orbital energies satisfying both the KS state equations 1 and the density-reproduction condition 2 simultaneously, together with the usual orbital orthonormality conditions and, where needed, a gauge-fixing condition for the potential.
This viewpoint may be framed formally as the feasibility problem \[\text{find }(\{\phi_i\},v) \quad\text{subject to}\quad \eqref{eq:ks95state} \;\text{and}\; \eqref{eq:density95reproduction},\] together with the usual orbital orthonormality conditions. The goal is not to minimize a preferred objective while enforcing one dominant constraint, but to find a jointly consistent solution of the coupled state and density relations.
Even in the absence of a primary outer objective, it is still natural to write down a constrained Lagrangian. Using the same orbital-level notation as above, let \(\{p_i\}_{i=1}^N\) denote multipliers associated with the KS state equations, let \(\mu(\mathbf{r})\) denote the multiplier associated with density reproduction, and let \(\Lambda=(\Lambda_{ij})\) enforce orbital orthonormality, while the orbital energies \(\{\varepsilon_i\}_{i=1}^N\) are treated as auxiliary unknowns of the explicit state-equation block. One is then led schematically to \[\begin{align} \mathcal{L}(\{\phi_i\},v,\{\varepsilon_i\},\{p_i\},\mu,\Lambda) :={}& \langle \mu, D[\{\phi_i\}] \rangle \\ &+ \sum_{i=1}^N \langle p_i, R_i[\phi_i,v,\varepsilon_i]\rangle \\ &+ \sum_{i,j=1}^N \Lambda_{ij}\bigl(\langle\phi_i,\phi_j\rangle-\delta_{ij}\bigr). \end{align}\] In this way, one obtains a KKT-like or saddle-system viewpoint in which both defining inverse relations appear explicitly as constraints.
This zero-objective viewpoint is the limiting feasibility realization of the same inverse-KS structure. Its role in the present work is not merely supplementary: it completes the optimization-theoretic classification by showing what remains when both defining relations are treated as constraints rather than as objectives or penalties.
Structural conclusion. The zero-objective feasibility viewpoint constitutes the two-constraint class in the optimization-theoretic taxonomy of inverse Kohn–Sham inversion.
We have developed a unified mathematical framework for the inverse Kohn–Sham problem by proceeding in two steps. First, we identified the density-constrained noninteracting constrained search associated with exact density functional theory as the natural variational anchor of inverse KS inversion. Second, building on this anchor, we developed an optimization-theoretic classification of inverse KS formulations according to how the KS state equations and density-reproduction condition are organized as objectives, constraints, penalties, or feasibility relations.
Within this framework, Wu–Yang appears as a reduced exact-multiplier realization, ZMP as a quadratic-penalty relaxation of the same density-constrained problem, and PDE-constrained approaches as explicit state-constraint formulations at the orbital level. This common viewpoint helps clarify why additive-constant ambiguity, asymptotic normalization, nondifferentiability, loss of dual attainment, and weak-gap instability recur across inversion methods: they reflect different manifestations of the same underlying variational and structural content rather than unrelated algorithmic pathologies [3], [21], [27]. The point of this classification is therefore structural: it organizes inverse KS theory at the level of optimization architecture rather than at the level of particular established algorithms.
The framework also naturally accommodates broader formulation classes. In particular, augmented-Lagrangian formulations provide an intermediate route between exact multiplier enforcement and penalty regularization, while enlarged all-at-once formulations offer a setting in which state consistency and density reproduction may be treated on more nearly equal footing. These extensions broaden the optimization-theoretic picture developed here, although their detailed analytical properties and numerical conditioning still merit systematic study.
More broadly, the inverse KS problem may be viewed not only as a practical tool for extracting reference exchange-correlation potentials, but also as a structured model problem at the intersection of density functional theory, convex analysis, nonsmooth optimization, and PDE-constrained control [19], [22], [23], [28], [29]. We hope that the present formulation helps sharpen these connections and provides a clearer language for comparing existing inversion formulations and methods and for developing new ones.