# Solution of the $v$-representability problem on a ring domain

Sarina M. Sutter  Theoretical Chemistry, Faculty of Exact Sciences, VU
University, Amsterdam, The Netherlands  Markus Penz<EMAIL_ADDRESS>Max
Planck Institute for the Structure and Dynamics of Matter and Center for Free-
Electron Laser Science, Hamburg, Germany  Department of Computer Science, Oslo
Metropolitan University, Oslo, Norway  Michael Ruggenthaler  Max Planck
Institute for the Structure and Dynamics of Matter and Center for Free-
Electron Laser Science, Hamburg, Germany  The Hamburg Center for Ultrafast
Imaging, Hamburg, Germany  Robert van Leeuwen  Department of Physics,
Nanoscience Center, University of Jyväskylä, Finland  Klaas J. H. Giesbertz
Theoretical Chemistry, Faculty of Exact Sciences, VU University, Amsterdam,
The Netherlands

###### Abstract

We provide a solution to the $v$-representability problem for a non-
relativistic quantum many-particle system on a ring domain in terms of Sobolev
spaces and their duals. Any one-particle density that is square-integrable,
has a square-integrable weak derivative, and is gapped away from zero can be
realized from the solution of a many-particle Schrödinger equation, with or
without interactions, by choosing a corresponding external potential. This
potential can contain a distributional contribution but still gives rise to a
self-adjoint Hamiltonian. Importantly, this allows for a well-defined
Kohn–Sham procedure but, on the other hand, invalidates the usual proof of the
Hohenberg–Kohn theorem.

###### Contents

  1. I Introduction
  2. II One-dimensional torus setting
  3. III Densities on the one-dimensional torus
  4. IV Convex formulation of density-functional theory
  5. V $v$-Representing potentials from non-empty subdifferentials
  6. VI KLMN theorem for distributional potentials
  7. VII Conclusions and outlook
  8. A Kinetic bounds for interactions

##  I Introduction

The origin of the $v$-representability problem can be found in the seminal
work of Hohenberg and Kohn [1, Footnote 12] that established the field of
density-functional theory (DFT). Therein, the universal density functional is
only defined for densities with the property of $v$-representability. The term
“$v$-representability” itself, according to Levy [2], is due to E. G. Larson
at the Boulder Theoretical Chemistry Conference, June 1975. The same paper
also presents the constrained-search approach as a way to extend the domain of
the universal functional to all $N$-representable densities (those that come
from $N$-particle wave functions). It was later realized by Levy [3] and Lieb
[4] that there exist densities that belong to a ground-state _ensemble_ of at
least $3$-fold degeneracy, but not to any pure ground state. Also, it was
found by Englisch and Englisch [5] that even this more general ensemble
$v$-representability is not enough to generate all possible densities. Since
it is also imaginable that Hamiltonians with and without interactions yield
entirely different ground-state densities, we arrive at the following types of
$v$-representability.

non-interacting pure-state $v$-representability | non-interacting ensemble $v$-representability
---|---
interacting pure-state $v$-representability | interacting ensemble $v$-representability

This distinction is indeed important when it comes to the widely-applied
method of Kohn and Sham [6]. The Kohn–Sham approach, where a non-interacting
auxiliary system is employed to determine the density of an interacting
physical system, relies on the assumption that the density is simultaneously
interacting and non-interacting $v$-representable. And so the original
$v$-representability problem is not solved just by the definition of
constrained-search density functionals. Accordingly, the $v$-representability
problem is one of the major open problems in the mathematical formulation of
density-functional theory [7; 8]. A special focus was given to the issue of
$v$-representability also by Mathieu Lewin in Teale  _et al._ [9, 4.5.2].

In this work, we will always consider the ensemble case and at the same time
allow for a whole class of different interaction operators $W$ in the
Hamiltonian. Since this class includes $W=0$, the difference between
interacting and non-interacting $v$-representability will not be relevant in
our treatment and the presented solution applies simultaneously to both cases.
The same general approach was followed in Chayes  _et al._ [10], where
$v$-representability for all fermionic densities $0<\rho(x)<1$ on an
(infinite) lattice was demonstrated. A simplification of the proof for finite
lattices can be found in Penz and van Leeuwen [11]. Lammert [12] was able to
provide $v$-representability for a coarse-grained formulation of DFT, and he
also provided crucial insights into the continuum case [13]. This work
demonstrated that the universal density functional, as formalized in Lieb [4],
is not functionally differentiable. Since functional differentiability would
imply $v$-representability, the counter-examples provided by Lammert [13]
highlight mathematical obstacles that have to be overcome in order to
guarantee that a density is indeed $v$-representable. In order to guarantee
functional differentiability of a density $\rho(x)>0$, he suggested to work
with the Sobolev space $H^{2}$ for densities instead of the Lebesgue spaces
employed by Lieb [4]. Sobolev spaces will play a crucial role in this work.

Apart from this paper, notable progress has been achieved by Aryasetiawan and
Stott [14, 15] and Aryasetiawan [16], where $v$-representability is studied
for a _non-interacting_ system on a one-dimensional bounded interval with zero
boundary conditions. In this setting, $v$-representability can be achieved for
any density that is strictly positive except at the boundary and has a finite
derivative everywhere. A key ingredient is the transformation of the decoupled
$N$-particle Schrödinger equation into $(N-1)$ coupled and non-linear
differential equations that depend on the given density $\rho$ instead of the
external potential $v$. While the issue of making sure that a solution also
gives back $\rho$ as a ground-state density that appears in higher dimensions
has been tackled in subsequent work by Chen and Stott [17, 18, 19], the issue
of solvability of the coupled non-linear differential equations was left open.
A recent work of Garrigue [20] approaches the issue of $v$-representability as
an inverse problem with discretized potential space and also discusses the
possibility of excited states. Nevertheless, the original $v$-representability
problem, especially in the case of interacting particles, must be still marked
as basically unsolved. But before we move any further, we will give our
formulation of the problem.

$v$-representability problem: Find an explicit set of densities, preferably as
large as possible, within the $N$-representable densities such that each
element $\rho$ from this set is the density of an (ensemble) ground state of a
self-adjoint Hamiltonian $H=H_{0}+V=-\frac{1}{2}\Delta+W+V$. While $H_{0}$,
the kinetic and interaction contribution, is considered fixed, the external
potential $V$ will depend on the chosen $\rho$.

From the above statement it follows that we also need to investigate for a
chosen set of densities the corresponding set of one-body potentials $v$ and
whether these potentials allow the definition of corresponding self-adjoint
Hamiltonians. Crucial insights in this regard can be obtained from an example
presented in Englisch and Englisch [5], which later has been discussed in more
detail in Chayes  _et al._ [10]. We will here quote some passages from the
latter reference since their discussion already points towards the set of
potentials that we will consider in this work. The example consists of the
density

$\rho(x)=(a+b|x|^{\alpha+1/2})^{2}$ (1)

with $a>b>0$, $1/2>\alpha>0$, that is defined locally around $x=0$, so it
could be on the real line as well as on the one-dimensional torus. Now for a
single particle with a ground-state wave function assumed to be real and
having $\Psi>0$ (this is justified by a result in Lieb and Loss [21, Th.
11.8]), we can take $\rho=\sqrt{\Psi}$. From the single-particle Schrödinger
equation it then follows $v=(\Delta\sqrt{\rho})/(2\sqrt{\rho})$, so we have
immediate access to the representing potential. But it appears that

$\Delta\sqrt{\rho(x)}=b\left(\alpha+\frac{1}{2}\right)|x|^{\alpha-1/2}\delta(x)+b\left(\alpha^{2}-\frac{1}{4}\right)|x|^{\alpha-3/2},$
(2)

so the resulting potential is actually distribution-valued. Chayes  _et al._
[10, p. 514f] write: “This $\rho$ is not the ground state of any Hamiltonian
$-\Delta+v$, having a real-valued potential, although it may be the ground
state of a Hamiltonian with a distributional potential.” And further that the
potential-density mapping “may map very singular potentials into very smooth
densities.” They suggest that the problem might be solved by including
generalized potentials, but that then “much of HK [Hohenberg–Kohn] theory […]
may break down.”

We will take this example and their insights as the starting point here,
indeed showing that by including such generalized potentials the
$v$-representability problem can be solved. But quite contrary to their
statement we would argue that only with this result, DFT has a real
theoretical foundation. The example and our results below demonstrate that
perfectly reasonable densities can only be realized with distributional
potentials. What parts of the theory must then eventually be changed will be
discussed in the conclusions (Section VII). With this, we formulate the main
result of this work.

######  Theorem 1.

Define the following set of $v$-representable densities on the one-dimensional
torus $\mathbb{T}$ and the dual space of (distributional) one-body potentials,
$\displaystyle\mathscr{X}_{>0}$ $\displaystyle=\\{\rho\in
L^{2}(\mathbb{T})\mid\nabla\rho\in L^{2}(\mathbb{T}),\smallint\rho=N,\forall
x\in\mathbb{T}:\rho(x)>0\\}$ (3) $\displaystyle\mathscr{X}^{*}$
$\displaystyle=\\{[v]=\\{v+c\mid c\in\mathbb{R}\\}\mid v=f+\nabla
g\;\text{with}\;f,g\in L^{2}(\mathbb{T})\\}.$ (4) Then for every
$\rho\in\mathscr{X}_{>0}$ there is an equivalence class
$[v]\in\mathscr{X}^{*}$ with corresponding potential
$V=\sum_{j=1}^{N}v(x_{j})$ acting on wave functions, such that $\rho$ is the
density of an (ensemble) ground state of the self-adjoint Hamiltonian
$H=H_{0}+V=-\frac{1}{2}\Delta+W+V$.

The interaction part $W$ is allowed to remain very general, it only needs to
fulfil a property discussed in Corollary 5. Further details about the
definitions in Theorem 1 will be given throughout the following sections.
These discuss the basic setting of the $N$-particle Schrödinger problem on the
one-dimensional torus in terms of quadratic forms (Section II), the considered
spaces for densities and potentials (Section III), the convex formulation of
DFT in the given setting (Section IV), how representing potentials follow from
a statement about subdifferentials (Section V), and the KLMN theorem for the
generalised potentials (Section VI). Finally, the proof of the main theorem
above is given in the concluding Section VII, building on all of the previous
sections. The same section also discusses implications of our findings and
future research directions. Appendix A shows the critical property of ‘kinetic
boundedness’ (Definition 4) for different classes of interactions.

##  II One-dimensional torus setting

Our spatial domain is the one-dimensional torus $\mathbb{T}$, i.e., without
loss of generality we take the interval $[0,1]$ where the points 0 and 1 are
identified and we can recover any other one-dimensional periodic domain by a
simple re-scaling. The reason to call it ‘one-dimensional torus’ and not
simply ‘ring’ or ‘circle’ is that this way we stay closer to the usual
periodic setting in $d>1$ dimensions employed for extended solid-state
systems. On it we consider a fixed number of $N$ spin-$\frac{1}{2}$ particles.
The Hilbert space is then the $N$-fold tensor product of the spaces
$L^{2}(\mathbb{T})\otimes\mathbb{C}^{2}$ where only anti-symmetric wave
functions are considered,

$\mathcal{H}=\bigl{(}L^{2}(\mathbb{T})\otimes\mathbb{C}^{2}\bigr{)}^{\wedge
N}.$ (5)

In what follows, we use the notations $\nabla_{j}$ and $\Delta_{j}$ for the
first and second order (weak) derivatives w.r.t. $x_{j}$, even though the
coordinate is one-dimensional. For any wave function $\Psi\in\mathcal{H}$ we
define the one-particle density on $\mathbb{T}$,

$\rho_{\Psi}(x)=N\sum_{\sigma_{1},\ldots,\sigma_{N}}\int\lvert\Psi(x\sigma_{1},x_{2}\sigma_{2},\dotsc,x_{N}\sigma_{N})\rvert^{2}\mathrm{d}x_{2}\dotsb\mathrm{d}x_{N}.$
(6)

Here, the sum over $\sigma_{1},\dotsc,\sigma_{N}$ sums over all spin
components, while all spatial particle coordinates except for one is
integrated out as well. This density will be our main object of interest. Note
that a $\Psi\in\mathcal{H}$ yields a $\rho_{\Psi}\in L^{1}(\mathbb{T})$ like
this, which does not allow for a point-wise evaluation $\rho_{\Psi}(x)$. This
already points to the beneficial use of Sobolev spaces and indeed a density
will later not only be pointwise well-defined but is even always a continuous
function.

The Hamiltonian of the system is

$H=H_{0}+V=-\frac{1}{2}\Delta+W+V=-\frac{1}{2}\sum_{j=1}^{N}\Delta_{j}+W+V,$
(7)

with the external potential $V=\sum_{j}v(x_{j})$ being defined from a one-body
potential $v$ and the interaction operator $W$ left as a general operator.
This setting will soon be generalized. The expectation value of $V$ with
respect to a wave function $\Psi$ is then given entirely in terms of the
density as

$\langle\Psi,V\Psi\rangle=N\langle\Psi,v\Psi\rangle=\langle
v,\rho_{\Psi}\rangle.$ (8)

Here, the first two angle brackets describe the inner product in $\mathcal{H}$
while the last one is the dual pairing between a density and a potential.
Instead of linear operators we will consider the more general case of
quadratic forms $V(\Phi,\Psi)$ and $W(\Phi,\Psi)$, usually still written as
$\langle\Phi,V\Psi\rangle$ and $\langle\Phi,W\Psi\rangle$.

######  Definition 2 (Reed and Simon [22, Sec. VIII.6]).

A _quadratic form_ is a map $A:Q(A)\times Q(A)\to\mathbb{C}$, where $Q(A)$ is
a dense subspace of $\mathcal{H}$ called the form domain. It is such that
$A(\Phi,\Psi)$ is linear in the second argument and conjugate linear in the
first. If $A(\Phi,\Psi)=\overline{A(\Psi,\Phi)}$ it is called _symmetric_ and
if for all $\Psi\in Q(A)$ it holds $A(\Psi,\Psi)\geq 0$ then it is _positive_.

A self-adjoint or symmetric operator defines a symmetric quadratic form by
$A(\Phi,\Psi)=\langle\Phi,A\Psi\rangle$ and equally a positive operator
defines a positive quadratic form. This makes it possible to take a
distribution-valued $v$ as an external one-body potential (Reed and Simon [22,
Sec. VIII.6, Ex. 1] also give such an example), and analogously to Eq. (8) we
simply have

$V(\Psi,\Psi)=\langle v,\rho_{\Psi}\rangle.$ (9)

Now, the angle bracket means the application of the distribution $v$ on
$\rho_{\Psi}$, which is the same as the dual pairing if a Banach space for
$\rho_{\Psi}$ and its topological dual for $v$ are considered. The full
quadratic form $V(\Phi,\Psi)$ can then be constructed with the polarization
identity. Adding such potentials to the Hamiltonian makes it a quadratic form
as well. This means switching from the usual operator domain of the self-
adjoint operator to the form domain [22, Sec. VIII.6]. In case of the
Laplacian $-\Delta=-\sum_{j}\Delta_{j}$, which is self-adjoint on all
$\Psi\in\mathcal{H}$ that have $\Delta_{j}\Psi$ also $L^{2}$-integrable and
where $\Psi$ and $\nabla_{j}\Psi$ are periodic, the resulting form domain is
relatively easy to see. By symmetry of the wave function and by partial
integration on the torus we have

$\langle\Psi,-\Delta\Psi\rangle=\sum_{j=1}^{N}\langle\Psi,-\Delta_{j}\Psi\rangle=\sum_{j=1}^{N}\langle\nabla_{j}\Psi,\nabla_{j}\Psi\rangle,$
(10)

where the boundary terms cancel due to periodicity. We note that
$\langle\nabla_{j}\Psi,\nabla_{j}\Psi\rangle=\|\nabla_{j}\Psi\rVert_{\mathcal{H}}^{2}$,
so the form domain is constructed from $L^{2}(\mathbb{T})$ functions that have
an $L^{2}$-integrable weak derivative. This gives exactly the Sobolev space
$H^{1}(\mathbb{T})$ with norm
$\|f\rVert^{2}_{H^{1}}=\|f\rVert^{2}_{L^{2}}+\|\nabla f\rVert^{2}_{L^{2}}$. We
find

$Q(-\Delta)=\bigl{(}H^{1}(\mathbb{T})\otimes\mathbb{C}^{2}\bigr{)}^{\wedge
N}.$ (11)

Under certain conditions, adding potentials in the form of quadratic forms to
the Hamiltonian still allows to connect it to a self-adjoint operator. This is
the content of the following classical theorem that is central to our work.

######  Theorem 3 (KLMN theorem, Reed and Simon [23], Th. X.17).

Let $A$ be a positive self-adjoint operator and $B$ a symmetric quadratic form
on $Q(A)\times Q(A)$. Further, let $0\leq a<1$ and $b\geq 0$ such that

$|B(\Psi,\Psi)|\leq a\langle\Psi,A\Psi\rangle+b\langle\Psi,\Psi\rangle$ (12)

for all $\Psi\in Q(A)$. Then there exists a unique self-adjoint operator $C$
with $Q(C)=Q(A)$, bounded below by $-b$, and

$\langle\Phi,C\Psi\rangle=\langle\Phi,A\Psi\rangle+B(\Phi,\Psi)$ (13)

for any $\Phi,\Psi\in Q(C)$.

The theorem also explains why we can write $B(\Phi,\Psi)$ and
$\langle\Phi,B\Psi\rangle$ interchangeably, since we could just set $B=C-A$
with the operators $A,C$ from above. We now get back to our case of
$A=-\Delta$ with the form domain $Q(-\Delta)$ given in Eq. (11). When
normalization is added, this gives us our basic space for wave functions,

$\mathscr{W}=\\{\Psi\in(H^{1}(\mathbb{T})\otimes\mathbb{C}^{2})^{\wedge
N}\mid\|\Psi\|_{\mathcal{H}}=1\\}.$ (14)

We also give the kinetic energy for a wave function $\Psi\in\mathscr{W}$,

$T(\Psi)=-\frac{1}{2}\langle\Psi,\Delta\Psi\rangle=\frac{1}{2}\sum_{j=1}^{N}\|\nabla_{j}\Psi\|_{\mathcal{H}}^{2},$
(15)

which is always finite because of the restriction to the Sobolev space
$H^{1}(\mathbb{T})$ in $\mathscr{W}$. The following property will be demanded
for the potentials that are added to the Hamiltonian in order to still have a
self-adjoint operator by means of the subsequent corollary.

######  Definition 4.

A quadratic form $B$ on $\mathscr{W}\times\mathscr{W}$ is called _kinetically
bounded_ , if there are $a,b\geq 0$ such that for all $\Psi\in\mathscr{W}$ it
holds $|B(\Psi,\Psi)|\leq aT(\Psi)+b$. The infimum of all such $a$ is called
the _relative kinetic bound_.

######  Corollary 5.

Let $W$ be a quadratic form that is positive and kinetically bounded with
relative kinetic bound $<1$. Then the interaction-only Hamiltonian $H_{0}$ is
a self-adjoint operator with form domain $Q(H_{0})=Q(-\Delta)$ that is
positive and for all $\Phi,\Psi\in Q(-\Delta)$ it holds

$\langle\Phi,H_{0}\Psi\rangle=-\frac{1}{2}\langle\Phi,\Delta\Psi\rangle+W(\Phi,\Psi).$
(16)

Take further $V$ as a quadratic form that is kinetically bounded with relative
kinetic bound $<1$. Then the full Hamiltonian $H$ is a self-adjoint operator
with form domain $Q(H)=Q(-\Delta)$ that is bounded below and for all
$\Phi,\Psi\in Q(-\Delta)$ it holds

$\langle\Phi,H\Psi\rangle=-\frac{1}{2}\langle\Phi,\Delta\Psi\rangle+W(\Phi,\Psi)+V(\Phi,\Psi).$
(17)

The proof is just a double application of Theorem 3 (KLMN theorem), first for
$W$, which still leaves the combined operator positive, then for $V$. That
$-\Delta$ is self-adjoint and positive itself was already noted after Eq.
(10). All considered potentials $V$ and $W$ will be such that the above
corollary holds. Lemma 16 shows that all potentials $V$ from the considered
class have the desired property and Appendix A shows it for a very broad
variety of possible interactions $W$.

##  III Densities on the one-dimensional torus

We define the set of _physical_ densities

$\mathscr{I}=\\{\rho\in L^{1}(\mathbb{T})\mid\nabla\sqrt{\rho}\in
L^{2}(\mathbb{T}),\rho\geq 0,\smallint\rho=N\\}.$ (18)

The special significance of this set is revealed by the following two results
that connect it to wave functions and their kinetic-energy content.

######  Lemma 6.

Let $\Psi\in\mathscr{W}$ then $\rho_{\Psi}\in\mathscr{I}$ is implied by the
estimate

$\int\left|\nabla\sqrt{\rho_{\Psi}(x)}\right|^{2}\mathrm{d}x\leq 2T(\Psi).$
(19)

###### Proof.

The properties $\rho_{\Psi}\in L^{1}(\mathbb{T})$, $\rho_{\Psi}\geq 0$ and
$\smallint\rho_{\Psi}=N$ follow straight from the definition of the density in
Eq. (6). The proof for the estimate is then the same as in Lieb [4, Th. 1.1].
With this estimate we also have $\nabla\sqrt{\rho_{\Psi}}\in
L^{2}(\mathbb{T})$ and consequently $\rho_{\Psi}\in\mathscr{I}$. ∎

######  Theorem 7.

There exist two constants $C^{T}_{1},C^{T}_{2}>0$ such that for all
$\rho\in\mathscr{I}$ there is a $\Psi\in\mathscr{W}$ that has
$\rho_{\Psi}=\rho$ ($N$-representability) and the following estimate holds

$T(\Psi)\leq
C^{T}_{1}+C^{T}_{2}\int\left|\nabla\sqrt{\rho(x)}\right|^{2}\mathrm{d}x.$ (20)

###### Proof.

We proceed as in Lieb [4, Th. 1.2]. However, due to boundedness and
periodicity of the domain there are some significant changes. As an ansatz for
$N$-representability, we define $f(x)=(2\pi/N)\int_{0}^{x}\rho(s)\mathrm{d}s$
and for $k=0,\dotsc,N-1$ the orbitals

$\phi_{k}(x)=\left(\frac{\rho(x)}{N}\right)^{1/2}\exp({\mathrm{i}kf(x)}).$
(21)

One can check in Lieb [4, Th. 1.2] that these orbitals are orthonormal and
that the Slater-determinant of all $\phi_{k}$ gives indeed the density $\rho$.
We focus here on estimating the kinetic energy that is the sum of all
$\|\nabla\phi_{k}\|_{L^{2}}^{2}$. We have

$\|\nabla\phi_{k}\rVert_{L^{2}}^{2}=\int\lvert\nabla\phi_{k}(x)\rvert^{2}\mathrm{d}x=\frac{1}{N}\int\left(\left\lvert\nabla\sqrt{\rho(x)}\right\rvert^{2}+k^{2}\rho(x)^{3}\right)\mathrm{d}x.$
(22)

We introduce $g=\sqrt{\rho}$ and define

$A=\int\left|\nabla
g(x)\right|^{2}\mathrm{d}x=\int\left|\nabla\sqrt{\rho(x)}\right|^{2}\mathrm{d}x.$
(23)

Since $\rho\geq 0$ and $\rho$ integrates to $N$, there needs to be one point
$x_{0}\in\mathbb{T}$ where $g(x_{0})\leq\sqrt{N}$ (a pointwise evaluation of
the function is possible because of an embedding into the continuous
functions, see Lemma 8). Then we have

$g(x)^{2}=g(x_{0})^{2}+2\int_{x_{0}}^{x}g(y)\nabla g(y)\mathrm{d}y\leq
N+2\left[\int g(y)^{2}\mathrm{d}y\right]^{1/2}\left[\int|\nabla
g(y)|^{2}\mathrm{d}y\right]^{1/2}=N+2\sqrt{NA},$ (24)

where we used Hölder’s inequality, extended the integration over the full
range of the torus, and that $g(x_{0})\leq\sqrt{N}$ as well as $g^{2}=\rho$.
Using $2\sqrt{NA}\leq N+A$ it then follows $g(x)^{4}\leq 3N^{2}+6NA$. The
second term in (22) can thus be bounded from above by

$\frac{k^{2}}{N}\int\rho(x)^{3}\mathrm{d}x=\frac{k^{2}}{N}\int
g(x)^{2}g(x)^{4}\mathrm{d}x\leq k^{2}(3N^{2}+6NA),$ (25)

and so the statement of Eq. (20) follows from summing $k$ from $0$ to $N-1$. ∎

Although the results above show that $\mathscr{I}$ gives an exhaustive set of
physical densities, i.e., every such density has a representing $N$-particle
wave function and every wave function in $\mathscr{W}$ gives a density in
$\mathscr{I}$, it is still unpractical for our purposes, because by the
condition $\nabla\sqrt{\rho}\in L^{2}(\mathbb{T})$ it does not yield a
_linear_ space. Our main space of interest is instead the _affine space_ of
codimension 1 (having one linear constraint, the normalization of the density)

$\mathscr{X}=\\{\rho\in H^{1}(\mathbb{T})\mid\smallint\rho=N\\}.$ (26)

Remember that a normed space $X\subseteq Y$ is called _continuously embedded_
in $Y$, written $X\hookrightarrow Y$, if there is a constant $C$ such that for
all $\xi\in X$ it holds $\|\xi\|_{Y}\leq C\|\xi\|_{X}$. The following lemma
shows a whole chain of (continuous) embeddings.

######  Lemma 8.

$\mathscr{I}\subset\mathscr{X}\subset H^{1}(\mathbb{T})\hookrightarrow
C^{0}(\mathbb{T})\hookrightarrow L^{\infty}(\mathbb{T})\hookrightarrow
L^{3}(\mathbb{T})\hookrightarrow L^{2}(\mathbb{T})\hookrightarrow
L^{1}(\mathbb{T})$.

###### Proof.

The continuous (even compact) embedding $H^{1}(\mathbb{T})\hookrightarrow
C^{0}(\mathbb{T})$ is a consequence of the Sobolev embedding theorem [24, Th.
8.8]. Here the space of continuous functions is equipped with the usual
maximum norm, so the next continuous embedding
$C^{0}(\mathbb{T})\hookrightarrow L^{\infty}(\mathbb{T})$ follows directly.
The remaining embeddings to the right hold since $\mathbb{T}$ is a finite
domain. Since all $\rho\in\mathscr{I}$ have $\sqrt{\rho}\in L^{2}(\mathbb{T})$
and $\nabla\sqrt{\rho}\in L^{2}(\mathbb{T})$, equivalent to $\sqrt{\rho}\in
H^{1}(\mathbb{T})$, we have for such densities that $\sqrt{\rho}\in
C^{0}(\mathbb{T})$ and consequently also $\rho$ is continuous. This allows the
following estimate,

$\|\nabla\rho\|_{L^{2}}^{2}=\int|\nabla\rho(x)|^{2}\mathrm{d}x=4\int\rho(x)\left|\nabla\sqrt{\rho(x)}\right|^{2}\mathrm{d}x\leq
4\|\rho\|_{L^{\infty}}\|\nabla\sqrt{\rho}\|_{L^{2}}^{2},$ (27)

which shows the first inclusion $\mathscr{I}\subset\mathscr{X}$. Finally,
$\mathscr{X}\subset H^{1}(\mathbb{T})$ trivially holds since $\mathscr{X}$ is
an affine space in $H^{1}(\mathbb{T})$. ∎

For later purposes we will also study the dual space of $\mathscr{X}$ denoted
as $\mathscr{X}^{*}$, later to be identified with the space of external
potentials. The dual space of the Sobolev space $H^{1}(\mathbb{T})$, denoted
$H^{-1}(\mathbb{T})$, is composed of all distributions
$v:H^{1}(\mathbb{T})\to\mathbb{R}$ that can be written as $v=f+\nabla g$ with
$f,g\in L^{2}(\mathbb{T})$. This dual space is endowed with the norm
$\|v\|^{2}_{H^{-1}}=\min\\{\|f\|^{2}_{L^{2}}+\|g\|^{2}_{L^{2}}\mid f,g\in
L^{2}(\mathbb{T}),v=f+\nabla g\\}$ [25, Th. 3.9]. An element $v\in
H^{-1}(\mathbb{T})$ acts on a function $\varphi\in H^{1}(\mathbb{T})$ as
$v(\varphi)=\langle v,\varphi\rangle=\langle f,\varphi\rangle-\langle
g,\nabla\varphi\rangle$. Additionally to $H^{1}(\mathbb{T})$, the affine space
$\mathscr{X}$ is defined by the normalization of its elements to the particle
number $N$, so if two elements of the dual space $\mathscr{X}^{*}$ just differ
by a constant, $v-v^{\prime}=c$, then their difference acts as $\langle
v-v^{\prime},\rho\rangle=cN$ on any arbitrary $\rho\in\mathscr{X}$.
Consequently, $v$ and $v^{\prime}$ cannot be distinguished and the dual space
$\mathscr{X}^{*}$ corresponds to $H^{-1}(\mathbb{T})$ modulo the addition of
constant functions. We thus have $\mathscr{X}^{*}$ given by equivalence
classes,

$\mathscr{X}^{*}=\\{[v]=\\{v+c\mid c\in\mathbb{R}\\}\mid v=f+\nabla
g\;\text{with}\;f,g\in L^{2}(\mathbb{T})\\}.$ (28)

In what follows, we commonly write $v$ instead of $[v]$ and by this mean any
representative of the equivalence class. Physically, if $v$ is a potential,
this form of equivalence is well known, since the addition of a constant only
shifts the energy but does not change any ground-state properties. This
equivalence is due to the _gauge freedom_ that we have when mathematically
representing an electromagnetic field in terms of potentials.

######  Example 1.

The delta distribution $\delta:H^{1}(\mathbb{T})\to\mathbb{R}$ is an element
of $H^{-1}(\mathbb{T})$. For this we have to show that there are $f,g\in
L^{2}(\mathbb{T})$ such that $\delta=f+\nabla g$. Take $g(x)=-x$ and $f(x)=1$
on $[0,1)\simeq\mathbb{T}$. Then
$\delta(\varphi)=\langle\delta,\varphi\rangle=\langle f,\varphi\rangle-\langle
g,\nabla\varphi\rangle$, or written out with integrals,

$\delta(\varphi)=\int_{0}^{1}\varphi(x)\mathrm{d}x+\int_{0}^{1}x\nabla\varphi(x)\mathrm{d}x.$
(29)

Then, by partial integration in the second integral we get

$\delta(\varphi)=\int_{0}^{1}\varphi(x)\mathrm{d}x-\int_{0}^{1}\varphi(x)\mathrm{d}x+x\varphi(x)\Big{|}_{0}^{1}=\varphi(1).$
(30)

By Lemma 8 the $\varphi\in H^{1}(\mathbb{T})$ is continuous and it is also
periodic, thus $\varphi(1)=\varphi(0)$, and we get
$\delta(\varphi)=\varphi(0)$, precisely the action of the delta distribution.

######  Example 2.

How does an element $v=f+\nabla g\in\mathscr{X}^{*}$ act as a potential on
wave functions or directly on densities in the dual pairing $\langle
v,\rho\rangle$ that gives the potential energy? We have $V=\sum_{k}v(x_{k})$
as the associated operator and thus for $V\Psi$ in a distributional sense for
any test function
$\Phi\in(C^{\infty}(\mathbb{T})\otimes\mathbb{C}^{2})^{\otimes N}$ it holds

$\langle\Phi,V\Psi\rangle=\sum_{k=1}^{N}\left(\langle\Phi,f(x_{k})\Psi\rangle-\langle\Phi,g(x_{k})\nabla\Psi\rangle-\langle\nabla\Phi,g(x_{k})\Psi\rangle\right).$
(31)

This has a certain resemblance to the action of the Hamiltonian in a magnetic
Schrödinger equation where terms like $\mathrm{i}\mathbf{A}\cdot\nabla\Psi$
appear. But the complex prefactor $\mathrm{i}$ and the missing
$\langle\nabla\Phi,\cdot\rangle$ term amount to a critical difference. In the
dual pairing the effect of the distributional potential is seen way simpler as

$\langle v,\rho\rangle=\langle f,\rho\rangle-\langle g,\nabla\rho\rangle,$
(32)

thus involving the gradient of the density in a semi-local fashion.

##  IV Convex formulation of density-functional theory

The ground-state problem consists of finding the eigenstate (or the whole
eigenspace in the case of degeneracy) corresponding to the lowest eigenvalue
of $H_{0}+V$. Here, we reformulate this eigenproblem into a more general
variational problem, by seeking the infimum of the spectrum. For this, we rely
on the variational principle for the lower bound of the spectrum of a self-
adjoint operator [26, Th. XIII.1 for $n=1$],

$E(v)=\inf\sigma(H)=\inf\sigma(H_{0}+V)=\inf_{\Psi\in\mathscr{W}}\langle\Psi,(H_{0}+V)\Psi\rangle<\infty.$
(33)

The notation $E(v)$, with the one-body potential $v$ as an argument, is due to
the fact that this $v$ will be the only variable component of the system,
while $H_{0}$ is considered fixed. The definition with an infimum directly
implies for $\lambda\in[0,1]$ that

$\displaystyle E(\lambda v+(1-\lambda)v^{\prime})$
$\displaystyle=\inf_{\Psi\in\mathscr{W}}\langle\Psi,(\lambda
H_{0}+(1-\lambda)H_{0}+\lambda V+(1-\lambda)V^{\prime})\Psi\rangle$ (34)
$\displaystyle\geq\lambda\inf_{\Psi\in\mathscr{W}}\langle\Psi,(H_{0}+V)\Psi\rangle+(1-\lambda)\inf_{\Psi\in\mathscr{W}}\langle\Psi,(H_{0}+V^{\prime})\Psi\rangle=\lambda
E(v)+(1-\lambda)E(v^{\prime}),$

so the functional $E$ is concave. It is possible to separate the variation
over wave functions in Eq. (33) into first a variation over all densities that
follow from wave functions in $\mathscr{W}$, i.e., the set $\mathscr{I}$ by
Lemma 6, and then vary over all wave functions that give such a density,
denoted as $\Psi\mapsto\rho$. We thus have

$E(v)=\inf_{\rho\in\mathscr{I}}\\{\tilde{F}(\rho)+\langle
v,\rho\rangle\\},\quad\quad\tilde{F}(\rho)=\inf_{\Psi\mapsto\rho}\langle\Psi,H_{0}\Psi\rangle.$
(35)

$\tilde{F}$ is called the _pure-state constrained-search functional_. Note
that any $\rho\in\mathscr{X}\setminus\mathscr{I}$ is either not a (non-
negative) density or it is connected to infinite kinetic energy by the
estimate Eq. (19) in Lemma 6. In both cases we set $\tilde{F}(\rho)=\infty$.
Consequently, the variation in Eq. (35) can be extended from $\mathscr{I}$ to
$\mathscr{X}$ without changing the result and we get

$E(v)=\inf_{\rho\in\mathscr{X}}\\{\tilde{F}(\rho)+\langle v,\rho\rangle\\}.$
(36)

Next, define a functional $F$ on $\mathscr{X}$ by means of convex conjugation
(Legendre–Fenchel transformation),

$F(\rho)=\sup_{v\in\mathscr{X}^{*}}\\{E(v)-\langle v,\rho\rangle\\}.$ (37)

This gives $F$ as the biconjugate of $\tilde{F}$ and as such $F$ is convex and
lower-semicontinuous (even weakly lower-semicontinuous) and $F\leq\tilde{F}$
[27, Prop. 2.19]. We can also arrive at $F$ by taking the lower-semicontinuous
convex envelope of $\tilde{F}$ [27, Cor. 2.23] or by generalizing the
constrained search in Eq. (35) from wave functions to density matrices. This
fact will later prove critical for showing the final step for
$v$-representability, in order to connect the density-functional formulation
back to quantum-mechanical states.

######  Theorem 9.

The functional $F$ on $\mathscr{I}$ is equal to the constrained search over
density matrices, i.e., for all $\rho\in\mathscr{I}$ it holds

$F(\rho)=F_{\mathrm{DM}}(\rho)=\inf_{\Gamma\mapsto\rho}\\{\operatorname{Tr}\Gamma
H_{0}\\}.$ (38)

Here, all considered density matrices are (possibly infinite) convex sums
$\Gamma=\sum_{k}\lambda_{k}\langle\Psi_{k},\cdot\rangle\Psi_{k}$,
$\lambda_{k}\geq 0$, $\sum_{k}\lambda_{k}=1$, with $\Psi_{k}\in\mathscr{W}$
all orthogonal. The density is then also a convex sum,
$\rho=\sum_{k}\lambda_{k}\rho_{k}$, with $\Psi_{k}\mapsto\rho_{k}$.

###### Proof.

For the first part, $F=F_{\mathrm{DM}}$, remember that Lieb [4, Th. 4.3] shows
that $F_{\mathrm{DM}}=F_{\text{Lieb}}$, with the latter also equal to the
lower-semicontinuous convex envelope of the same $\tilde{F}$ but w.r.t. the
$L^{1}\cap L^{3}$-topology. Now with the embeddings from Lemma 8 we have that
the topology of $H^{1}(\mathbb{T})$ is finer than those of $L^{1}(\mathbb{T})$
and $L^{3}(\mathbb{T})$ and thus also of $L^{1}(\mathbb{T})\cap
L^{3}(\mathbb{T})$. This implies that a function that is continuous (or lower-
semicontinuous) w.r.t. the $L^{1}\cap L^{3}$-norm is also continuous (or
lower-semicontinuous) w.r.t. the $H^{1}$-norm. In other words, the set of all
functions continuous (or lower-semicontinuous) w.r.t. the $H^{1}$-norm is
larger than the set of functions continuous (or lower-semicontinuous) w.r.t.
the coarser topology. We conclude that the convex envelope (as the pointwise
supremum over all convex, lower-semicontinuous functions that lie below)
w.r.t. the $H^{1}$-topology is always larger or equal to the convex envelope
w.r.t. the $L^{1}\cap L^{3}$-topology, thus $F\geq
F_{\text{Lieb}}=F_{\mathrm{DM}}$. The inequality in the other direction
follows easily by noting that the ground-state energy can be equally well
defined by considering density matrices instead of pure states,

$E(v)=\inf_{\rho\in\mathscr{I}}\\{F_{\mathrm{DM}}(\rho)+\langle
v,\rho\rangle\\},$ (39)

and so for all $\rho\in\mathscr{I}$ we have $E(v)\leq
F_{\mathrm{DM}}(\rho)+\langle v,\rho\rangle$. But by definition

$F(\rho)=\sup_{v\in\mathscr{X}^{*}}\\{E(v)-\langle
v,\rho\rangle\\}\leq\sup_{v\in\mathscr{X}^{*}}F_{\mathrm{DM}}(\rho)=F_{\mathrm{DM}}(\rho),$
(40)

which shows $F=F_{\mathrm{DM}}=F_{\text{Lieb}}$ on $\mathscr{I}$.
That the density matrix is always a sum of orthogonal and normalized states
$\Psi_{k}\in\mathcal{H}$ is a standard result [4, Sec. 4.B]. That they are all
$\Psi_{k}\in\mathscr{W}$ is clear since the kinetic energy expectation value
would be infinite otherwise. ∎

A less technical view on the different density functionals can be found in a
recent review together with a detailed discussion on the relevance of the
Hohenberg–Kohn theorem [28]. A recent comprehensive, mathematical study is
also available [8].

Note that something critical happened in Eq. (37) above. The potential $v$ is
now from the space $\mathscr{X}^{*}$ that could in general be different from
the class of allowed potentials as demanded by Corollary 5. This means one
potentially loses the connection to a self-adjoint operator formalism and
switches entirely to a convex analysis treatment of the ground-state problem.
In Section VI we will take special care of this and make the connection back
to a self-adjoint Hamiltonian with potential $v\in\mathscr{X}^{*}$.

Finally, by biconjugation it is also possible to go back to $E$ from $F$ with
the opposite convex conjugation,

$E(v)=\inf_{\rho\in\mathscr{X}}\\{F(\rho)+\langle v,\rho\rangle\\}.$ (41)

This shows the common feature of all formulations of DFT: The previous ground-
state problem given by a _linear_ many-particle Schrödinger equation with its
high-dimensional configuration space is translated into a variational problem
incorporating a _non-linear_ convex functional on just the one-particle
density space. But note that there is a critical difference between this way
of defining the so-called _universal density functional_ $F$ and how this is
done in most of the DFT literature. This difference lies in our restriction to
the affine space $\mathscr{X}$ where the normalization to the particle number
is already included. If $F$ is defined on a larger space like
$H^{1}(\mathbb{T})$ then arbitrarily close to any physical density there lie
elements that are not even $N$-representable. This can be avoided if we limit
ourselves to the affine space $\mathscr{X}$ as we will see in the next
section. One important ingredient is the following bound on $F(\rho)$ for
$\rho\in\mathscr{I}$.

######  Lemma 10.

There exist two constants $C_{1}^{F},C_{2}^{F}>0$ such that for all
$\rho\in\mathscr{I}$ it holds

$\frac{1}{2}\|\nabla\sqrt{\rho}\|_{L^{2}}^{2}\leq F(\rho)\leq
C_{1}^{F}+C_{2}^{F}\|\nabla\sqrt{\rho}\|_{L^{2}}^{2}.$ (42)

If $\rho\notin\mathscr{I}$ then $F(\rho)=\infty$.

###### Proof.

Define $G(\rho)=\frac{1}{2}\|\nabla\sqrt{\rho}\|_{L^{2}}^{2}$ for
$\rho\in\mathscr{I}$ and $G(\rho)=\infty$ otherwise. Then for the lower bound
use that for all $\rho\in L^{1}\cap L^{3}$ Lieb [4, Th. 3.8] already showed
that $G(\rho)\leq F(\rho)$. This already implies that for
$\rho\notin\mathscr{I}$ we have $G(\rho)=F(\rho)=\infty$.
For the upper bound we employ Eq. (20) of Theorem 7 that states that for
$\rho\in\mathscr{I}$ there is a $\Psi\in\mathscr{W}$ with $\rho_{\Psi}=\rho$
and $T(\Psi)\leq C_{1}^{T}+C_{2}^{T}\|\nabla\sqrt{\rho}\|_{L^{2}}^{2}$. For
the interaction term we have due to the kinetic boundedness according to
Definition 4 $|\langle\Psi,W\Psi\rangle|\leq aT(\Psi)+b$. Both estimates
together give an upper bound for $F(\rho)$,

$F(\rho)\leq\tilde{F}(\rho)\leq\langle\Psi,H_{0}\Psi\rangle\leq
T(\Psi)+\langle\Psi,W\Psi\rangle\leq
C_{1}^{F}+C_{2}^{F}\|\nabla\sqrt{\rho}\|_{L^{2}}^{2},$ (43)

which concludes the proof. ∎

##  V $v$-Representing potentials from non-empty subdifferentials

In the previous section we found the ground-state energy and density from the
variational problem

$E(v)=\inf_{\rho\in\mathscr{X}}\\{F(\rho)+\langle v,\rho\rangle\\},$ (44)

where $F$ is a convex and lower-semicontinuous functional on $\mathscr{X}$. A
convex functional generally allows the definition of the _subdifferential_ at
$\rho\in\mathscr{X}$ as a set in the dual space $\mathscr{X}^{*}$,

$\partial
F(\rho)=\\{u\in\mathscr{X}^{*}\mid\forall\rho^{\prime}\in\mathscr{X}:F(\rho)-F(\rho^{\prime})\leq\langle
u,\rho-\rho^{\prime}\rangle\\}.$ (45)

An element $u\in\partial F(\rho)$ is then called a _subgradient_. In general,
the subdifferential may be empty. If $F$ is (Gateaux) differentiable at some
point $\rho\in\mathscr{X}$ then $\partial F(\rho)$ contains just a single
element that gives exactly the derivative. Note that while the gradient is a
local property for general functionals, the subgradient contains _global_
information of the convex functional through the inequality $F(\rho)-\langle
u,\rho\rangle\leq F(\rho^{\prime})-\langle u,\rho^{\prime}\rangle$ which holds
for any $\rho^{\prime}\in\mathscr{X}$ if $u\in\partial F(\rho)$. If we just
replace $u=-v$ then we see that this is exactly the condition for minimizing
the variational problem of Eq. (44). This means the variational problem of
finding the density that minimizes the functional $\rho\mapsto F(\rho)+\langle
v,\rho\rangle$ can equivalently be replaced by the condition $-v\in\partial
F(\rho)$. This in turn means that if $\partial F(\rho)$ is non-empty, we know
that at least one $v\in\mathscr{X}^{*}$ exists for which $\rho$ is a minimizer
in Eq. (44).

The critical ingredients in order to get a non-empty subdifferential are then
the following definition and two lemmata.

######  Definition 11.

The _effective domain_ of a convex functional $F$ on $\mathscr{X}$ consists of
all points $\rho\in\mathscr{X}$ where $F(\rho)<\infty$.

######  Lemma 12 (Barbu and Precupanu [27], Th. 2.14).

Let $F$ be a convex functional on $\mathscr{X}$ then $F$ is continuous on the
whole interior of its effective domain if and only if $F$ is uniformly bounded
from above on a neighborhood of any interior point of its effective domain.

######  Lemma 13 (Barbu and Precupanu [27], Prop. 2.36).

If a convex functional $F$ on $\mathscr{X}$ is continuous at some
$\rho\in\mathscr{X}$ then it has a non-empty subdifferential at this point.

The strategy is then clear. We need to identify a set in the interior of the
effective domain of $F$ and show the necessary bound from above. For this we
define a set of densities in $\mathscr{X}$ that are gapped away from zero,

$\mathscr{X}_{>0}=\\{\rho\in\mathscr{X}\mid\forall
x\in\mathbb{T}:\rho(x)>0\\}.$ (46)

######  Lemma 14.

The set $\mathscr{X}_{>0}$ is open in $\mathscr{X}$.

###### Proof.

Note that the set $\\{f\in L^{\infty}(\mathbb{T})\mid\forall
x\in\mathbb{T}:f(x)>0\\}$ is open in $L^{\infty}(\mathbb{T})$. Because of the
continuous embedding $H^{1}(\mathbb{T})\hookrightarrow L^{\infty}(\mathbb{T})$
from Lemma 8 the set $\\{f\in H^{1}(\mathbb{T})\mid\forall
x\in\mathbb{T}:f(x)>0\\}$ is then open in $H^{1}(\mathbb{T})$. Hence the
restriction to the affine subspace $\mathscr{X}$ is also open. ∎

We then have the desired result that guarantees us a non-empty
subdifferential.

######  Theorem 15.

For any $\rho\in\mathscr{X}_{>0}$ it holds $\rho\in\mathscr{I}$ and the
subdifferential $\partial F(\rho)\subset\mathscr{X}^{*}$ is non-empty.

###### Proof.

Any $\rho\in\mathscr{X}_{>0}$ is continuous on $\mathbb{T}$ by Lemma 8, so we
can choose a minimum $\eta=\min_{x}\\{\rho(x)\\}>0$. We then have the estimate

$\|\nabla\sqrt{\rho}\|_{L^{2}}^{2}=\int\left|\nabla\sqrt{\rho(x)}\right|^{2}\mathrm{d}x=\frac{1}{4}\int\frac{|\nabla\rho(x)|^{2}}{\rho(x)}\mathrm{d}x<\frac{1}{4\eta}\int|\nabla\rho(x)|^{2}\mathrm{d}x=\frac{1}{4\eta}\|\nabla\rho\|_{L^{2}}^{2}\leq\frac{1}{4\eta}\|\rho\|_{H^{1}}^{2}.$
(47)

This inequality already establishes that $\rho\in\mathscr{I}$. Together with
Lemma 10 this also means that $F$ is bounded above at any such
$\rho\in\mathscr{X}_{>0}$. Thus every $\rho\in\mathscr{X}_{>0}$ belongs to the
effective domain of $F$ and since by Lemma 14 $\mathscr{X}_{>0}$ is open, the
whole set must belong to the interior of the effective domain of $F$. By the
continuous embedding $H^{1}(\mathbb{T})\hookrightarrow L^{\infty}(\mathbb{T})$
from Lemma 8, we can choose a neighborhood for any $\rho\in\mathscr{X}_{>0}$
that is contained within an $\varepsilon$-ball around $\rho$ in
$L^{\infty}$-norm. Let $\varepsilon<\eta$, then an estimate like in Eq. (47)
holds uniformly on the whole neighborhood of $\rho\in\mathscr{X}_{>0}$. This
carries over to a uniform bound for $F(\rho)$. Thus, Lemma 12 becomes
applicable and $F$ is continuous on $\mathscr{X}_{>0}$. Lemma 13 then
concludes the proof. ∎

This proves that every $\rho\in\mathscr{X}_{>0}$ is ‘representable’ by a
$v\in-\partial F(\rho)\subset\mathscr{X}^{*}$ in the sense that
$E(v)=F(\rho)+\langle v,\rho\rangle$, i.e., $\rho$ is a minimizer in the
variational problem of Eq. (44). Yet, we do not consider this already as
sufficient for showing real $v$-representability, since it should also be
guaranteed that this $v$ can be the external potential of a valid Hamiltonian.

##  VI KLMN theorem for distributional potentials

In order to show that $H=H_{0}+V$ with $V$ defined by a $v\in\mathscr{X}^{*}$
as described in Section II is connected to a valid self-adjoint operator we
have to look to Corollary 5. The critical condition on $V$ for this corollary
to hold is that it is kinetically bounded with relative kinetic bound $<1$.
Here we will show that the relative kinetic bound is $0$ (which does not mean
that the $T(\Psi)$ in the estimate vanishes, just that the prefactor can be
made arbitrarily small).

######  Lemma 16.

Let $v\in\mathscr{X}^{*}$ then this defines a kinetically bounded form $V$ on
$\mathscr{W}\times\mathscr{W}$ with relative kinetic bound $0$.

###### Proof.

This proof is inspired by the treatment of distributional potentials in the
context of the KLMN theorem by Herczyński [29]. According to Eq. (8) and
Definition 4 we have to test $|\langle v,\rho_{\Psi}\rangle|\leq\varepsilon
T(\Psi)+b_{\varepsilon}$ for all $\Psi\in\mathscr{W}$ where $\varepsilon$ can
be taken arbitrarily small. Taking $\varepsilon$ smaller will typically mean
that the other bound $b_{\varepsilon}$ grows. According to Eq. (28) any
$v\in\mathscr{X}^{*}$ can we written as $v=f+\nabla g$ where $f,g\in
L^{2}(\mathbb{T})$. We use the triangle inequality,

$|\langle v,\rho_{\Psi}\rangle|\leq|\langle
f,\rho_{\Psi}\rangle|+|\langle\nabla g,\rho_{\Psi}\rangle|,$ (48)

and we will continue with each term individually. We choose a sequence
$\\{g_{n}\\}_{n}$ in $C^{\infty}(\mathbb{T})$ such that
$\|g-g_{n}\|_{L^{2}}\to 0$. We take the supremum of $\nabla g_{n}$ out of the
inner product and get by the definition of the distributional part of the
potential that

$|\langle\nabla g,\rho_{\Psi}\rangle|\leq|\langle\nabla
g_{n},\rho_{\Psi}\rangle|+|\langle\nabla(g-g_{n}),\rho_{\Psi}\rangle|\leq
N\|\nabla g_{n}\|_{L^{\infty}}+|\langle g-g_{n},\nabla\rho_{\Psi}\rangle|\leq
N\|\nabla
g_{n}\|_{L^{\infty}}+\|g-g_{n}\|_{L^{2}}\|\nabla\rho_{\Psi}\|_{L^{2}}.$ (49)

From Eq. (27) we have $\|\nabla\rho_{\Psi}\|_{L^{2}}\leq
2\|\sqrt{\rho_{\Psi}}\|_{L^{\infty}}\|\nabla\sqrt{\rho_{\Psi}}\|_{L^{2}}$ and
from the continuous embedding $H^{1}(\mathbb{T})\hookrightarrow
L^{\infty}(\mathbb{T})$ of Lemma 8 that
$\|\sqrt{\rho_{\Psi}}\|_{L^{\infty}}\leq C\|\sqrt{\rho_{\Psi}}\|_{H^{1}}$.
Since $\|\nabla\sqrt{\rho_{\Psi}}\|_{L^{2}}\leq\|\sqrt{\rho_{\Psi}}\|_{H^{1}}$
this combines to $\|\nabla\rho_{\Psi}\|_{L^{2}}\leq
2C\|\sqrt{\rho_{\Psi}}\|_{H^{1}}^{2}=2C(N+\|\nabla\sqrt{\rho_{\Psi}}\|_{L^{2}}^{2})\leq
2C(N+2T(\Psi))$ by Eq. (19) in Lemma 6. So the estimate is now

$|\langle\nabla g,\rho_{\Psi}\rangle|\leq N\|\nabla
g_{n}\|_{L^{\infty}}+2C\|g-g_{n}\|_{L^{2}}(N+2T(\Psi)).$ (50)

For $|\langle f,\rho_{\Psi}\rangle|$ the estimate is simpler since no partial
integration is necessary. Again, we introduce a sequence $\\{f_{n}\\}_{n}$ in
$C^{\infty}(\mathbb{T})$ such that $\|f-f_{n}\|_{L^{2}}\to 0$ and have like in
Eq. (49) that

$|\langle f,\rho_{\Psi}\rangle|\leq
N\|f_{n}\|_{L^{\infty}}+\|f-f_{n}\|_{L^{2}}\|\rho_{\Psi}\|_{L^{2}}.$ (51)

We then use
$\|\rho_{\Psi}\|_{L^{2}}\leq\|\sqrt{\rho_{\Psi}}\|_{L^{\infty}}\|\sqrt{\rho_{\Psi}}\|_{L^{2}}$
and from Lemma 8 also the continuous embedding
$H^{1}(\mathbb{T})\hookrightarrow L^{2}(\mathbb{T})$ to get
$\|\rho_{\Psi}\|_{L^{2}}\leq C^{\prime}\|\sqrt{\rho_{\Psi}}\|_{H^{1}}^{2}\leq
C^{\prime}(N+2T(\Psi))$. Then the second estimate is

$|\langle f,\rho_{\Psi}\rangle|\leq
N\|f_{n}\|_{L^{\infty}}+C^{\prime}\|f-f_{n}\|_{L^{2}}(N+2T(\Psi))$ (52)

and we can combine this to

$|\langle v,\rho_{\Psi}\rangle|\leq N\left(\|f_{n}\|_{L^{\infty}}+\|\nabla
g_{n}\|_{L^{\infty}}\right)+\left(2C\|g-g_{n}\|_{L^{2}}+C^{\prime}\|f-f_{n}\|_{L^{2}}\right)(N+2T(\Psi)).$
(53)

But by increasing the index $n$ the $\|f-f_{n}\|_{L^{2}}$ and
$\|g-g_{n}\|_{L^{2}}$ can be taken arbitrarily small such that the whole
prefactor for $T(\Psi)$ drops below any $\varepsilon>0$, while the remaining
terms increase to some large but finite value. ∎

This result makes Theorem 3 (KLMN theorem) applicable to potentials
$v\in\mathscr{X}^{*}$. We will go one step further with the next theorem,
showing also the corresponding Hamiltonian always allows for a ground state in
the appropriate sense. As a preparatory step we show the following equivalence
of norms.

######  Lemma 17.

Let $v\in\mathscr{X}^{*}$ then the corresponding potential form $V$ defines a
self-adjoint Hamiltonian $H=H_{0}+V$ that is bounded below. With an
appropriate shift of the potential $v$ by adding a constant, the square root
of $\langle\Psi,H\Psi\rangle$ is equivalent to the $H^{1}$-norm of any
$\Psi\in\mathscr{W}$.

###### Proof.

That any $v\in\mathscr{X}^{*}$ defines a self-adjoint Hamiltonian that is
bounded below is proven if Lemma 16 and Corollary 5 are combined. By Lemma 16
for $V$ and by assumption for $W$, both quadratic forms have relative kinetic
bounds $<1$, so we can find $0\leq a_{V},a_{W}<1$ and $b_{V},b_{W}\geq 0$ like
in Definition 4 where $a_{V}$ can even be taken arbitrarily small, so
$a_{V}+a_{W}<1$. Now reordering the expectation value
$\langle\Psi,H\Psi\rangle=T(\Psi)+\langle\Psi,W\Psi\rangle+\langle\Psi,V\Psi\rangle$
gives the estimate

$T(\Psi)\leq\langle\Psi,H\Psi\rangle+|\langle\Psi,W\Psi\rangle|+|\langle\Psi,V\Psi\rangle|\leq\langle\Psi,H\Psi\rangle+aT(\Psi)+b,$
(54)

where we introduced $0\leq a=a_{V}+a_{W}<1$ and $b=b_{V}+b_{W}\geq 0$. Thus we
have

$(1-a)T(\Psi)-b\leq\langle\Psi,H\Psi\rangle$ (55)

as a lower bound for $\langle\Psi,H\Psi\rangle$. Next, remember that a
potential $v\in\mathscr{X}^{*}$ is only defined up to an additive constant as
an equivalence class, so we are actually free to choose an arbitrary additive
constant for $H$. Let $H^{\prime}=H+(1-a)/2+b$ then

$\frac{1-a}{2}(2T(\Psi)+1)\leq\langle\Psi,H^{\prime}\Psi\rangle.$ (56)

Since $2T(\Psi)+1$ is just the square of the $H^{1}$-norm for a (normalized)
$\Psi\in\mathscr{W}$, this gives the lower bound for the equivalence. The
upper bound is seen more directly from the same kinetic bounds for the
expectation values of $V$ and $W$.

$\langle\Psi,H^{\prime}\Psi\rangle=\langle\Psi,H\Psi\rangle+\frac{1-a}{2}+b\leq(1+a)T(\Psi)+\frac{1-a}{2}+2b\leq\max\left\\{\frac{1+a}{2},\frac{1-a}{2}+2b\right\\}(2T(\Psi)+1).$
(57)

This establishes the equivalence between the square root of
$\langle\Psi,H^{\prime}\Psi\rangle$ and the $H^{1}$-norm of $\Psi$. ∎

######  Theorem 18.

Let $v\in\mathscr{X}^{*}$ then the corresponding potential form $V$ defines a
self-adjoint Hamiltonian $H=H_{0}+V$ that is bounded below. This $H$ has a
minimizer $\Psi_{0}\in\mathscr{W}$ in Eq. (33),
$E(v)=E_{0}=\langle\Psi_{0},(H_{0}+V)\Psi_{0}\rangle$, that is also a solution
to the Schrödinger equation $H\Psi_{0}=E_{0}\Psi_{0}$ (in a distributional
sense).

###### Proof.

Like in Lemma 17 we know that any $v\in\mathscr{X}^{*}$ defines a self-adjoint
Hamiltonian that is bounded below. Next, take a sequence
$\\{\Psi_{n}\\}_{n}\in\mathscr{W}$ that realizes the infimum in Eq. (33), so
$\langle\Psi_{n},H\Psi_{n}\rangle\to E_{0}$. Consequently, for any
$\varepsilon>0$ there is a $n_{0}\in\mathbb{N}$ such that for all $n>n_{0}$ it
holds $\langle\Psi_{n},H\Psi_{n}\rangle<E_{0}+\varepsilon$. By Lemma 16 for
$V$ and by assumption for $W$, both quadratic forms have relative kinetic
bounds $<1$, so we can find $a_{V},a_{W}<1$ like in Definition 4 where $a_{V}$
can even be taken arbitrarily small, so $a_{V}+a_{W}<1$. With the estimate
from the kinetic bound it holds for all $n>n_{0}$ that

$E_{0}+\varepsilon>\langle\Psi_{n},H\Psi_{n}\rangle=T(\Psi_{n})+\langle\Psi_{n},W\Psi_{n}\rangle+\langle\Psi_{n},V\Psi_{n}\rangle\geq(1-a_{W}-a_{V})T(\Psi_{n})-b_{W}-b_{V}.$
(58)

But for all (normalized) $\Psi_{n}\in\mathscr{W}$ the $T(\Psi_{n})$ is
equivalent to the $H^{1}$-norm squared, so the above inequality means that all
$\Psi_{n}$ with $n>n_{0}$ are in a $H^{1}$-bounded set. Since the Sobolev
space $H^{1}$ is reflexive, this guarantees that $\\{\Psi_{n}\\}_{n}$ has a
subsequence that converges weakly in $H^{1}$ [30, Cor. 11.9]. Further, since
the torus is a bounded domain, by the Rellich–Kondrachov theorem [25, Th. 6.3]
$H^{1}$ is compactly embedded in $L^{2}$. This means the previously chosen
subsequence permits another subsequence, again called $\\{\Psi_{n}\\}_{n}$,
that converges strongly in $L^{2}$ to a limit that we call $\Psi_{0}$. Because
of the weak convergence in $H^{1}$ and the strong convergence in $L^{2}$ we
know that $\Psi_{0}\in\mathscr{W}$ which implies $\|\Psi_{0}\|_{L^{2}}=1$.
From Lemma 17 we now know that with an appropriate constant shift we can
define $H^{\prime}=H+c$, where the square root of the expectation value
$\langle\Psi,H^{\prime}\Psi\rangle$ defines a norm equivalent to the
$H^{1}$-norm. Since the norm is non-increasing under weak limits [30, Thm.
$11.2$] and any expectation value of $H$ with respect to
$\Psi_{0}\in\mathscr{W}$ must be bounded below by $E_{0}$, we have

$E_{0}+c\leq\langle\Psi_{0},(H+c)\Psi_{0}\rangle=\langle\Psi_{0},H^{\prime}\Psi_{0}\rangle\leq\lim_{n\to\infty}\langle\Psi_{n},H^{\prime}\Psi_{n}\rangle=E_{0}+c$
(59)

and thus $E_{0}=\langle\Psi_{0},H\Psi_{0}\rangle$.
For showing that this $\Psi_{0}$ solves the Schrödinger equation in a
distributional sense we follow Lieb and Loss [21, Th. 11.8]. Take any test
function $\Phi$ and note with $\alpha\in\mathbb{R}$ from the variational
principle of Eq. (33) that

$E_{0}\|\Psi_{0}+\alpha\Phi\|^{2}_{\mathcal{H}}\leq\langle(\Psi_{0}+\alpha\Phi),H(\Psi_{0}+\alpha\Phi)\rangle,$
(60)

where the norm of the vector $\Psi_{0}+\alpha\Phi$ on the left appears due to
normalization. Writing the squares all out and using that
$E_{0}=\langle\Psi_{0},H\Psi_{0}\rangle$ we have

$0\leq
2\alpha\mathop{\mathsf{missing}}{Re}\langle\Phi,(H-E_{0})\Psi_{0}\rangle+\alpha^{2}(\langle\Phi,H\Phi\rangle-
E_{0}\|\Phi\|^{2}_{\mathcal{H}}).$ (61)

Now this equation needs to hold for $\alpha$ arbitrarily small and so it
follows

$0\leq
2\alpha\mathop{\mathsf{missing}}{Re}\langle\Phi,(H-E_{0})\Psi_{0}\rangle$ (62)

and since $\alpha$ can have any sign we conclude that
$\mathop{\mathsf{missing}}{Re}\langle\Phi,(H-E_{0})\Psi_{0}\rangle=0$.
Replacing $\Phi$ by $\mathrm{i}\Phi$ gives
$\mathop{\mathsf{missing}}{Im}\langle\Phi,(H-E_{0})\Psi_{0}\rangle=0$ and
consequently $\langle\Phi,(H-E_{0})\Psi_{0}\rangle=0$ for arbitrary test
functions $\Phi$ which concludes the proof. ∎

##  VII Conclusions and outlook

To prove our main result from the introduction is now just a matter of
collecting the results from the previous sections.

###### Proof of Theorem 1.

Let $\rho\in\mathscr{X}_{>0}$, then Theorem 15 shows that $\rho\in\mathscr{I}$
and $\partial F(\rho)$ is non-empty. Take $-v\in\partial
F(\rho)\subset\mathscr{X}^{*}$ then $\rho$ is a minimizer in the variational
problem Eq. (44) and $E(v)=F(\rho)+\langle v,\rho\rangle$ (see Section V). In
order to connect to quantum-mechanical ground states we need the constrained-
search functional. We cannot use $\tilde{F}$, since $F\neq\tilde{F}$ in
general, but we have $F=F_{\mathrm{DM}}$ on $\mathscr{I}$ by Theorem 9. By
this result,

$E_{0}=E(v)=F_{\mathrm{DM}}(\rho)+\langle
v,\rho\rangle=\inf_{\Gamma\mapsto\rho}\\{\operatorname{Tr}\Gamma
H_{0}\\}+\langle
v,\rho\rangle=\inf_{\Gamma\mapsto\rho}\\{\operatorname{Tr}\Gamma H\\},$ (63)

where we simply recombined $H_{0}+V=H$. Now as in Theorem 18 we next want to
show that this infimum allows for a minimizer $\Gamma_{0}\mapsto\rho$. This is
achieved by Lieb [4, Cor. 4.5(ii)] that can be applied in the current setting
since $F_{\mathrm{DM}}$ is defined in exactly the same manner (with the only
difference of a one-dimensional torus domain instead of $\mathbb{R}^{3}$ which
only simplifies things) and every $\rho\in\mathscr{X}_{>0}$ is also in
$\mathscr{I}$ by Theorem 15, the set considered in the reference. Now, by
Theorem 9 such a $\Gamma_{0}\mapsto\rho$ can be constructed as a (possibly
infinite) convex sum
$\Gamma_{0}=\sum\lambda_{k}\langle\Psi_{0,k},\cdot\rangle\Psi_{0,k}$ with
$\Psi_{0,k}\in\mathscr{W}$ all orthogonal. Assume that there is one
$\Psi_{0,l}$ in this set with $\lambda_{l}>0$ that has an energy expectation
value that is above the ground-state energy, i.e.,
$\langle\Psi_{0,l},H\Psi_{0,l}\rangle=E_{0}+\varepsilon$ with $\varepsilon>0$.
Then it follows immediately from Eq. (63) that
$E_{0}=\sum_{k}\lambda_{k}E_{0}+\lambda_{l}\varepsilon=E_{0}+\lambda_{l}\varepsilon$,
a clear contradiction. Consequently, all $\Psi_{0,k}$ that appear in
$\Gamma_{0}$ with $\lambda_{k}>0$ have
$\langle\Psi_{0,k},H\Psi_{0,k}\rangle=E_{0}$ and thus are ground-state
solutions for the Schrödinger equation $H\Psi_{0,k}=E_{0}\Psi_{0,k}$ in a
distributional sense as demonstrated in Theorem 18. Since it also holds from
$\Gamma_{0}\mapsto\rho$ that $\sum\lambda_{k}\rho_{k}=\rho$ with
$\Psi_{0,k}\mapsto\rho_{k}$, $\rho$ is indeed $v$-representable. In case we
need more than a single $\Psi_{0,k}$ to arrive at $\rho$, $v$ yields a
degenerate ground state and $\rho$ is the density of an ensemble state. ∎

That this result only holds for one-dimensional periodic domains can of course
be seen as a severe limitation, yet the considered torus geometries are
customarily employed in condensed matter physics in the form of Born–von-
Kármán periodic boundary conditions. But also one-dimensional models are
widespread in this field, e.g., the Kronig–Penney model [31] that even
frequently employs a periodic _distributional_ potential in the form of delta
peaks. This means the presented theory of $v$-representability, apart from
solving this important problem of DFT, is also immediately relevant for models
of condensed matter physics. With our main result proven, we want to conclude
with a list of follow-up questions and ideas how to generalize our findings.

  * •

A first open question would be that of the ‘maximality’ of the proposed
$v$-representable density set $\mathscr{X}_{>0}$. So, is any density with
$\rho(x)=0$ either not $v$-representable at all or can it only be represented
by a $v\notin\mathscr{X}^{*}$? This is equivalent to the statement that any
$v\in\mathscr{X}^{*}$ yields a ground-state density that is gapped away from
zero, i.e., $\rho\in\mathscr{X}_{>0}$. Such a result would close the circle
and allow for a well-defined mapping between the potential and density sets
$\mathscr{X}^{*}$ and $\mathscr{X}_{>0}$ (yet not necessarily a one-to-one
mapping due to potential non-unique $v$-representability and the possible
occurrence of degeneracies).

  * •

Another obvious open question is the possibility of a generalization to a
$(d\geq 2)$-dimensional torus domain or other more-dimensional compact
manifolds. The current method critically uses the embedding
$H^{1}=W^{1,2}\hookrightarrow L^{\infty}$ from Lemma 8 which only holds for
$d=1$. This is then employed for showing that $\mathscr{X}_{>0}$ is open
(Lemma 14) and also to have a kinetic bound for potentials from
$\mathscr{X}^{*}$ (Lemma 16). Other Sobolev embeddings are available, and in
$d\geq 2$ dimensions $W^{1,p}\hookrightarrow L^{\infty}$ for $p>d$ could be
applied. But Lemma 16 also critically relies on an estimate of the
$W^{1,2}$-norm by the kinetic energy, which is again needed to show
applicability of the KLMN theorem, and the same is not achievable for
$W^{1,p}$ with $p>2$. Consequently, a generalization to higher dimensions
requires a different approach. In any case, we would not expect that the basic
finding of this work would change in higher dimensions: That in order to
represent a reasonable class of densities with finite kinetic energy,
distributional potentials have to be considered.

  * •

A crucial step in the proof of Hohenberg and Kohn [1] is to show that for two
potentials that differ by more than a constant, the ground states cannot be
possibly equal. In Penz  _et al._ [28] this was termed “HK2” theorem and
formulated equivalently: If two potentials share any common eigenstate then
the potentials are equal up to a constant. But Hohenberg and Kohn [1] do not
give any details about this part of the proof, which actually relies on
dividing $(V-V^{\prime})\Psi=\text{const.}\,\Psi$ by the ground-state wave
function $\Psi$. This firstly requires the wave function to be non-zero almost
everywhere, which is secured by the unique continuation property in certain
settings [32], but also it only works if the potential is a multiplication
operator. (Dreizler and Gross [33], after Eq. (2.24) who repeat the example of
Englisch and Englisch [5], comment on precisely this problem, but they put the
focus on a possible failure of the unique continuation property instead.) As
we showed here, in order to have $v$-representability, we must extend the
domain of potentials to distributions, hence this step in the proof of the
Hohenberg–Kohn theorem will generally fail. While this is not a falsification
of Hohenberg–Kohn it is still an important indication that something might go
wrong here. Similarly, in finite-lattice systems it is _known_ that there are
counterexamples to the Hohenberg–Kohn theorem while one has full
$v$-representability for all densities $0<\rho<1$ [11]. So while we gain
$v$-representability we actually might loose the unique mapping from densities
to potentials. In this context it is interesting to remember the observation
of Simen Kvaal in Teale  _et al._ [9, 2.1.17] that neither the convex
formulation of exact DFT nor the Kohn–Sham approach do in any way rely on the
Hohenberg–Kohn theorem. Rather they depend on the constraint-search functional
and the Legendre–Fenchel transformation and, in the case of Kohn–Sham,
specifically on interacting and non-interacting $v$-representability. This
means that the given result is imperative for a well-defined Kohn–Sham
procedure, since it proves the existence of the $v_{\mathrm{s}}$ potential of
the auxiliary system for the first time (cf. Mathieu Lewin in Teale  _et al._
[9, 4.5.2]). Another application is the adiabatic connection where indeed
$v$-representability for different interaction strengths needs to be assumed,
a property that can be achieved in the setting given here. In this sense it
can be argued that the achievement of $v$-representability outweighs the
potential loss of the validity of the Hohenberg–Kohn theorem.

  * •

The previous question on the validity of the Hohenberg–Kohn theorem for the
class $\mathscr{X}^{*}$ of potentials even connects directly to functional
differentiability by another result of convex analysis. Since $F(\rho)$ is
continuous at any $\rho\in\mathscr{X}_{>0}$ we have a non-empty
subdifferential $\partial F(\rho)$ (Theorem 15) and as a consequence $\rho$ is
$v$-representable. If $\partial F(\rho)$ contains just a _single_ element,
which would be guaranteed by the Hohenberg–Kohn theorem, then this implies
that $F$ is also Gateaux differentiable at $\rho$, and vice versa [27, Prop.
2.40]. Most formulations of DFT rely implicitly on this functional
differentiability. The $v$-representability is a necessary but not a
sufficient criterion for functional differentiability and the previous point
even gives some indications for a possible failure of the Hohenberg–Kohn
theorem, so the question about functional differentiability still remains
unanswered. Nevertheless we note the following interesting corollary.

######  Corollary 19.

The two statements are equivalent:

    1. (i)

Every $\rho\in\mathscr{X}_{>0}$ is represented by a _unique_ potential
$v\in\mathscr{X}^{*}$ (Hohenberg–Kohn theorem).

    2. (ii)

$F$ is Gateaux differentiable at all $\rho\in\mathscr{X}_{>0}$.

  * •

Example 2 in Section III already showed how potentials from $\mathscr{X}^{*}$
act on wave functions and densities. But it would be interesting to find ways
how to implement such potentials numerically. Since the example showed that
the potential has a semi-local effect, this brings the $v$-representing
potential close to the realm of gradient expansions [34; 35] and generalized
gradient approximations [36; 37; 38].

  * •

The main difference to other formulations of DFT is the switch from the
density space $L^{1}\cap L^{3}$ as in Lieb [4] to one with a finer topology,
the Sobolev space $H^{1}$. This means the meaning of ‘closeness’ changes, and
densities that have been nearby in the former topology, and led to non-
differentiability, are now remote, thus allowing for possible
differentiability of the functional. The reason is that the new topology
measures the energy content of the associated wave function more accurately,
as expressed through Lemma 10 and the embeddings of Lemma 8. This can have
important consequences for other problems in DFT such as the question of
convergence of the Kohn–Sham iteration scheme [39; 40; 41] where, after all,
it is always convergence with respect to a certain topology. Further and
interestingly, these Sobolev spaces appeared similarly in a recent
reformulation of the Zhao–Morrison–Parr method of density-potential inversion
in terms of Moreau–Yosida regularization Penz  _et al._ [42], although exactly
in reverse. There, the density space is $H^{-1}$ and the potentials are taken
from $H^{1}$, as it would be according to the Poisson equation $-\Delta
v=4\pi\rho$.

  * •

If one considers a Kohn–Sham system then the potential from a density-
potential inversion is usually decomposed as
$v=v_{\mathrm{ext}}+v_{\mathrm{H}}+v_{\mathrm{xc}}$. Here, the Hartree
potential fulfils (in a 3d continuous setting) the Poisson equation $-\Delta
v_{\mathrm{H}}=4\pi\rho$ and so for $\rho\in H^{1}$ it holds
$v_{\mathrm{H}}\in H^{3}$. Analogous to that the exchange-correlation
potential is sometimes considered as the effect of a fictitious ‘xc density’,
$-\Delta v_{\mathrm{xc}}=4\pi\rho_{\mathrm{xc}}$ [43]. But assuming a regular
$v_{\mathrm{ext}}$ we thus have from a general $v\in H^{-1}$ that also
$v_{\mathrm{xc}}\in H^{-1}$ which means $\rho_{\mathrm{xc}}\in H^{-3}$ which
is anything but a usual density distribution. So the current analysis
indicates that the exchange-correlation potential cannot be generally assumed
to be the effect of a fictitious charge density.

  * •

The formulation can be extended to general particle numbers by extending the
constrained search over the full Fock space. Densities integrating to a non-
integer number would not be $N$-representable anymore in a strict sense, since
there is not an integer-particle state generating that density. Yet, this
definition is also readily extended to Fock space. To differentiate from the
integer case, perhaps $\mathcal{F}$(ock)-representable density would be a
suitable name. A notable difference compared to the presented formulation is
that the constant of the potential determines the particle number to some
extent: to each integer density a closed interval of constants would be
associated and for non-integer densities a unique constant. So the potential
space would not be a simple equivalence class of $H^{-1}$, but have a more
involved equivalence class structure. This new complication could be avoided
by going all the way to a grand canonical setting with a finite temperature.
Since the entropy can by bounded by the kinetic energy [44, Appendix E], one
could still work with $\rho\in\mathscr{I}$. The entropy works as a smoothener,
similar to the Moreau–Yosida regularization [42], so might even lead to some
form of differentiability as in its one-body reduced density matrix
counterpart [45; 46].

###### Acknowledgements.

The initial ideas for this paper were born at a joint research visit of M.
Penz, M. Ruggenthaler and R. van Leeuwen in November 2021 at the Oberwolfach
Research Institute for Mathematics as a part of their “Research in Pairs”
program. S. M. Sutter and K. J. H. Giesbertz acknowledge support by the
Netherlands Organisation for Scientific Research (NWO) under Vici grant
724.017.001 and KLEIN-1 grant OCENW.KLEIN.434. R. van Leeuwen acknowledges the
Academy of Finland grant under project number 356906. M. Ruggenthaler
acknowledges the Cluster of Excellence “CUI: Advanced Imaging of Matter” of
the Deutsche Forschungsgemeinschaft (DFG), EXC 2056, project ID 390715994.

##  Appendix A Kinetic bounds for interactions

######  Lemma 20.

The two-body multiplication operator $W=\sum_{i<j}w(x_{i}-x_{j})$ with $w\in
L^{1}([-1,1])$ is kinetically bounded with relative kinetic bound $0$.

###### Proof.

Without loss of generality we only consider the spinless case and limit
ourselves to the $w(x_{1}-x_{2})$ term. The idea is to rotate the coordinates
such that the two-body interaction becomes a one-body operator and then
proceed as for the non-distributional potentials. This is achieved with the
substitution $y_{1}=(x_{1}-x_{2})/\sqrt{2}$ and $y_{2}=(x_{1}+x_{2})/\sqrt{2}$
which corresponds to the change of coordinates
$(x_{1},\ldots,x_{N})=r(y_{1},\ldots,y_{N})=((y_{1}+y_{2})/\sqrt{2},(y_{2}-y_{1})/\sqrt{2},y_{3},\ldots,y_{N})$.
The integration domain $\mathbb{T}^{N}\simeq[0,1]^{N}$ for $(x_{i})_{i}$ then
changes to $[-1/\sqrt{2},1/\sqrt{2}]\times[0,1/\sqrt{2}]\times[0,1]^{N-2}$ for
$(y_{i})_{i}$ (notation for this will be suppressed in the integrals), which
is displayed in Fig. 1.

$x_{1}$$x_{2}$$y_{2}$$y_{1}$$1$$1$

Figure 1: Change of the periodic integration domain when the $(x_{i})_{i}$
coordinates (black square) are changed to $(y_{i})_{i}$ (tilted blue
rectangle). The two domains are equivalent, since the two hatched blue (red)
areas are equivalent due to periodicity of the torus domain
$\mathbb{T}^{2}\simeq[0,1]^{2}$.

We then have

$\displaystyle\langle\Psi,w(x_{1}-x_{2})\Psi\rangle$
$\displaystyle=\int\mathrm{d}x_{1}\ldots\mathrm{d}x_{N}\lvert\Psi\rvert^{2}w(x_{1}-x_{2})$
(64) $\displaystyle=\int\mathrm{d}y_{1}\ldots\mathrm{d}y_{N}\lvert(\Psi\circ
r)(y_{1},\ldots,y_{N})\rvert^{2}w(\sqrt{2}y_{1})=\langle
w(\sqrt{2}\cdot),\rho_{r}\rangle,$

where we introduced $\rho_{r}$ as the $y_{1}$-marginal of $|\Psi\circ r|^{2}$.
Since $w\in L^{1}([-1,1])$ there is a sequence $w_{n}\in C^{\infty}([-1,1])$
such that $\|w-w_{n}\rVert_{L^{1}}\to 0$. By inserting $w=w_{n}+(w-w_{n})$,
Eq. (64) can be further estimated using the Sobolev embedding
$H^{1}(\mathbb{T})\hookrightarrow L^{\infty}(\mathbb{T})$ from Lemma 8 for
$\|\rho_{r}\|_{L^{\infty}}=\|\sqrt{\rho_{r}}\|_{L^{\infty}}^{2}$ and the
estimate in terms of the wave function from Lieb [4, Th. 1.1].

$\displaystyle|\langle w(\sqrt{2}\cdot),\rho_{r}\rangle|$
$\displaystyle\leq\|w_{n}\rVert_{L^{\infty}}\|\rho_{r}\|_{L^{1}}+\frac{1}{\sqrt{2}}\|w-w_{n}\rVert_{L^{1}}\|\rho_{r}\|_{L^{\infty}}$
(65)
$\displaystyle\leq\|w_{n}\rVert_{L^{\infty}}\|\rho_{r}\|_{L^{1}}+C\|w-w_{n}\rVert_{L^{1}}\left(\|\sqrt{\rho_{r}}\|_{L^{2}}^{2}+\|\nabla\sqrt{\rho_{r}}\|_{L^{2}}^{2}\right)$
$\displaystyle\leq
N\left(\|w_{n}\rVert_{L^{\infty}}+C\|w-w_{n}\rVert_{L^{1}}\right)\|\Psi\circ
r\rVert_{L^{2}}^{2}+CN\|w-w_{n}\rVert_{L^{1}}\|\nabla_{y_{1}}(\Psi\circ
r)\rVert_{L^{2}}^{2}$ $\displaystyle\leq
N\left(\|w_{n}\rVert_{L^{\infty}}+C\|w-w_{n}\rVert_{L^{1}}\right)\|\Psi\rVert_{L^{2}}^{2}+\frac{CN}{2}\|w-w_{n}\rVert_{L^{1}}(\|\nabla_{x_{1}}\Psi\rVert_{L^{2}}^{2}+\|\nabla_{x_{2}}\Psi\rVert_{L^{2}}^{2})$

In the last line, we just used
$\nabla_{y_{1}}=(\nabla_{x_{1}}-\nabla_{x_{2}})/\sqrt{2}$ and the triangle
inequality. We have $\|\nabla_{x_{j}}\Psi\rVert_{L^{2}}^{2}\leq 2T(\Psi)$ and
$n$ can be chosen large enough such that $\|w-w_{n}\rVert_{L^{1}}$ becomes
arbitrarily small. Thus, we conclude the zero kinetic bound. ∎

This result covers all typical interactions that depend on the distance
between two particles. The next one even extends this to distributional
potentials of the form $\mathscr{X}^{*}$.

######  Lemma 21.

A distributional interaction of the form
$W=\sum_{i<j}\nabla_{x_{i}}g(x_{i}-x_{j})$ with $g\in L^{2}([-1,1])$ is
kinetically bounded with relative kinetic bound $0$.

###### Proof.

As in the previous proof we only need to consider the spinless case and use
the same coordinate transformation
$r:(y_{1},\ldots,y_{N})\mapsto(x_{1},\ldots,x_{N})$ to get

$\nabla_{x_{1}}g(x_{1}-x_{2})=\frac{1}{2}\left(\nabla_{x_{1}}g(x_{1}-x_{2})-\nabla_{x_{2}}g(x_{1}-x_{2})\right)=\frac{1}{\sqrt{2}}\nabla_{y_{1}}g(\sqrt{2}y_{1})=:\nabla_{y_{1}}\tilde{g}(y_{1}).$
(66)

We can thus write
$\langle\Psi,(\nabla_{x_{1}}g(x_{1}-x_{2}))\Psi\rangle=\langle\nabla\tilde{g},\rho_{r}\rangle$
and then continue exactly like in the proof of Lemma 16. By introducing a
sequence $\\{\tilde{g}_{n}\\}_{n}$ in $C^{\infty}([-1,1])$ that has
$\|\tilde{g}-\tilde{g}_{n}\|_{L^{2}}\to 0$ we get an estimate

$|\langle\nabla\tilde{g},\rho_{r}\rangle|\leq\|\nabla\tilde{g}_{n}\|_{L^{\infty}}\|\rho_{r}\|_{L^{1}}+2C\|\tilde{g}-\tilde{g}_{n}\|_{L^{2}}\left(\|\rho_{r}\|_{L^{1}}+\|\nabla\sqrt{\rho_{r}}\|_{L^{2}}^{2}\right).$
(67)

From here on, the same estimates as in Eq. (65) apply which concludes the
kinetic boundedness. ∎

A well-known example for such a distributional interaction is the
$\delta$-interaction. On the torus the distribution $\delta(x)+1$ is given by
$\nabla g(x)$ where $g(x)=1-x$ which is then periodically extended to $[-1,1]$
for a proper domain of the interaction. Such a system is studied in Gaudin
[47].

## References

  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn,  _Inhomogeneous electron gas,_Phys. Rev. 136, B864 (1964).
  * Levy [1979] M. Levy,  _Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem,_Proc. Natl. Acad. Sci. 76, 6062 (1979).
  * Levy [1982] M. Levy,  _Electron densities in search of Hamiltonians,_Phys. Rev. A 26, 1200 (1982).
  * Lieb [1983] E. H. Lieb,  _Density functionals for Coulomb-systems,_Int. J. Quantum Chem. 24, 243 (1983).
  * Englisch and Englisch [1983] H. Englisch and R. Englisch,  _Hohenberg–Kohn theorem and non-V-representable densities,_Physica A Stat. Mech. Appl. 121, 253 (1983).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham,  _Self-consistent equations including exchange and correlation effects,_Phys. Rev. 140, A1133 (1965).
  * Wrighton  _et al._ [2023] J. Wrighton, A. Albavera-Mata, H. F. Rodr\́mathrm{i}guez, T. S. Tan, A. C. Cancio, J. Dufty, and S. Trickey,  _Some problems in density functional theory,_Lett. Math. Phys. 113, 41 (2023).
  * Lewin  _et al._ [2023] M. Lewin, E. H. Lieb, and R. Seiringer, in  _Density Functional Theory: Modeling, Mathematical Analysis, Computational Methods, and Applications_, edited by E. Cancès and G. Friesecke (Springer, 2023) pp. 115–182.
  * Teale  _et al._ [2022] A. M. Teale, T. Helgaker, A. Savin, C. Adamo, B. Aradi, A. V. Arbuznikov, P. W. Ayers, E. J. Baerends, V. Barone, P. Calaminici, _et al._ ,  _DFT exchange: sharing perspectives on the workhorse of quantum chemistry and materials science,_Phys. Chem. Chem. Phys. 24, 28700 (2022).
  * Chayes  _et al._ [1985] J. T. Chayes, L. Chayes, and M. B. Ruskai,  _Density functional approach to quantum lattice systems,_J. Stat. Phys. 38, 497 (1985).
  * Penz and van Leeuwen [2021] M. Penz and R. van Leeuwen,  _Density-functional theory on graphs,_J. Chem. Phys. 155, 244111 (2021).
  * Lammert [2010] P. E. Lammert,  _Well-behaved coarse-grained model of density-functional theory,_Phys. Rev. A 82, 012109 (2010).
  * Lammert [2007] P. E. Lammert,  _Differentiability of Lieb functional in electronic density functional theory,_Int. J. Quantum Chem. 107, 1943 (2007).
  * Aryasetiawan and Stott [1986] F. Aryasetiawan and M. Stott,  _Density-functional theory for two noninteracting spinless fermions,_Physical Review B 34, 4401 (1986).
  * Aryasetiawan and Stott [1988] F. Aryasetiawan and M. Stott,  _Effective potentials in density-functional theory,_Physical Review B 38, 2974 (1988).
  * Aryasetiawan [1989] F. Aryasetiawan,  _The v-representability problem in one dimension,_ personal communication (1989).
  * Chen and Stott [1991a] J. Chen and M. Stott,  _v-representability for systems of a few fermions,_Physical Review A 44, 2809 (1991a).
  * Chen and Stott [1991b] J. Chen and M. Stott,  _v-representability for systems with low degeneracy,_Physical Review A 44, 2816 (1991b).
  * Chen and Stott [1993] J. Chen and M. Stott,  _v-representability for noninteracting systems,_Physical Review A 47, 153 (1993).
  * Garrigue [2022] L. Garrigue,  _Building Kohn–Sham potentials for ground and excited states,_Arch. Ration. Mech. Anal. 245, 949 (2022).
  * Lieb and Loss [2001] E. H. Lieb and M. Loss,  _Analysis_ , 2nd ed. (American Mathematical Society, 2001).
  * Reed and Simon [1980] M. Reed and B. Simon,  _Methods of Modern Mathematical Physics. Vol. I: Functional analysis_ , 2nd ed., Vol. 1 (Academic Press, 1980).
  * Reed and Simon [1975] M. Reed and B. Simon,  _Methods of Modern Mathematical Physics. Vol. II: Fourier Analysis, Self-Adjointness_ , Vol. 2 (Academic Press, 1975).
  * Brezis [2011] H. Brezis,  _Functional Analysis, Sobolev Spaces and Partial Differential Equations_ (Springer, 2011).
  * Adams and Fournier [2003] R. A. Adams and J. J. Fournier,  _Sobolev spaces_ , 2nd ed. (Elsevier, 2003).
  * Reed and Simon [1978] M. Reed and B. Simon,  _Methods of Modern Mathematical Physics. Vol. IV: Analysis of Operators_ , Vol. 4 (Academic Press, 1978).
  * Barbu and Precupanu [2012] V. Barbu and T. Precupanu,  _Convexity and Optimization in Banach Spaces_ , 4th ed. (Springer, 2012).
  * Penz  _et al._ [2023a] M. Penz, E. I. Tellgren, M. A. Csirik, M. Ruggenthaler, and A. Laestadius,  _The structure of density-potential mapping. Part I: Standard density-functional theory,_ACS Phys. Chem. Au 3, 334 (2023a).
  * Herczyński [1989] J. Herczyński,  _On Schrödinger operators with distributional potentials,_J. Operator Theory 21, 273 (1989).
  * Clason [2020] C. Clason,  _Introduction to Functional Analysis_ (Springer Nature, 2020).
  * Singh [1983] S. Singh,  _Kronig–Penney model in reciprocal lattice space,_Am. J. Phys. 51, 179 (1983).
  * Garrigue [2018] L. Garrigue,  _Unique continuation for many-body Schrödinger operators and the Hohenberg-Kohn theorem,_Math. Phys. Anal. Geom. 21, 27 (2018).
  * Dreizler and Gross [2012] R. M. Dreizler and E. K. Gross,  _Density functional theory: an approach to the quantum many-body problem_ (Springer, 2012).
  * Langreth and Mehl [1983] D. C. Langreth and M. J. Mehl,  _Beyond the local-density approximation in calculations of ground-state electronic properties,_Phys. Rev. B 28, 1809 (1983).
  * Svendsen and von Barth [1996] P. S. Svendsen and U. von Barth,  _Gradient expansion of the exchange energy from second-order density response theory,_Phys. Rev. B 54, 17402 (1996).
  * Perdew [1985] J. P. Perdew,  _Accurate density functional for the energy: Real-space cutoff of the gradient expansion for the exchange hole,_Phys. Rev. Lett. 55, 1665 (1985).
  * Perdew [1986] J. P. Perdew,  _Density-functional approximation for the correlation energy of the inhomogeneous electron gas,_Phys. Rev. B 33, 8822 (1986).
  * Becke [1988] A. D. Becke,  _Density-functional exchange-energy approximation with correct asymptotic behavior,_Phys. Rev. A 38, 3098 (1988).
  * Penz  _et al._ [2019] M. Penz, A. Laestadius, E. I. Tellgren, and M. Ruggenthaler,  _Guaranteed convergence of a regularized Kohn–Sham iteration in finite dimensions,_Phys. Rev. Lett. 123, 037401 (2019).
  * Penz  _et al._ [2020] M. Penz, A. Laestadius, E. I. Tellgren, M. Ruggenthaler, and P. E. Lammert,  _Erratum: Guaranteed convergence of a regularized Kohn–Sham iteration in finite dimensions,_Phys. Rev. Lett. 125, 249902 (2020).
  * Lammert [2020] P. E. Lammert,  _A bivariate view of Kohn–Sham iteration and the case for potential mixing,_(5 Jul 2020) arXiv e-prints [cond-mat.mtrl-sci] 1807.06125v2 (2020), accessed 2023-11-10.
  * Penz  _et al._ [2023b] M. Penz, M. A. Csirik, and A. Laestadius,  _Density-potential inversion from Moreau–Yosida regularization,_Electron. Struct. 5, 014009 (2023b).
  * Andrade and Aspuru-Guzik [2011] X. Andrade and A. Aspuru-Guzik,  _Prediction of the derivative discontinuity in density functional theory from an electrostatic description of the exchange and correlation potential,_Phys. Rev. Lett. 107, 183002 (2011).
  * Seidl  _et al._ [2017] M. Seidl, S. Di Marino, A. Gerolin, L. Nenna, K. J. H. Giesbertz, and P. Gori-Giorgi,  _The strictly-correlated electron functional for spherically symmetric systems revisited,_(16 Feb 2017) arXiv e-prints [cond-mat.str-el] 1702.05022 (2017), accessed 2023-11-22.
  * Giesbertz and Ruggenthaler [2019] K. J. H. Giesbertz and M. Ruggenthaler,  _One-body reduced density-matrix functional theory in finite basis sets at elevated temperatures,_Phys. Rep. 806, 1 (2019).
  * Sutter and Giesbertz [2023] S. M. Sutter and K. J. H. Giesbertz,  _One-body reduced density-matrix functional theory for the canonical ensemble,_Phys. Rev. A 107, 022210 (2023).
  * Gaudin [1967] M. Gaudin,  _Un systeme a une dimension de fermions en interaction,_ Physics Letters A 24, 55 (1967).


Annamaria Ficara, Giacomo Fiumara, Pasquale De Meo and Antonio Liotta
11institutetext: University of Palermo, Palermo, Italy,
11email<EMAIL_ADDRESS>22institutetext: University of Messina, Messina,
Italy,
22email: {gfiumara<EMAIL_ADDRESS>33institutetext: University of Bozen-
Bolzano, Bolzano, Italy,
33email<EMAIL_ADDRESS>
# Correlation analysis of node and edge centrality measures in artificial
complex networks

Annamaria Ficara${}^{\textrm{\Letter}}$(0000-0001-9517-4131) 1122 Giacomo
Fiumara(0000-0003-1528-7203) 22 Pasquale De Meo(0000-0001-7421-216X) 22
Antonio Liotta(0000-0002-2773-4421) 33

###### Abstract

The role of an actor in a social network is identified through a set of
measures called centrality. Degree centrality, betweenness centrality,
closeness centrality and clustering coefficient are the most frequently used
metrics to compute the node centrality. Their computational complexity in some
cases makes unfeasible, when not practically impossible, their computations.
For this reason we focused on two alternative measures, WERW-Kpath and Game of
Thieves, which are at the same time highly descriptive and computationally
affordable. Our experiments show that a strong correlation exists between
WERW-Kpath and Game of Thieves and the classical centrality measures. This may
suggest the possibility of using them as useful and more economic replacements
of the classical centrality measures.

###### keywords:

Complex Networks $\cdot$ Social Network Analysis $\cdot$ Centrality Measures
$\cdot$ Correlation Coefficients $\cdot$ K-path

##  1 Introduction

Social network analysis (SNA) uses graph theory and networks to understand
social structures and has diverse fields of application such us psychology,
economics, criminology [1, 2, 3, 4], organizational studies and information
science.

One of the most effective SNA tools to measure social interactions has been a
simple graph which consists of a set of nodes (i.e. individuals or
organizations), often called actors in the SNA tradition, with edges between
them, also called links or connections (i.e. relational ties, e.g., friendship
relationships).

The term node centrality indicates a family of measures aimed at
characterizing the importance or the position of an actor in a network [5]. An
actor might be important for being connected to a large number of nodes or for
his power in the control over the information flowing through the network. His
absence could also result in a social network made of many isolated
components. Despite an abundance of methods for measuring centrality of
individual nodes, there are by now only a few metrics to measure centrality of
individual edges. De Meo et al. [6, 7, 8] presented a novel measure called the
$K$-path to compute link centrality. The advantage of using this metric is
that it can be computed with a near-linear time algorithm called WERW-Kpath.
Most recently, Mocanu et al. [9] developed an algorithm called Game of Thieves
which is able to compute actors and links centrality in a polylogarithmic
time. An often asked, but rarely answered, question is: are these centrality
measures correlated? [10] If there exists a high correlation between the
centrality metrics, they will have a similar behavior in statistical analyses
and for this reason the development of multiple measures could be redundant.
If there is not high correlation, these measures are unique and they can be
associated with different outcomes. Many researchers carried out studies on
the correlations between centrality measures. Valente et al. [10] investigated
the correlation among the most commonly used centrality measures and they
identified degree as the measure with the strongest overall correlations. Shao
et al. [11] uses degree to approximate closeness, betweenness, and eigenvector
in artificial and real networks. Oldham et al. [12] used $212$ different real
networks and calculated correlations between $17$ different centrality
measures showing a positive correlation. In this work, as in a previous study
[13], we try to answer to a different question: are these centrality measures
correlated with Game of Thieves? If so, we can use Game of Thieves thus
considerably reducing the execution time in the computation of node and edge
centrality in very large networks. We have done a correlation analysis using
the most common types of coefficients which are Pearson, Spearman and Kendall
and three kinds of artificial complex networks: scale-free, small-world and
Erdös-Rényi random networks. In a previous work [13], we analyzed the
correlation coefficients among node centrality metrics when the number of
nodes in both artificial and real networks increased using GoT only to measure
an actor centrality. In this study, we have made an analysis on artificial
networks taking into account the increase of the number of links when the
network size does not change. We also used GoT to compute edges centrality
comparing it with the WERW-Kpath. The paper is planned as follows. Sect. 2
presents a brief description of the main theoretical definitions of the
considered centrality metrics and correlation coefficients. Sect. 3 contains a
discussion of our methods and the obtained results. Sect. 4 shows the
conclusions and future developments of our study.

##  2 Background

###  2.1 Basic Definitions

A graph $G$ is an ordered pair $G=(N,E)$ where $N$ and $E$ are the sets of
nodes and edges respectively. A graph is called undirected if edges are
bidirectional. The main node centrality measures are:

  * •

degree centrality (DC) [14], which is defined as:

$DC_{i}=\frac{d_{i}}{n-1}\,,$ (1)

where $d_{i}$ is the number of edges of a node $i$ (i.e. the degree of $i$)
and $n$ is the network size.

  * •

betweenness centrality (BC) [15], which is given by the following formula:

$BC_{i}=\sum\limits_{h,k}\frac{v^{i}_{hk}}{g_{hk}}\,,$ (2)

where $v^{i}_{hk}$ is the number of shortest paths from a node $h$ to a node
$k$ by passing through $i$ and $g_{hk}$ is the total number of shortest paths
from $h$ to $k$.

  * •

closeness centrality (CL) [14], which is defined as:

$CL_{i}=\frac{n}{\sum\limits_{j}d_{ij}}\,,$ (3)

where $d_{ij}$ is the distance between $i$ and $j$ and $n$ is the network
size.

  * •

clustering coefficient (CC) [16], which is defined as the number of complete
triangles $T_{i}$ in which a node $i$ participates divided by the maximum
possible number of these triangles:

$CC_{i}=\frac{2T_{i}}{d_{i}(d_{i}-1)}\,.$ (4)

A more recent centrality measure is called Game of Thieves (GoT) [9]. The game
proceeds in epochs. When it begins, each node has a certain number of thieves
and virtual diamonds or vdiamonds. At each epoch $e$, a thief located on a
node $i$ randomly picks a neighbor of $i$. He moves to this new node and, if
he finds a vdiamond, he fetches it. Then, he brings the vdiamond back to his
home node. After $T$ epochs, the centrality of each node $i$ is computed as
the average number of vdiamonds present on $i$:

$\Phi_{T}^{i}=\frac{1}{T}\sum\limits_{e=0}^{T}\Phi_{e}^{i}\,.$ (5)

The centrality of each link $l$ is computed as the average number of thieves
who carry a vdiamond passing through $l$:

$\Psi^{l}_{T}=\frac{1}{T}\sum\limits_{e=0}^{T}\Psi_{e}^{l}\,.$ (6)

An other recent measure of link centrality for social networks is the K-path
which is defined as:

$L^{k}(l)=\sum\limits_{i}\frac{\sigma_{i}^{k}(l)}{\sigma_{i}^{k}}\,,$ (7)

where $\sigma_{i}^{k}(l)$ is the number of $k$-paths originating from a node
$i$ and traversing the link $l$ and $\sigma_{i}^{k}$ is the number of
$k$-paths which originate from $i$. A near linear time algorithm called WERW-
Kpath (Weighted Edge Random Walks – $K$ Path) [6, 7, 8] is able to compute
this centrality index.

Table 1: Centrality measures. $n$ is the cardinality of $N$ and $m$ is the cardinality of $E$ in a graph $G(N,E)$. Measure | Centrality | Computational Complexity |
---|---|---|---
DC | Nodes | $O(m)$ |
BC | Nodes |  $O(mn)$ or $O(n^{3})$ |
CL | Nodes | $O(n^{3})$ |
CC | Nodes | $O(n^{2})$ |
GoT | Nodes and Edges |  $O(log^{2}n)$ or $O(log^{3}n)$ |
WERW-Kpath | Edges | $O(km)$ |

Table 1 shows how GoT represents a great step forward in terms of
computational complexity with respect to classical algorithms of node
centrality such us BC, CL and CC which have a quadratic computational
complexity. The DC has a linear time complexity but GoT still remains a better
option because it is also able to compute the centrality of edges.

###  2.2 Correlation Coefficients

Given two random variables $a$ and $b$, the Pearson’s $r$ correlation
coefficient [17] is defined as:

$r=\frac{cov(a,b)}{\sigma_{a}\sigma_{b}}\,,$ (8)

where $cov(a,b)$ is the covariance of $a$ and $b$ and $\sigma_{a}\sigma_{b}$
is the product of their standard deviations.

The Spearman’s $\rho$ rank correlation coefficient [18] is defined as the
Pearson’s $r$ between the rank values of the two variables. For a sample of
size $s$, the $s$ raw scores $a$ and $b$ are converted to ranks $rg_{a}$ and
$rg_{b}$, and $\rho$ is computed as:

$\rho=\frac{cov(rg_{a},rg_{b})}{\sigma_{rg_{a}}\sigma_{rg_{b}}}\,,$ (9)

where $cov(rg_{a},rg_{b})$ is the covariance of the rank variables and
$\sigma_{rg_{a}}\sigma_{rg_{b}}$ are the standard deviations of the rank
variables.

Given two samples $a$ and $b$, where each sample size is $s$, Kendall’s $\tau$
rank correlation coefficient [19] is defined according to the following
formula:

$\tau=\frac{s_{c}-s_{d}}{\frac{1}{2}s(s-1)}\,,$ (10)

where $s_{c}$ and $s_{d}$ are the numbers of concordant and discordant pairs
respectively and $\frac{1}{2}s(s-1)$ is the total number of pairings with $a$
and $b$.

##  3 Methodology and Results

We investigated the correlations among the node and edge centrality measures,
described in Subsect. 2.1, in different kinds of artificial networks. The
network models include the scale-free (SF) [20], the small-world (SW) [21] and
the Erdös-Rényi (ER) [22] random networks. For each class, we randomly
generated undirected and unweighted networks with $10,000$ nodes and between
$5,000$ and $50,000$ edges.

We used the model proposed by Holme and Kim [23] to generate the SF networks.
In each experiment, we chose the number of nodes $n=10,000$, we added
$m=\\{5,15,25,35,50\\}$ random edges for each new node $i$, and we picked a
probability $p=0.3$ of adding a triangle after adding a random edge.

SW graphs are generated using Newman-Watts-Strogatz (NWS) small-world model
[21]. In each experiment, we chose the number of nodes $n=10,000$,
$k=\\{6,18,32,64\\}$ neighbors with which connect each node $i$ in the ring
topology, and a probability $p=0.6$ of rewiring each edge.

ER networks are generated choosing the network size $n=10,000$, and for edge
creation the probability values $p=\\{0.001,0.003,0.005,0.010\\}$. For the
implementation of the artificial networks and the centrality metrics such as
DC, BC, CL, and CC, we used Python and NetworkX module [24]. To run the WERW-
Kpath111Available at http://www.emilio.ferrara.name/code/werw-kpath/ [6, 7, 8]
we set the value of the random walk to $k=10$. For GoT222Available at
http://github.com/dcmocanu/centrality-metrics-complex-networks we set $1$
thief per node and the initial amount of vdiamonds per node equal to the
network size. We let GoT to run for $T=\log^{3}n$ epochs. In Fig. 1 are shown
the results of the Pearson’s $r$ correlation coefficient (left column), the
Spearman’s $\rho$ (middle column) and the Kendall’s $\tau$ rank correlation
coefficients (right column) for SF networks (upper row), for SW networks
(middle row) and ER random networks (lower row). In SF networks (see Fig. 1
upper row), GoT and DC have the strongest negative correlation. GoT exhibits a
strong negative correlation with the BC and slightly weaker with the CL. There
is a positive correlation between GoT and the CC which becomes stronger for
the rank correlation coefficients. In particular, we can observe that the
Sperarman’s $\rho$ is always larger than the Kendall’s $\tau$. Except for the
DC, we can observe a small deviation of the rank correlation coefficients when
the number of edges in the analyzed networks is smaller which is not visible
anymore when edges grow big enough. In SW networks (see Fig. 1 middle row),
there is an almost identical strong negative correlation among GoT and DC, BC
and CL. We can also still observe a strong positive correlation between GoT
and the CC. All of these correlations become weaker when the number of edges
grows. In ER networks (see Fig. 1 lower row), the correlation among GoT and
DC, BC and CL is the same as in SW networks. But, in this kind of networks,
there is no correlation between GoT and the CC even if we can observe a very
little deviation when the number of edges is small. In all networks, we can
observe a strong positive correlation among GoT and WERW-Kpath which remains
constant when the number of edges grows.

Figure 1: Pearson’s $r$ (a), Spearman’s $\rho$ (b) and Kendall’s $\tau$ (c)
correlation coefficients.

##  4 Conclusions

In this work we used artificial networks, i.e. SF, SW and ER networks, to
examine the correlation between well known and more recently proposed
centrality measures. A strong correlation implies the possibility of
approximating the metric with the highest computational complexity using the
other. We have done an analysis observing the Pearson, Spearman and Kendall
correlation coefficients on artificial networks by changing the amount of
links given a fixed network size. Our results show a strong correlation among
GoT and the most known centrality algorithms. For this reason, it can
substitute them in the computation of node centrality in large networks. GoT
can also replace the CC as a property of SW and SF networks. A future
development of our work will be to focus more on the application of WERW-
Kpath, changing the values of $k$, and GoT to the case of edge centrality. The
two measures seems to be equivalent and it could be also interest to make a
deeper analysis on their computational costs on large networks.

## References

  * [1] Ficara, A., Cavallaro, L., De Meo, P., Fiumara, G., Catanese, S., Bagdasar, O., Liotta, A.: Social network analysis of Sicilian Mafia interconnections.  In: H. Cherifi, S. Gaito, J.F. Mendes, E. Moro, L.M. Rocha (eds.) Complex Networks and Their Applications VIII, pp. 440–450. Springer International Publishing, Cham (2020).  10.1007/978-3-030-36683-4_36
  * [2] Calderoni, F., Catanese, S., De Meo, P., Ficara, A., Fiumara, G.: Robust link prediction in criminal networks: A case study of the Sicilian Mafia.  Expert Systems with Applications 161, 113,666 (2020).  10.1016/j.eswa.2020.113666
  * [3] Cavallaro, L., Ficara, A., De Meo, P., Fiumara, G., Catanese, S., Bagdasar, O., Song, W., Liotta, A.: Disrupting resilient criminal networks through data analysis: The case of Sicilian Mafia.  PLOS ONE 15(8), 1–22 (2020).  10.1371/journal.pone.0236476
  * [4] Cavallaro, L., Ficara, A., Curreri, F., Fiumara, G., De Meo, P., Bagdasar, O., Liotta, A.: Graph comparison and artificial models for simulating real criminal networks.  In: R.M. Benito, C. Cherifi, H. Cherifi, E. Moro, L.M. Rocha, M. Sales-Pardo (eds.) Complex Networks and Their Applications IX, pp. 286–297. Springer International Publishing, Cham (2021)  10.1007/978-3-030-65351-4_23
  * [5] Wasserman, S., Faust, K., Granovetter, M., Iacobucci, D.: Social Network Analysis: Methods and Applications.  Structural Analysis in the Social Sciences. Cambridge University Press (1994)
  * [6] De Meo, P., Ferrara, E., Fiumara, G., Ricciardello, A.: A novel measure of edge centrality in social networks.  Knowledge Based Systems 30, 136–150 (2012).  10.1016/j.knosys.2012.01.007
  * [7] De Meo, P., Ferrara, E., Fiumara, G., Provetti, A.: Enhancing community detection using a network weighting strategy.  Information Sciences 222, 648–668 (2013).  10.1016/j.ins.2012.08.001
  * [8] De Meo, P., Ferrara, E., Fiumara, G., Provetti, A.: Mixing local and global information for community detection in large networks.  Journal of Computer and System Sciences 80(1), 72–87 (2014).  10.1016/j.jcss.2013.03.012
  * [9] Mocanu, D.C., Exarchakos, G., Liotta, A.: Decentralized dynamic understanding of hidden relations in complex networks.  Scientific Reports 8(1), 1571 (2018).  10.1038/s41598-018-19356-4
  * [10] Valente, T.W., Coronges, K., Lakon, C., Costenbader, E.: How correlated are network centrality measures?  Connections (Toronto, Ont.) 28(1), 16–26 (2008)
  * [11] Shao, C., Cui, P., Xun, P., Peng, Y., Jiang, X.: Rank correlation between centrality metrics in complex networks: An empirical study.  Open Physics 16(1), 1009–1023 (2018).  10.1515/phys-2018-0122
  * [12] Oldham, S., Fulcher, B., Parkes, L., Arnatkeviciute, A., Suo, C., Fornito, A.: Consistency and differences between centrality measures across distinct classes of networks.  PLOS ONE 14(7), 1–23 (2019).  10.1371/journal.pone.0220061
  * [13] Ficara, A., Fiumara, G., De Meo, P., Liotta, A.: Correlations among Game of Thieves and other centrality measures in complex networks.  In: G. Fortino, A. Liotta, R. Gravina, A. Longheu (eds.) Data Science and Internet of Things. Springer International Publishing (2021).  10.1007/978-3-030-67197-6_3
  * [14] Freeman, L.C.: Centrality in social networks conceptual clarification.  Social Networks 1(3), 215–239 (1978).  10.1016/0378-8733(78)90021-7
  * [15] Brandes, U.: On variants of shortest-path betweenness centrality and their generic computation.  Social Networks 30(2), 136–145 (2008).  10.1016/j.socnet.2007.11.001
  * [16] Saramäki, J., Kivelä, M., Onnela, J.P., Kaski, K., Kertész, J.: Generalizations of the clustering coefficient to weighted complex networks.  Phys. Rev. E 75, 027,105 (2007).  10.1103/PhysRevE.75.027105
  * [17] Chen, P., Popovich, P.: Correlation: parametric and nonparametric measures.  Sage university papers series. No. 07-139. Sage Publications (2002)
  * [18] Spearman, C.: General intelligence, objectively determined and measured.  The American Journal of Psychology 15(2), 201–292 (1904).  10.2307/1412107
  * [19] Kendall, M., Gibbons, J.: Rank Correlation Methods.  Charles Griffin Book. E. Arnold (1990)
  * [20] Barabási, A.L., Albert, R.: Emergence of scaling in random networks.  Science 286(5439), 509–512 (1999).  10.1126/science.286.5439.509
  * [21] Newman, M., Watts, D.: Renormalization group analysis of the small-world network model.  Physics Letters A 263(4), 341–346 (1999).  10.1016/S0375-9601(99)00757-4
  * [22] Erdös, P., Rényi, A.: On Random Graphs I.  Publicationes Mathematicae Debrecen 6, 290 (1959)
  * [23] Holme, P., Kim, B.J.: Growing scale-free networks with tunable clustering.  Phys. Rev. E 65, 026,107 (2002).  10.1103/PhysRevE.65.026107
  * [24] Hagberg, A.A., Schult, D.A., Swart, P.J.: Exploring network structure, dynamics, and function using networkx.  In: G. Varoquaux, T. Vaught, J. Millman (eds.) Proceedings of the 7th Python in Science Conference, pp. 11–15. Pasadena, CA USA (2008)