is referred to, e.g., the essay by the famous mathematical physicist John Baez
[121].

###### The resolution: regularization and renormalization

Nevertheless, we need to “cure” this inconsistency for our theory and the
usual way to do so is by a so-called _regularization_ procedure for short
distances. This means in practice that we consider a small sphere instead of a
point charge, which directly removes the divergence in Eq. (1.23). Still, we
need to decide about the exact shape and size of this sphere, which is
formally done by introducing a _form-factor_ in the equations.484848There are
different ways to do this and certain connected issues. See [9, part 2/part 4]
for details By the introduction of the cutoff, we remove small length scales
from the description that “come from a more refined theory” [9, p. 15]. The
regularization removes the divergence of the self-energy, i.e., electrostatic
energy of the particles corresponding to the Coulomb force (1.23). But it
still contributes to the energy-momentum relation of the coupled theory, which
consequently differs from the one of a free particle. Since the latter is an
experimental fact, we need to “fix” the coupled energy-momentum relation by
performing a _renormalization_. We see here that although mostly discussed in
the area of quantum field theory, the technique of renormalization is already
required on the classical level. In the case of Lorentz-Maxwell theory, this
means that we choose the particle’s mass as a parameter that we utilize to fit
the theory to the experiment [9, part 6 and part 7].

We see that the description of coupled light-matter systems is a very
challenging topic, even on the “simplest,” i.e., the classical, level and
there are many other subtleties and interesting issues.

###### The quantization of the classical theory: gauge dependence

Starting from the classical Lorentz-Maxwell theory, we now introduce NR-QED by
the correspondence principle. The standard way to do this is to employ the
canonical quantization procedure, which requires the _Hamiltonian_ formulation
of the classical theory. Regarding the matter part, this is straightforward,
but for the Hamiltonian formulation of Maxwell’s theory, care has to be taken.
To describe the electromagnetic field in terms of canonical variables, one
usually introduces the vector and scalar potential $\mathbf{A}$ and $\phi$.
These are connected to the fields $\mathbf{E},\mathbf{B}$ by

$\displaystyle\mathbf{E}(\mathbf{r},t)$
$\displaystyle=-\tfrac{1}{c}\partial_{t}\mathbf{A}(\mathbf{r},t)-c\nabla\phi(\mathbf{r},t)$
(1.24a) $\displaystyle\mathbf{B}(\mathbf{r},t)$
$\displaystyle=\tfrac{1}{c}\nabla\times\mathbf{A}(\mathbf{r},t),$ (1.24b)

which is not a unique correspondence but leaves a so-called gauge freedom.
This means that the quantization procedure _depends_ on the gauge [10]. In the
nonrelativistic description, it is very convenient to employ the so-called
_Coulomb gauge_

$\displaystyle\nabla\cdot\mathbf{A}(\mathbf{r},t)=0.$ (1.25)

The Coulomb gauge removes the longitudinal part of the vector potential, which
by means of Gauss’s law (Eq. (1.20c)) also fixes the scalar potential
$\phi(\mathbf{r},t)$.494949To see this, compare Eq. (1.20c) with the
definition of $\mathbf{E}$ in terms of $bA,\phi$, cf. Eq. (1.24a):
$\rho=\nabla\cdot\mathbf{E}=\nabla\cdot(-1/c\partial_{t}\mathbf{A}-c\nabla\phi)=-c\nabla^{2}\phi$.
The remaining two transversal degrees of freedom of
$\mathbf{A}=\mathbf{A}_{\perp}$ describe the two known polarizations of the
electromagnetic field. This means that in the Coulomb gauge, there is no
“superfluous” field component and thus it can be seen as kind of maximal
gauge. That is why basically all derivations of NR-QED are done in the Coulomb
gauge. However, there are also other gauges that are maximal in the same
sense, e.g., the Poincare gauge and employing these might be advantageous in
certain scenarios [11].

###### The Hamiltonian of NR-QED: how to construct a well-defined mathematical
theory of light-matter interaction

However, to discuss the quantization of the Maxwell-Lorentz theory we stick to
the standard case. We employ Coulomb gauge and follow the prescription of the
canonical quantization, which leads the so-called _Pauli-Fierz Hamiltonian_.
Since this is a standard procedure that can be found in several textbooks, we
only summarize the most important aspects. We again follow closely the book of
Herbert Spohn [9]. The Pauli-Fierz Hamiltonian describes the quantum dynamics
of $N$ charged particles that are coupled to the (vacuum) electromagnetic
field. Since we are interested in molecular systems, we further divide
$N=N_{e}+N_{n}$ in $N_{e}$ electrons and $N_{n}$ nuclei. The quantum
formulation does not remedy the inconsistencies of the classical theory, but
adds further problems. The most important property of a Hamiltonian of a
quantum theory is hermiticity (which in mathematical literature is usually
called _self-adjointness_). Importantly, it is possible to prove the self-
adjointness of the Pauli-Fierz Hamiltonian under the following two conditions
[9, part 13.3].

  1. 1.

We need to remove the self-energy of the particles analogously to the
classical case.

  2. 2.

We need to introduce a further “suitable” ultraviolet cutoff. This comes _on
top_ of the classical cutoff that we already included by quantizing extended
instead of point charges. Still, the new cutoff just further renormalizes the
mass and in most situations, we can include it as in the classical case by
using the physical instead of the bare mass of the charge carriers. However,
although the proof for self-adjointness provides a definition of the cutoff,
it remains arbitrary up to a certain degree.

A Hamiltonian that obeys these conditions is the Pauli-Fierz Hamiltonian of
Def. 2

List of defenvs 2 Pauli-Fierz Hamiltonian

A system consisting of $N_{e}$ electrons and $N_{n}$ nuclei with masses
$m_{i}$ and charges $q_{i},i=1,...,N_{n}+N_{e}$ that are minimally coupled to
the electromagnetic field in the Coulomb gauge are described in the
nonrelativistic energy regime by the _Pauli-Fierz Hamiltonian_
$\displaystyle\begin{split}\hat{H}_{PF}=&\sum_{i}^{N_{e}+N_{n}}\tfrac{1}{2m_{i}}\left(-\mathbf{i}\hbar\nabla_{i}-\frac{q_{i}}{c}\hat{\mathbf{A}}_{\perp}(\mathbf{r}_{i})\right)^{2}+\frac{1}{8\pi\epsilon_{0}}\sum_{i\neq
j}^{N_{e}+N_{n}}\frac{q_{i}q_{j}}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}\\\
&+\frac{1}{2}\sum_{\mathbf{k},s}\left((-\mathbf{i}\partial_{q_{\mathbf{k},s}})^{2}+\omega_{k}^{2}q_{\mathbf{k},s}^{2}\right)\end{split}$
(1.26) Here $\mathbf{i}$ is the imaginary unit, $\hbar$ the Planck constant
and $\epsilon_{0}$ the vacuum permittivity. We denoted the field operator of
the (transversal) vector potential with $\hat{\mathbf{A}}_{\perp}$ that reads
$\displaystyle\hat{\mathbf{A}}_{\perp}(\mathbf{r})=\sqrt{\tfrac{\hbar
c^{2}}{\epsilon_{0}}}\sum_{\mathbf{k},s}\frac{1}{\omega_{k}}\left(\hat{a}_{\mathbf{k},s}^{\vphantom{\dagger}}\mathbf{S}_{\mathbf{k},s}(\mathbf{r})+\hat{a}_{\mathbf{k},s}^{\dagger}\mathbf{S}^{*}_{\mathbf{k},s}(\mathbf{r})\right),$
(1.27) where
$\displaystyle\mathbf{S}_{\mathbf{k},s}(\mathbf{r})=\tfrac{1}{\sqrt{V}}\mathbf{n}_{\mathbf{k},s}e^{\mathbf{i}\mathbf{k}\cdot\mathbf{r}}$
(1.28) are the mode functions for a plane wave with wave-vector $\mathbf{k}$
and frequency $\omega_{k}=ck$ ($k\equiv|\mathbf{k}|$). We denote with
$V=a^{3}$ the volume of the quantization box with side-length $a$.
Furthermore, $\mathbf{n}_{\mathbf{k},s}$ denotes the two polarization vectors
($s=1,2$) for every mode. To arrive at the form (1.26), we have defined the
displacement coordinates $\displaystyle q_{\mathbf{k},s}=$
$\displaystyle\sqrt{\tfrac{\hbar}{2\omega_{k}}}(\hat{a}_{\mathbf{k},s}^{\dagger}+\hat{a}_{\mathbf{k},s}^{\vphantom{\dagger}})$
(1.29) $\displaystyle-\mathbf{i}\partial_{q_{\mathbf{k},s}}=$
$\displaystyle\mathbf{i}\sqrt{\tfrac{\hbar\omega_{k}}{2}}(\hat{a}_{\mathbf{k},s}^{\dagger}-\hat{a}_{\mathbf{k},s}^{\vphantom{\dagger}}).$
(1.30)

###### A remark on the validity of Pauli-Fierz theory

We want to conclude this subsection with a small remark on the interpretation
Pauli-Fierz theory. This paragraph is not necessary to follow the course of
this text and can be skipped by the fast reader.

Pauli-Fierz theory as defined in Def. 2 is mathematically a well-defined
quantum theory, i.e., the Hamiltonian of definition 2 is self-adjoint and it
has eigenstates that can be described by a many-particle wave function.
However, to remove the occurring divergences of the theory, we have to
introduce a cutoff that is in principle arbitrary. Thus, for high-energy
processes the accuracy of the theory is somewhat unclear. Nevertheless, the
range of validity, on which these mathematical issues are not problematic, is
huge:

> The claimed range of validity of the Pauli-Fierz Hamiltonian is
> flabbergasting… As the bold claim goes, any physical phenomenon in between
> [gravity on the Newtonian level and nuclear- and high-energy physics],
> including life on Earth, is accurately described through the Pauli-Fierz
> Hamiltonian.  (Spohn [9, p. 157])

The conclusion that Pauli-Fierz theory provides a useful description of, e.g.,
molecules as complex as DNA or whole organisms is indeed “flabbergasting.”
Thus, it seems that Dirac’s vision has become true (see the quotation in the
beginning of Sec. 1.3.1): we have found a mathematical theory that describes
“a large part of physics and the whole of chemistry.” We just need to solve
this theory by “approximate practical methods of applying quantum mechanics
[…], which can lead to an explanation of the main features of complex atomic
systems without too much computation.”

Interestingly, although Pauli-Fierz theory is known since the 1930s, no
practical method has yet been developed that is capable to solve the theory
for any relativistic scenarios. Spohn remarks on this reductionist’s
perspective:

> Of course, our trust is not based on strict mathematical deductions from the
> Pauli-Fierz Hamiltonian. This is too difficult a program. Our confidence
> comes from well-studied limit cases. […] In the static limit we imagine
> turning off the interaction to the quantized part of the Maxwell field. This
> clearly results in Schrödinger particles interacting through a purely
> Coulombic potential, for which many predictions are accessible to
> experimental verification. But beware, even there apparently simple
> questions remain to be better understood. For example, the size of atoms as
> we see them in nature remains mysterious if only the Coulomb interaction and
> the Pauli exclusion principle are allowed.  (Spohn [9, p. 157])

Another example for such a limit is the standard form of cavity QED that we
will derive in the next section from Pauli-Fierz theory by performing the
long-wavelength approximation. We will see in the course of this thesis, that
only finding the (approximate) ground state of cavity QED is an utmost
difficult task. Without the great knowledge of the many-electron problem that
we have nowadays, this would probably be an impossible task. This explicitly
shows the limitation of the very common (reductionist’s) point of view that
defines a hierarchy between more general theories and their limit cases.
Although Pauli-Fierz theory is more fundamental in the sense that it is more
general than, e.g., many-electron Schrödinger theory, it is the predictions of
the lower member that are the most important justification for the higher
lying member.

This illustrates how long the way from writing down a Hamiltonian until
describing nature really is. And most importantly, it is not “just
mathematics” that is required to go this way as the reductionist’s hypothesis
suggest. This has been pointed out already in 1972 by Anderson [122] with
several examples from molecular and solid-state physics. The impressive
success of the cavity-QED models is another examples that prove this far too
simple idea wrong: the power of these models lies especially in their validity
in settings, where the explicit mathematical derivation from first-principles,
i.e., the Pauli-Fierz Hamiltonian is not possible (see Sec. 1.2). Spohn
summarizes this in the following way:

> [If t]he Hamiltonian is a self-adjoint linear operator, [… o]f course this
> does not mean that we have solved any physical problem. It just assures us
> of a definite mathematical framework within which consequences can be
> explored.  (Spohn [9, p. 158])

This is the reason for the long discussion of the mathematical foundations of
Pauli-Fierz theory in this subsection. To develop new methods that accurately
describe electron-photon interaction from first principles, we need a
(mathematically and physically) firm ground. With the corresponding methods,
we will however not “solve” Pauli-Fierz theory in the sense that we
diagonalize $\hat{H}_{PF}$. There are no analytical solutions known and even
the simplest case of one particle, coupled to the Maxwell-field is basically
numerically inaccessible. But we need such a “definite mathematical framework”
to derive approximate theories.

#####  1.3.3 NR-QED in the Long-Wavelength Limit

In this subsection, we derive the cavity-QED Hamiltonian from Pauli-Fierz
theory. The physical setting, defined by this limit is the basis for our
analysis in the following chapters. Since we aim to understand strong-coupling
phenomena that require the strong mode confinement of cavities, we constrain
the electromagnetic field description to only a _few modes_. These modes are
assumed to have long wavelengths in comparison to the spatial extension of the
considered matter systems, such that the _dipole approximation_ is
valid.505050Note that the dipole approximation is accurate also in more
general settings in the context of cavity QED, e.g., nanoplasmonic
environments. Regarding the matter-description, we only consider the
(adiabatic) electronic structure in the electrostatic potential of
(quasi)static or clamped nuclei. This means that we constrain our description
to settings, where the _Born-Oppenheimer_ approximation can be applied and
disregard the dynamics of the nuclei and their quantum nature. This setting is
however sufficient to describe a large range of chemical processes
[123].515151See also Ch. 2.

Despite these many simplifications, we are still confronted with an
extraordinarily difficult problem. Also in this limit, there are currently
only two special cases, for which analytical solutions are known.525252The
first considers one electron in a harmonic potential [9, part13.7] and the
second $N$ electrons without interaction in a box with periodic boundary
conditions [124].

###### The starting point: the macroscopic description of Pauli-Fierz theory

The starting point of the derivation is the non-relativistic theory of quantum
electrodynamics, which is described by the Pauli-Fierz Hamiltonian
$\hat{H}_{PF}$ (Def. 2). Since we are interested in processes with subatomic
length and energy scales, we use from now on a suitable unit system, so-called
atomic units (see Def. 3)

List of defenvs 3 Atomic units

The system of _atomic units_ (denoted by _[a.u.]_) is defined by the following
base units: mass $m_{e}$ the electron mass charge $e$ the electronic charge
action $\hbar$ the reduced Planck constant permittivity $4\pi\epsilon_{0}$ the
inverse Coulomb constant This means that we measure mass in multiples of the
electron mass etc. This is formally done by setting these four fundamental
constants to one. Derived from these base units, we have the following derived
units (we list only the most important ones) length $$\mathrm{\text{ $a_{0}$
}}$=4\pi\epsilon_{0}\hbar^{2}/(m_{e}e^{2})$ $1$\mathrm{\text{ $a_{0}$
}}$=0.529$\mathrm{\SIUnitSymbolAngstrom}$$ “bohr” energy $$\mathrm{\text{
$E_{\textup{ h }}$ }}$=\hbar^{2}/(m_{e}a_{0}^{2})$ $1$\mathrm{\text{
$E_{\textup{ h }}$ }}$=27.211$\mathrm{eV}$$ “Hartree” Furthermore, we can
derive the numerical value of the speed of light $c=1/\alpha\approx 137$
_[a.u.]_ , where we denoted with $\alpha$ the fine-structure constant
$\alpha=e^{2}/(4\pi\epsilon_{0}\hbar c)$.

We start the derivation with the _macroscopic description_ of Maxwell’s
equations, which separates the charge current

$\displaystyle\mathbf{j}(\mathbf{r},t)=\mathbf{j}^{b}(\mathbf{r},t)+\mathbf{j}^{f}(\mathbf{r},t)$

into its bound part $\mathbf{j}^{b}$ and the free current $\mathbf{j}^{f}$.
The former is generated by the magnetization $\mathbf{M}$ and polarization
$\mathbf{P}$ of the matter system as

$\displaystyle\mathbf{j}^{b}(\mathbf{r},t)=\nabla\times\mathbf{M}(\mathbf{r},t)+\partial_{t}\mathbf{P}(\mathbf{r},t).$

We can then define the displacement field
$\mathbf{D}=\epsilon_{0}\mathbf{E}+\mathbf{P}$ and the magnetization field
$\mathbf{H}=\frac{1}{\mu_{0}}\mathbf{B}-\mathbf{M}$, which are generated by
the free current as

$\displaystyle\mathbf{j}^{f}(\mathbf{r},t)=-\nabla\times\mathbf{H}(\mathbf{r},t)+\partial_{t}\mathbf{D}(\mathbf{r},t).$

This construction is in a sense artificial, because the actual forces that are
exerted by the field are only due to the total electric and magnetic fields
$\mathbf{E},\mathbf{B}$. But it is very suited for the following theoretical
considerations, because it separates out the part of the electromagnetic field
that is exclusively produced by the (bound) charges of our system.

The connection to the macroscopic form of (1.26) is given by the unitary
Power-Zienau-Woolley (PZW) transformation

$\displaystyle\hat{U}_{PZW}=\exp\left(-\mathbf{i}\alpha\int\mathrm{d}\mathbf{r}\hat{\mathbf{P}}_{\perp}(\mathbf{r})\cdot\hat{\mathbf{A}}_{\perp}(\mathbf{r})\right),$

where we only need to transform the transversal fields, because of the Coulomb
gauge. With the choice of our matter system, we can explicitly define the
polarization

$\displaystyle\mathbf{P}(\mathbf{r})=-\sum_{i=1}^{N_{e}}\mathbf{r}_{i}\int_{0}^{1}\delta^{3}(\mathbf{r}-s\mathbf{r}_{i})\mathrm{d}s+\sum_{i=1}^{N_{n}}Z_{i}\mathbf{R}_{i}\int_{0}^{1}\delta^{3}(\mathbf{R}-s\mathbf{R}_{i})\mathrm{d}s,$
(1.31)

where we separated electronic $\mathbf{r}_{i}$ and nuclear coordinates
$\mathbf{R}_{i}$ and denoted the charge of the $i$-th nucleus by $Z_{i}$.
Since we included all physical charges in the definition of $\mathbf{P}$, we
trivially have $\mathbf{j}^{f}=0$.535353This follows, because we treat a
closed system here. In a more general setting, the free current could, e.g.,
represent an external current that pumps the cavity mode.

###### The dipole approximation

We can now introduce the dipole approximation by disregarding the integration
over s in the definition of $\mathbf{P}$ (Eq. (1.31)) or equivalently by
setting

$\displaystyle\mathbf{A}(\mathbf{r})\approxeq\mathbf{A}(0),$ (cf. 1.6)

where we assumed that the coordinate frame origin is located at the center of
charge. Specifically, this leads to the transformation

$\displaystyle\hat{U}_{PZW,dip}=\exp\left(-\mathbf{i}\sqrt{4\pi}\sum_{\mathbf{k},s}S_{\mathbf{k},s}(0)\mathbf{n}_{\mathbf{k},s}\cdot\mathbf{R}\right)q_{\mathbf{k},s},$
(1.32)

where $\mathbf{n}_{\mathbf{k},s}$ is the polarization vector and
$S_{\mathbf{k},s}(0)\propto 1/\sqrt{V}$ the (absolute value of the) mode
function of the photon mode with wave vector $\mathbf{k}$ and polarization $s$
(see Def. 2). Additionally, we have defined the total dipole of the matter
system

$\displaystyle\mathbf{R}=-\sum_{i=1}^{N_{e}}\mathbf{r}_{i}+\sum_{i=1}^{N_{n}}Z_{i}\mathbf{R}_{i}.$
(1.33)

Further, we constrain the description to a finite number of $M$ photon modes,
where $M$ cannot be chosen arbitrarily, since the dipole approximations breaks
down for very large $M$ [124]. In the following, we will usually assume that
$M$ is very small (corresponding to the enhanced modes of an optical cavity)
and thus, this issue will not be problematic. Note that Eq. (1.32) is also
called _length-gauge_ transformation and it is the generalized form of the
transformation $U_{d}$ (defined above Eq. (1.2.1)) that was employed in the
derivation of the cavity-QED models. Furthermore, we subsume
$\alpha={\mathbf{k},s}$ and apply a canonical transformation that exchanges
the role of the photon coordinates and momenta

$\displaystyle-\mathbf{i}\partial_{q_{\alpha}}$
$\displaystyle\rightarrow-\omega_{\alpha}p_{\alpha}$ $\displaystyle
q_{\alpha}$
$\displaystyle\rightarrow-\mathbf{i}\omega_{\alpha}^{-1}\partial_{p_{\alpha}}.$

The non-relativistic QED Hamiltonian in the long-wavelength limit reads then

$\displaystyle\hat{H}_{LW}=\hat{H}_{n}+\hat{H}_{e}+\hat{H}_{ne}+\hat{H}_{p}+\hat{H}_{ep}+\hat{H}_{np},$
(1.35)

where we defined the nuclear Hamiltonian

$\displaystyle\hat{H}_{n}=$
$\displaystyle\sum_{i}^{N_{n}}\frac{1}{2M_{i}}\nabla_{\mathbf{R}_{i}}^{2}+\frac{1}{2}\sum_{i\neq
j}^{N_{n}}\frac{Z_{i}Z_{j}}{|\mathbf{R}_{i}-\mathbf{R}_{j}|},$ with the
nuclear masses $M_{i}$. The electronic Hamiltonian reads
$\displaystyle\hat{H}_{e}=$
$\displaystyle\sum_{i}^{N_{e}}\frac{1}{2}\nabla_{\mathbf{r}_{i}}^{2}+\frac{1}{2}\sum_{i\neq
j}^{N_{e}}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}$ and the interaction
Hamiltonian between both matter degrees of freedom is given by
$\displaystyle\hat{H}_{ne}=$ $\displaystyle\frac{1}{2}\sum_{i\neq
j}^{N_{n}}\frac{Z_{j}}{|\mathbf{r}_{i}-\mathbf{R}_{j}|}.$ The Hamiltonians for
the contributions of all the $M$ photon modes can be collected together by
$\displaystyle\hat{H}_{p}+\hat{H}_{ep}+\hat{H}_{np}=$
$\displaystyle\frac{1}{2}\sum_{\alpha}^{M}\left[-\partial_{p_{\alpha}}^{2}+\omega_{\alpha}^{2}\left(p_{\alpha}-\frac{\bm{\lambda}_{\alpha}}{\omega_{\alpha}}\cdot\mathbf{R}\right)^{2}\right],$

where the light matter interaction becomes explicit by the appearance of the
total dipole $\mathbf{R}$. We also introduced the _fundamental light-matter
coupling constant_ for mode $\alpha$

$\displaystyle\bm{\lambda}_{\alpha}=\sqrt{4\pi}S_{\alpha}(0)\bm{\epsilon}_{\alpha}.$
(1.36)

###### The Born-Oppenheimer approximation

The Hamiltonian (1.35) describes the electrons and nuclei of a matter systems,
which are basically _all_ degrees of freedom molecules and
solids.545454According to the standard microscopic model of matter, that
neglects the inner structure of the nuclei. It is thus a suited basis for a
first-principles description of, e.g., the phenomena of polaritonic chemistry.
Nevertheless, the accurate description of these degrees of freedom is very
challenging and in many cases, the electronic and nuclear dynamics can be
separated (Born-Oppenheimer approximation). The electrons are then described
in the electrostatic potential of the nuclei with fixed coordinates
$\mathbf{R}_{i}$ (clamped nuclei). The electronic energy as a function of the
$\mathbf{R}_{i}$ (potential energy surface) defines then an effective
potential for the nuclear dynamics. In this work, we focus on the electronic
part of the matter-description (electronic structure, see also Ch. 2).
However, since electrons and nuclei couple in exactly the same way (with their
dipole operator) to the photon modes, a large part of our discussion can be
generalized straightforwardly to the nuclear theory. For further details about
additional implications of the approximation for the coupled light-matter
systems and possible improvements, the reader is referred to Ref. [125].

To derive the electron-photon theory formally, we let $M_{i}\rightarrow\infty$
for $i=1,...,N_{n}$. Since $m_{e}\ll M_{i}$, this approximation has a wide
range of validity.555555However, there are also many situations, when the
Born-Oppenheimer approximation fails, but these cases are beyond the scope of
this work. The interested reader is referred to, e.g., the review by Worth and
Cederbaum [123]. Thus, the kinetic energy of the nuclei
$\sum_{i}\tfrac{1}{2M_{i}}\nabla^{2}_{\mathbf{R}_{i}}\rightarrow 0$ and we can
treat them as fixed classic charges with coordinates $\mathbf{R}_{i}$. Under
the Born-Oppenheimer approximation (abbreviated by BOA in the equations), the
nuclear contribution

$\displaystyle\hat{H}_{n}\overset{BOA}{\longrightarrow}$
$\displaystyle\frac{1}{2}\sum_{i\neq
j}^{N_{n}}\frac{Z_{i}Z_{j}}{|\mathbf{R}_{i}-\mathbf{R}_{j}|}$ is merely a
constant that we can neglect in many scenarios. We will consider the term only
for comparisons of different nuclear configurations, e.g., to study simple
chemical reactions. The electron-nuclear interaction has then a prescribed
local potential
$\hat{H}_{ne}\overset{BOA}{\longrightarrow}v(\mathbf{r})=\sum_{i=1}^{N_{n}}\frac{1}{|\mathbf{R}_{i}-\mathbf{r}|}$
that we include in the electronic Hamiltonian
$\displaystyle\hat{H}_{e}\overset{BOA}{\longrightarrow}$
$\displaystyle\hat{H}_{e}+v(\mathbf{r}).$

We also neglect the constant contribution of the nuclear charge to the total
dipole
$\mathbf{R}\overset{BOA}{\longrightarrow}-\sum_{i=1}^{N}\mathbf{r}_{i}$, where
we have removed the index of $N_{e}\rightarrow N$ for convenience. We note
that the Hamiltonian (1.35) is symmetric with respect to the simultaneous
inversion of $\mathbf{r}\rightarrow-\mathbf{r}$ and $p_{\alpha}\rightarrow-
p_{\alpha}$. Additionally, there is the ambiguity of the sign of the electric
charge [126]. That is why one can find different definitions of the bilinear
coupling term in the literature. Importantly both choices,
$p_{\alpha}\pm\frac{\bm{\lambda}_{\alpha}}{\omega_{\alpha}}\cdot\mathbf{R}$,
lead to the _same_ expectation values of the observables. We make use of this
freedom and redefine

$\displaystyle\mathbf{R}=+\sum_{i=1}^{N}\mathbf{r}_{i},$

which is the more common choice in the literature.

Collecting all the terms, we arrive at the coupled electron-photon Hamiltonian
in the long-wavelength limit or the _cavity-QED Hamiltonian_ that is
explicitly given in Def. 4. It will be the basis for (most of) the discussions
of the rest of this thesis.

List of defenvs 4 Cavity-QED Hamiltonian

A system consisting of $N$ nonrelativistic electrons that are coupled via
their dipole to $M$ modes of the electromagnetic field in the length gauge are
described by the _cavity-QED Hamiltonian_
$\displaystyle\hat{H}=\sum_{i}^{N}\left[\frac{1}{2}\nabla_{\mathbf{r}_{i}}^{2}+v(\mathbf{r}_{i})\right]+\frac{1}{2}\sum_{i\neq
j}^{N}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}+\frac{1}{2}\sum_{\alpha}^{M}\left[-\partial_{p_{\alpha}}^{2}+\omega_{\alpha}^{2}\left(p_{\alpha}+\frac{\bm{\lambda}_{\alpha}}{\omega_{\alpha}}\cdot\mathbf{R}\right)^{2}\right]$
(1.37) Here $\mathbf{R}=+\sum_{i=1}^{N}\mathbf{r}_{i}$ is the total dipole of
the matter systems and for each mode $\alpha$, $\omega_{\alpha}$ denotes the
frequency and
$\displaystyle\bm{\lambda}_{\alpha}=\sqrt{4\pi}S_{\alpha}\bm{\epsilon}_{\alpha}$
(1.38) is the light-matter coupling constant, that includes the polarization
vector $\epsilon_{\alpha}$ and the (effective) mode function $S_{\alpha}$. In
free space, the modes are plane waves and thus, $S_{\alpha}$ is given by the
simple expression $S_{\alpha}=\tfrac{1}{\sqrt{V}}$ (cf. Eq. (1.28)). However,
inside a cavity the mode functions may have a nontrivial form [127]. For later
convenience, we further differentiate the contributions
$\displaystyle\hat{H}=\hat{T}+\hat{V}+\hat{W}+\hat{H}_{ph}+\hat{H}_{int}+\hat{H}_{self},$
where
$\displaystyle\hat{T}=\sum_{i}^{N}\frac{1}{2}\nabla_{\mathbf{r}_{i}}^{2},\quad\hat{V}=\sum_{i}^{N}v(\mathbf{r}_{i}),\quad\hat{W}=\frac{1}{2}\sum_{i\neq
j}^{N}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|},$
$\displaystyle\hat{H}_{ph}=\frac{1}{2}\sum_{\alpha}^{M}\left[-\partial_{p_{\alpha}}^{2}+\omega_{\alpha}^{2}p_{\alpha}^{2}\right],\quad\hat{H}_{int}=\sum_{i}^{N}\sum_{\alpha}^{M}\omega_{\alpha}p_{\alpha}\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{i},\quad\hat{H}_{self}=\frac{1}{2}\sum_{i,j=1}^{N}\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{i}\right)\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{j}\right)$

###### A remark on the gauge choice

We conclude this subsection with a remark on the specific gauge choice in
$\hat{H}$. In the length (PZW) gauge, we describe the photon degrees of
freedom with the (transversal) displacement field

$\displaystyle\hat{\mathbf{D}}_{\perp}=$
$\displaystyle\sum_{\alpha=1}^{M}\frac{\omega_{\alpha}}{4\pi}\bm{\lambda}_{\alpha}p_{\alpha}$
as canonical variable with corresponding momentum
$\displaystyle\hat{\mathbf{B}}=$
$\displaystyle\sum_{\alpha=1}^{M}(-\mathbf{i}\partial_{p_{\alpha}})\bm{\lambda}_{\alpha}\times\epsilon_{\alpha}.$
We recover the electric field as $\displaystyle\hat{\mathbf{E}}=$
$\displaystyle
4\pi\left(\hat{\mathbf{D}}_{\perp}-\hat{\mathbf{P}}_{\perp}\right)$ with the
(transversal) polarization $\displaystyle\hat{\mathbf{P}}_{\perp}=$
$\displaystyle\sum_{\alpha}\frac{1}{4\pi}\bm{\lambda}_{\alpha}(\bm{\lambda}_{\alpha}\cdot\mathbf{R}).$

The original creation and annihilation operators for the electromagnetic field
modes have now the form565656Note that
$-\mathbf{i}\partial_{p_{\alpha}}^{+}=-\mathbf{i}\partial_{p_{\alpha}}$ is the
canonical momentum operator that is hermitian.

$\displaystyle\hat{a}_{\alpha}=$
$\displaystyle\frac{1}{\sqrt{2\omega_{\alpha}}}\left(-\mathbf{i}\partial_{p_{\alpha}}-\mathbf{i}\omega_{\alpha}p_{\alpha}+\mathbf{i}\bm{\lambda}_{\alpha}\cdot\mathbf{R}\right)$
$\displaystyle\hat{a}_{\alpha}^{+}=$
$\displaystyle\frac{1}{\sqrt{2\omega_{\alpha}}}\left(-\mathbf{i}\partial_{p_{\alpha}}+\mathbf{i}\omega_{\alpha}p_{\alpha}-\mathbf{i}\bm{\lambda}_{\alpha}\cdot\mathbf{R}\right)$

This has consequences for the interpretation of photon-observables. Most
importantly, the expectation value

$\displaystyle\tilde{E}_{ph}=$
$\displaystyle\left\langle\hat{H}_{ph}\right\rangle=\left\langle\tfrac{1}{2}\sum_{\alpha}^{M}\left[-\partial_{p_{\alpha}}^{2}+\omega_{\alpha}^{2}p_{\alpha}^{2}\right]\right\rangle$
(1.40) does not correspond to the energy of the electromagnetic field that is
defined as
$\tfrac{1}{8\pi}\int\mathrm{d}\mathbf{r}(\mathbf{E}^{2}+\mathbf{B}^{2})$ in
terms of the electric field, but is instead connected to the displacement-
field $\mathbf{D}$. To calculate the photon energy (in its standard
definition), we have to consider instead $\displaystyle E_{ph}=$
$\displaystyle\left\langle\hat{H}_{ph}+\hat{H}_{int}+\hat{H}_{self}\right\rangle$
(1.41) Correspondingly, the photon number of mode $\alpha$ is calculated as
$N_{ph}=\frac{1}{\omega_{\alpha}}E_{ph}-\frac{1}{2}=\sum_{\alpha}\braket{\hat{N}_{ph,\alpha}}$,
where we defined the photon number operator for mode $\alpha$
$\displaystyle\hat{N}_{ph,\alpha}=$
$\displaystyle\frac{1}{\omega_{\alpha}}\left(-\frac{1}{2}\frac{\partial^{2}}{\partial
p_{\alpha}^{2}}+\frac{\omega_{\alpha}^{2}}{2}\left(p_{\alpha}-\frac{\bm{\lambda}_{\alpha}}{\omega_{\alpha}}\cdot\mathbf{R}\right)^{2}\right)-\frac{1}{2},$
(1.42) For the sake of completeness, let us also define to “photon-number”
operator with respect to Eq. (1.40), i.e.,
$\displaystyle\tilde{N}_{ph,\alpha}=$
$\displaystyle\frac{1}{2\omega_{\alpha}}\left[-\partial_{p_{\alpha}}^{2}+\omega_{k}^{2}p_{\alpha}^{2}\right]-\frac{1}{2},$
(1.43)

that we call the _mode occupation_ to differentiate it from $N_{ph,\alpha}$.

This illustrates that in the length gauge, the electronic and photonic degrees
of freedom are (partly) mixed. Consequently, a part of the electron-photon
interaction is described only by matter quantities as, e.g., the occurrence of
the dipole-self-interaction $\hat{H}_{self}$ shows. This is a convenient
starting point to develop new first-principles methods for light-matter
interaction, since the description of matter is much better understood. For
instance, one could approximate the light-matter interaction by only taking
$\hat{H}_{self}$ into account and would end up with a matter-only theory,
which provides a good (qualitative) description of certain scenarios [128,
125]. Without the PWZ transformation, cf. Eq. (1.32), such a mixing does not
occur and we describe the electromagnetic field by the vector potential and
the electric field exactly as in the Coulomb gauge of NR-QED. Within the long-
wavelength description, this is called the _velocity gauge_. Instead of
$\hat{H}_{self}$, we would find in this gauge the _diamagnetic term_ ,
mentioned in Sec. 1.2, that is proportional to $\hat{A}^{2}$.

### Chapter 2 Matter from first principles: electronic-structure theory

Before we investigate how the coupled light-matter problem can be approached
from a first-principles perspective (which is the topic of Sec. 3), we discuss
in this chapter the “simpler” problem of matter-only systems. How to
accurately describe the microscopic details of matter is a topic so well
studied that many techniques have already become basic textbook knowledge.
Nevertheless, we discuss this topic in the following in great detail and with
a particular focus on the conceptual level. This serves on the one hand to
provide a complete picture also for the unacquainted reader. On the other
hand, we aim at highlighting for each method the essential ingredients that
makes it numerically efficient, yet accurate. We then show that this
efficiency is lost, once we consider the standard form of the light-matter
problem based on individual electrons and photons (Sec. 3). An analysis of
this drawback provides already the physical rationale why a polariton picture
would be preferable (part II).

The “keystone” [123] of the first-principles description of matter is the
Born-Oppenheimer approximation that we also have employed in the derivation of
the cavity-QED Hamiltonian in the previous section. Within the Born-
Oppenheimer approximation, the dynamics of the electrons and nuclei are
decoupled, which allows to separate the description of matter into methods
that focus on nuclear effects (such as vibrations or rotations) and so-called
_electronic-structure theory_. The latter comprises all approaches to describe
the electronic degrees of freedom of matter systems including the majority of
all known first-principles methods.111Note that there are also first-
principles methods that, e.g., focus explicitly on the description of the
coupled electron-nuclei dynamics [129] or only the nuclei [130]. A very
prominent role in electronic-structure theory is played by the ground state,
that determines fundamental properties of matter systems such as the
equilibrium geometry or the effective screening of the Coulomb interaction
inside a solid. It also plays a key-part in understanding chemical processes
like reactions that often can be described in the quasistatic picture of the
Born-Oppenheimer approximation. Many first-principles methods exclusively aim
to accurately describe electronic ground states and we put the focus on such
methods in the following.

The accurate description of the ground state of a many-electron system is an
extraordinary challenging problem that is a research topic since the
beginnings of quantum mechanics. Consequently, many approaches have been
developed for this task and it depends strongly on the specific scenario,
which one is the best choice. However, there are certain generic features that
play a role in all approaches.

The most important example for that is the complexity that arises from the
_many-body problem_ and that basically prohibits the exact description of most
many-electron systems. The only explicitly known way for such a description is
to solve the Schrödinger equation for the many-body wave function. For a
system consisting of $N$ electrons, the wave function depends on $4N$
coordinates,222Three spatial and one spin coordinate for each electron. and
thus, the configuration space, i.e., the space of all these wave functions
grows exponentially with $N$. For large $N$, this makes the concept of the
many-body wave function in practice useless. The complexity of the many-body
problem occurs in some way or the other in all electronic-structure approaches
and it is the reason, why practical methods always rely on approximations.

Another example of a generic feature of all electronic-structure methods is
the _particle-exchange symmetry_. Electrons are fermions and a fermionic wave
functions is antisymmetric with respect to the exchange of the coordinates of
any two electrons. As a consequence, fermions adhere to the Pauli principle,
i.e., two fermions cannot have the same quantum number. Accounting for the
Pauli principle is fundamental for an accurate description of the electronic
structure.

The third example is the _variational principle_ that defines the ground state
as the minimizer of the total energy. This suggests to determine the ground
state by an explicit minimization of the energy, which is the starting point
of every electronic-structure method. To carry out this minimization, it is
necessary to _characterize the state_ of the system, which in turn defines the
energy functional. In fact, there are several ways to do this characterization
and each of these defines a class of electronic-structure methods. In the
following, we present three different electronic-structure approaches that
adhere to the above principles:

  1. 1.

Hartree-Fock (HF) theory (Sec. 2.2)

  2. 2.

Density functional theory (DFT, Sec. 2.3)

  3. 3.

Reduced density matrices (RDMs) that are the basis of variational RDM theory
(Sec. 2.4.1) or one-body RDM functional theory (RDMFT, Sec. 2.4.2)

Let us now briefly explain these approaches, highlighting the just defined
generic principles.

HF theory is one of the simplest examples of the class of wave-function
methods, that directly approximate the many-body wave function. The HF wave
function is one Slater determinant (single-reference ansatz), which is the
simplest many-body wave function that is antisymmetric and already provides a
good qualitative description of many systems. In general, wave function
methods guarantee the particle-exchange symmetry by employing Slater
determinants. One can improve on the accuracy of HF by considering wave
functions built from more than one Slater determinant (multi-reference
ansatz). Important examples for this approach are _configuration interaction_
, _coupled cluster theory_ 333Note however that in the most common practical
coupled cluster methods, the wave function is not entirely reconstructable
[131]. or _multiconfigurational self-consistent field theory_ [132, part 5].
More recently, tensor-network methods such as _density matrix renormalization
group_ [133] theory have additionally enriched the spectrum of this approach.
Such _wave-function methods_ yield until today the most precise results for
many systems. However, their high precision usually comes with equivalently
high numerical costs. Systems with larger particle numbers of say $N\gtrapprox
100$ can rarely be described accurately with such methods. This is the
manifestation of the many-body problem in wave function methods.

In DFT instead, we drop the concept of the many-body wave function and
describe the state of a system with the (one-body) electron density
$\rho(\mathbf{r})$. The many-body problem manifests here more indirectly than
for wave-function methods. The reformulation of many-body quantum mechanics in
terms of $\rho$ leads to an unknown quantity, which is called the _exchange-
correlation functional_. Constructing approximations for this functional is a
very difficult task and severely limits the applicability of functional
methods to certain system classes. The most common starting point for the
development of functionals is the Kohn-Sham (KS) construction that describes
the state of the system with one Slater determinant exactly as in HF theory.
An important reason for this choice is that the Slater determinant accounts
for the fermionic antisymmetry. Usual approximate DFTs have numerical costs
that scale much better with the particle number than wave-function methods and
thus can in principle achieve accurate results for very large systems.444State
of the art DFT methods can describe systems with more than $1000$ particles.
See for instance Ref. [134], where the authors present results for more then
$8000$ electrons.

Employing RDMs as a basic variable provides yet another perspective on many-
electron systems. For instance, we can describe the expectation value of the
energy of an $N$-electron system exactly in terms of the 2-body RDM (2RDM).
The configuration space of such a description does not grow with the particle
number and thus a variational minimization in terms of the 2RDM is in
principle feasible. However, the many-body problem also arises in this
description in the form of conditions that determine the exact configuration
space, i.e., the set of all 2RDMs. The number of these so-called
$N$-representability conditions grows exponentially with the particle number
$N$. Thus, in practical variational minimizations of the energy with respect
to the 2RDM (variational 2RDM theory) only a subset of conditions is
considered.

A special role is here taken by the 1-body RDM (1RDM) $\gamma$ that has
comparatively simple $N$-representability conditions. These reflect the
exchange symmetry of the particle species and thus provide an alternative to
employing Slater-determinants to enforce the antisymmetry of electrons.
Although $\gamma$ is not sufficient to describe the energy of a many-electron
system in a linear way, it carries all information on the system by virtue of
a generalized Hohenberg-Kohn theorem (Gilbert’s theorem). This establishes
RDMFT as an alternative to DFT that employs $\gamma$ instead of $\rho$ as the
basic variable. Usual approximate RDMFTs are numerically more expensive than
DFT, but conversely, they can account for multi-reference effects in an easier
way.

We want to remark on a fourth and last common feature of all electronic-
structure methods, which is their computational complexity. Independently of
the details of a certain method, we always exchange a linear problem on a huge
configuration space, i.e., solving the many-body Schrödinger equation, with a
nonlinear problem on a comparatively small configuration space. Importantly,
the equations of electronic-structure (and in general many-body) methods
involve usually nonlinear operators and some form of self-consistency, which
is a challenging combination. To do this in practice, it is not sufficient to
employ, e.g., a solver for partial differential equations from a standard
library. But we usually need specific algorithms that have been designed to
solve the equations of a specific theory. Thus, electronic-structure methods
are fundamentally connected to computational mathematics. That we are nowadays
able to apply a method like HF to problems that involve many hundred electrons
is a consequence, more of algorithmic improvements than of the increase of
computer power.

####  2.1 The wavefunction, symmetry and the many-body problem

We start our survey on (equilibrium) electronic-structure theory with the
proper definition of the problem that we want to solve. We aim to accurately
describe the ground state of a system, consisting of $N$ electrons that repel
each other via the Coulomb interaction exposed to some local potential
$v(\mathbf{r})$. Usually $v$ is thought as the electrostatic potential of some
classic charges at fixed positions that represent the nuclei in the Born-
Oppenheimer approximation (see Sec. 1.3.3), but we do not have to specify this
for the general theory. This setting is described by the (electronic-
structure) Hamiltonian

$\displaystyle\hat{H}=\sum_{i}^{N}\left[\frac{1}{2}\nabla_{\mathbf{r}_{i}}^{2}+v(\mathbf{r}_{i})\right]+\frac{1}{2}\sum_{i\neq
j}^{N}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|},$ (2.1)

that can be derived in several ways, e.g., by the canonical quantization of
the corresponding classical theory or as the static limit of the Pauli-Fierz
Hamiltonian (Def. 2), i.e., by turning off the interaction to the quantized
part of the Maxwell field.555Equivalently, we can derive Hamiltonian (2.1)
from the cavity-QED Hamiltonian, cf. Eq. (1.37). Eq. (2.1) defines a linear
operator666We generally indicate operator by the ’$\hat{\phantom{a}}$’ symbol.
However, because we work in the real-space picture, we will often drop the
symbol for operators that obviously only depend on the position, such as
$\hat{\mathbf{r}}\equiv\mathbf{r}$ or $\hat{v}(\hat{\mathbf{r}})\equiv
v(\mathbf{r})$. The reason is that they act multiplicative and thus are
represented by a function. that has a well-defined spectrum and is bounded
from below.777Both follows from the proof for self-adjointness of $\hat{H}$
that is valid for a very large class of local potentials $v(\mathbf{r})$
(Kato-Rellich theorem). See for example Ref. [part 6][135]. Hence, we can
define the ground state $\Psi_{0}$ that by a variational minimization of the
energy expectation value $E=\braket{\Psi}{\hat{H}\Psi}$ over all “physical”
wave functions $\Psi$, i.e.,

$\displaystyle
E_{0}=\braket{\Psi_{0}}{\hat{H}\Psi_{0}}=\inf_{\Psi}\braket{\Psi}{\hat{H}\Psi},$
(2.2)

where we will discuss the term “physical” in the following paragraphs.888Note
that for certain potentials $v$, it can be proven that the ground state exists
and thus we can exchange the $\inf$ by a $\min$ [135, part 6.2]. This includes
the important case where
$v(\mathbf{r})=\sum_{i}^{N_{n}}1/|\mathbf{r}-\mathbf{R}_{i}|$ describes the
electrostatic potential of $N_{n}$ nuclei at the positions $\mathbf{R}_{i}$.
This _variational principle_ is of central importance in electronic-structure
theory and it will play a fundamental role for all methods that we present in
the following. Sloppily we can say that the goal of (equilibrium) electronic-
structure theory is to (approximately) solve the minimization problem (2.2)
with Hamiltonian (2.1).

#####  2.1.1 The ground state of a two-electron problem

The general definition of the electronic-structure problem by Eqs. (2.1) and
(2.2) is quite abstract and for the unacquainted reader, it might be hard to
imagine how difficult solving this problem really is. To make things more
tangible, we will put Hamiltonian (2.1) aside for a moment and consider the
problem of two electrons in one spatial dimension and without spin.

###### The one-particle problem

We start with one particle in a “box,” that we define simply by constraining
the domain of all quantities to an interval $L\subset\mathbb{R}$ and zero-
boundary conditions. The energy of this system is described by the Hamiltonian

$\displaystyle\hat{H}=\hat{H}(x)=-\tfrac{1}{2}\partial_{x}^{2},$ (2.3)

that is just the kinetic energy operator
$\hat{t}(x)=-\tfrac{1}{2}\partial_{x}^{2}$ on the domain. According to the
prescription of standard many-body quantum mechanics, we can describe the
state of the system by a wave function

$\displaystyle\Psi(x)\in\mathfrak{h}^{1}[L],$ (2.4)

where we denote with $\mathfrak{h}^{1}[L]$ the Hilbert space of square-
integrable functions.999A function $f:[a,b]\rightarrow\mathbb{C}$ is square-
integrable if its ($L^{2}$-)norm is finite, i.e.,
$\int_{a}^{b}|f(x)|^{2}\mathrm{d}x<\infty$. This is a necessary condition for
the interpretation of $|\psi(x)|^{2}$ as probability amplitude.

To determine the ground state $\Psi_{0}$ of this system, we employ the
variational principle

$\displaystyle
E_{0}=\braket{\Psi_{0}}{\hat{H}\Psi_{0}}=\inf_{\Psi\in\mathcal{C}}\braket{\Psi}{\hat{H}\Psi}.$
(cf. 2.2)

This means (assuming no degeneracy) that $\Psi_{0}$ is the one (normalized)
wave function101010Note that if $\Psi_{0}$ exists and is not degenerate, it
must be an eigenstate of $\hat{H}$, because $\hat{H}$ is self-adjoint on
$\mathfrak{h}^{1}[L]$ and thus has according to the _spectral theorem_ [135,
part 1] an eigenrepresentation that is bounded from below, i.e.,
$\hat{H}=\sum_{i=0}E_{i}\ket{\tilde{\Psi}_{i}}\bra{\tilde{\Psi}_{i}}$ with
$E_{0}<E_{1}\leq E_{2}...$. If we expand
$\Psi_{0}=\sum_{i}c_{i}\tilde{\Psi}_{i}$, we see that the lowest energy is
given by $\Psi_{0}\equiv\tilde{\Psi}_{0}$. If there is degeneracy, i.e.,
$E_{0}=E_{1}=...+E_{m}$ for some finite $m$, then
$\Psi_{0}=\sum_{i=1}^{m}c_{i}\tilde{\Psi}_{i}$. However, it might be that
$\Psi_{0}\notin\mathfrak{h}^{1}[L]$ and thus, the system cannot reach its
infimum. For instance, this happens in the free-space case, where
$L\rightarrow\infty$. with the lowest energy expectation value chosen from the
set $\mathcal{C}\subset\mathfrak{h}^{1}[L]$ of all _allowed_ wave functions,
which we call the _configuration space_.

We can solve this minimization problem simply by diagonalizing the kinetic
energy operator, i.e.,

$\displaystyle-\tfrac{1}{2}\nabla_{x}^{2}\psi_{i}(x)=k_{i}^{2}/2\psi_{i}(x),$
(2.5)

where the $\psi_{i}(x)=\sqrt{2/L}\sin(xk_{i})$ are the “box states” with the
momenta $k_{i}=i\pi/L$. Thus, the configuration space is
$\mathcal{C}=\text{span}(\psi_{1},\psi_{2},...)$ and the ground state
$\Psi_{0}=\psi_{1}(x)$ is just the box-state with the lowest allowed momentum
$k_{1}=\pi/L$.

###### The two-particle problem on a first glance

Let us now consider two particles in a box. Their energy is described by the
Hamiltonian

$\displaystyle\hat{H}=\hat{H}(\mathbf{r}_{1},\mathbf{r}_{2})=\sum_{i=1}^{2}-\tfrac{1}{2}\partial_{x}^{2},$
(2.6)

that is just the sum of the kinetic energy operators for each particle (with
index $i$) on the domain. To describe the state of the system, we need a wave
function that depends on two coordinates, i.e.,

$\displaystyle\Psi(x_{1},x_{2})\in\mathfrak{h}^{2}.$ (2.7)

The linear structure of the theory suggests now that we define the two-body
Hilbert space

$\displaystyle\mathfrak{h}^{2}[L]=\mathfrak{h}^{1}[L]\otimes\mathfrak{h}^{1}[L]$
(2.8)

simply as the product of the one-particle space $\mathfrak{h}^{1}[L]$. This is
equivalent to employing the wave function ansatz

$\displaystyle\Psi(x_{1},x_{2})=\sum_{i,j=1}^{\infty}c_{ij}\psi_{i}(x_{1})\psi_{j}(x_{2}),$
(2.9)

where we introduced the two-index expansion-coefficient $c_{ij}\in\mathbb{C}$.
To guarantee the normalization $|\Psi|^{2}=1$ of the wave function, the
coefficients need to satisfy the sum rule $\sum_{i}|c_{i}|^{2}=1$. Thus, if we
perform the minimization over $\mathcal{C}=\mathfrak{h}^{2}[L]$, we find the
ground state

$\displaystyle\Psi_{0}(x_{1},x_{2})=\psi_{1}(x_{1})\psi_{1}(x_{2}),$

with eigenvalue $E_{0}=2\pi/L$. Obviously, this state does not describe a
valid configuration of electrons, because it does _not adhere to the Pauli-
principle_.

###### What differentiates many from one: exchange symmetry

The reason is that we have not performed the minimization over the correct
configuration space. Characterizing this space is a crucial step to solve the
many-electron minimization problem. In fact, the physically correct
configuration space is _considerably smaller_ than $\mathfrak{h}^{2}$. The
true ground state of the “two electrons in a box” problem is the lowest
eigenstate of $\hat{H}$ that is also _antisymmetric_ under the exchange of the
coordinates, i.e.,

$\displaystyle\Psi(x_{1},x_{2})=-\Psi(x_{2},x_{1}).$ (2.10)

We call this the _exchange symmetry_ , which plays a fundamental role all
theoretical descriptions of quantum many-body systems. This is especially true
in the realm of first-principles methods, where we do not only utilize the
variational principle as a formal definition, but actually perform the
minimizations of the form (2.2) with large-scale numerical algorithms. To do
so, we need to _parametrize_ $\mathcal{C}$ in some way and for that, the
antisymmetry of electrons is one of the most important tools as we will see in
the following.

###### How to parametrize the antisymmetric space: Slater determinants

Let us therefore regard again our wave function ansatz (2.9) with the
expansion coefficients $c_{ij}$. Since the indices correspond explicitly to
the particle coordinates, we can simply transfer the condition (2.10) to the
coefficients by demanding

$\displaystyle c_{ij}=-c_{ji}.$

We can subsume this condition into the expansion and get the more common
form111111We redefine here the coefficients $c_{ij}\rightarrow c_{ij}\sqrt{2}$
slightly such that the normalization conditions / sum rules do not change.

$\displaystyle\begin{split}\Psi(x_{1},x_{2})=&\tfrac{1}{\sqrt{2}}\sum_{i,j}c_{ij}\left(\psi_{i}(x_{1})\psi_{j}(x_{2})-\psi_{j}(x_{1})\psi_{i}(x_{2})\right)\\\
=&1/\sqrt{2}\sum_{i,j}c_{ij}|\psi_{i}(x_{1})\psi_{j}(x_{2})|_{-}\end{split}$
(2.11)

of an expansion in terms of the (two-body) _Slater determinants_
$\Psi_{ij}(x_{1},x_{2})=1/\sqrt{2}\left(\psi_{i}(x_{1})\psi_{j}(x_{2})-\psi_{j}(x_{1})\psi_{i}(x_{2})\right)=1/\sqrt{2}|\psi_{i}\psi_{j}|_{-}$.
The set of all these two-body Slater determinants spans the antisymmetric two-
particle Hilbert space

$\displaystyle\begin{split}\mathfrak{h}^{2}_{A}=&\mathfrak{h}^{1}\wedge\mathfrak{h}^{1}\\\
\equiv&\mathcal{A}\mathfrak{h}^{2},\end{split}$ (2.12)

where the operator $\mathcal{A}$ turns the tensor product $\otimes$ into an
antisymmetric tensor product $\wedge$ or “antisymmetrizes” a given tensor
[136].

Thus, we can parametrize the correct configuration space by employing the
ansatz (2.11). This means that we calculate the matrix elements
$H_{ij}=\braket{\psi_{i}}{\hat{H}\phi_{j}}$ and reformulate the original
minimization problem (2.2) in terms of the coefficients:

$\displaystyle E[c_{ij}]=\braket{\Psi[c_{ij}]}{\hat{H}\Psi[c_{ij}]}.$ (2.13)

To calculate the 2-body matrix elements $H_{ij}$, let us see how $\hat{H}$
acts on one Slater determinant, built by the box states. We have

$\displaystyle\hat{H}\left(\psi_{i}(x_{1})\psi_{j}(x_{2})-\psi_{j}(x_{1})\psi_{i}(x_{2})\right)=$
$\displaystyle\left[\sum_{i=1}^{2}-\tfrac{1}{2}\partial_{x_{i}}^{2}\right]\left(\vphantom{\sum}\psi_{i}(x_{1})\psi_{j}(x_{2})-\psi_{j}(x_{1})\psi_{i}(x_{2})\right)$
$\displaystyle=$
$\displaystyle\left[-\tfrac{1}{2}\partial_{x_{1}}^{2}\psi_{i}(x_{1})\right]\psi_{j}(x_{2})-\left[-\tfrac{1}{2}\partial_{x_{1}}^{2}\psi_{j}(x_{1})\right]\psi_{i}(x_{2})$
$\displaystyle+\left[-\tfrac{1}{2}\partial_{x_{2}}^{2}\psi_{j}(x_{2})\right]\psi_{i}(x_{1})-\left[-\tfrac{1}{2}\partial_{x_{2}}^{2}\psi_{i}(x_{2})\right]\psi_{j}(x_{1})$
$\displaystyle=$ $\displaystyle
2(k_{i}^{2}/2+k_{j}^{2}/2)\left(\vphantom{\sum}\psi_{i}(x_{1})\psi_{j}(x_{2})-\psi_{j}(x_{1})\psi_{i}(x_{2})\right).$

We see that Slater determinants, constructed from the eigenstates of
$\hat{t}(x)$, are the eigenfunctions of the many-body Hamiltonian
$\hat{H}(x_{1},x_{2})=\hat{t}(x_{1})+\hat{t}(x_{2})$. Thus, we know that the
ground state is given by $c_{ij}=1$ for $i=0,j=1$ (or vice versa) and
$c_{kl}=0$ for all other combinations of indices: the ground state of
$\hat{H}$ is a single Slater determinant, constructed from the lowest two
eigenstates of $\hat{t}(x)$.

###### Including a local potential

Let us now generalize the description and consider a local potential
$\hat{v}(x)$. Since the potential now should describe the geometry of the
problem, we can formally assume $L\rightarrow\infty$. The many-body
Hamiltonian reads then

$\displaystyle\hat{H}=\sum_{i=1}^{2}\hat{t}(x_{i})+\hat{v}(x_{i}).$ (2.14)

We can straightforwardly generalize our solution strategy from before, if
instead of Eq. (2.5), we consider the eigenvalue problem

$\displaystyle[\hat{t}(x)+\hat{v}(x)]\psi_{i}(x)=e_{i}\psi_{i}(x).$ (2.15)

The physical ground state is simply given as the Slater determinant
$\Psi_{0}(x_{1},x_{2})=1/\sqrt{2}|\psi_{0}\psi_{1}|_{-}$ with ground state
energy $E_{0}=e_{0}+e_{1}$.

###### Solving a minimization problem in practice: the basis set

For almost all choices of $v$, Eq. (2.15) cannot be solved analytically and
requires a numerical solver. This means that we need to parametrize the
single-particle states (orbitals) $\psi_{i}$ in a _known basis set_ with say
$B$ elements, i.e., $\psi_{i}=\sum_{j=1}^{B}b_{i}^{j}\xi_{j}$. An obvious
choice would be a (sufficiently large) discretized interval in real space,
such that $\xi_{j}(x)=\theta_{j}(x)$, where $j=1,\dotsc,B$ denote the $B$ grid
points and $\theta_{j}(x)$ defines the discretization.121212A possible choice
would be $\theta_{x_{j}}=d_{x}^{-1/2}\Theta(d_{x}/2-|x_{i}-x|)$, where $d_{x}$
is the spacing and $\Theta$ denotes the Heaviside step-function. The
differential operators would then be approximated by, e.g., _finite-
differences_ (see Ch. 7). This will introduce an error, which however can be
controlled, i.e., with increasing $B$, the approximate solutions will converge
to the exact results. We call such a procedure convergence with respect to the
basis set.

Of course, we can also choose a different basis and in fact, discretizing the
real-space is often not the best choice. For example, we might need very large
$B$ to converge with respect to the basis set. Another option is to employ the
eigenfunctions of the kinetic energy operator, which for $L\rightarrow\infty$
are _plane waves_ [137, part 2]. Employing this basis is equivalent to
describing the whole problem in Fourier space. The obvious advantage of such a
choice is the analytically known form of the basis for all $x\in\mathbb{R}$
and thus, the description is not constrained by a box of finite volume.
However, the accurate description of small distances will be inefficient in
this basis. In fact, plane waves are a standard tool to describe the bulk of
many condensed-matter systems, which are modelled as infinitely extended
periodic systems [137]. These two examples demonstrate how strongly the basis
choice depends on our _knowledge_ of the problem.

Plane waves and discretized real-space bases are in this sense the most
generic bases. To describe for instance finite-systems, usually neither of the
two is used, but so-called _atomic orbitals_. These are especially designed
bases with the aim to reduce the necessary number of orbitals to a minimum.
The basic idea is to make use of the analytically known eigenstates of the
Hydrogen atom and the most straightforward approach is thus to use the lowest
say $B_{a}$ of these eigenfunctions centered at the positions of all the
$N_{n}$ nuclei of the system. The basis set consists then of $B=B_{a}N_{n}$ of
such atomic orbitals, which is in most cases considerably smaller then any
basis set obtained from a parametrization of the discretized real space.
Another very common approach is to employ Gaussian functions instead of the
Hydrogen orbitals because of their convenient mathematical properties. These
are the two most employed strategies to obtain quantum chemical basis-
sets.131313The interested reader is referred to, e.g., the corresponding
chapters in the textbook of Helgaker et al. [132].

However, we want to stress that also real-space bases are used in
practice141414Also in this work, see Ch. 6. because they allows for a
description of systems with less foreknowledge and a considerably easier
visualization [138].

###### The interacting problem

To finalize our survey, there is one last missing piece, which is the
inclusion of the interaction. The repulsion between two electrons in 3d is
described by the Coulomb operator that reads (employing the coordinates
$\mathbf{r}_{1},\mathbf{r}_{2}$)

$\displaystyle\hat{w}(\mathbf{r}_{1},\mathbf{r}_{2})=\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}.$
(2.16)

This is a two-body operator, because it depends on two coordinates, but in
contrast to, e.g., Eq. (2.14), it _cannot_ be expressed as a sum of two one-
body operators. We call this the _order_ of an operator: $\hat{w}$ has the
order two, but for instance $\sum_{i=1}^{2}[\hat{t}(x_{i})+\hat{v}(x_{i})]$
has the order one. It is a fundamental property of interaction operators that
their order is larger than one. Importantly, this spoils our solution strategy
from above.

To demonstrate this, let us consider a generic interaction operator in 1d
$w(x_{1},x_{2})$151515Note that the one-dimensional version of the Coulomb
operator $w(x_{1},x_{2})=\frac{1}{|x_{1}-x_{2}|}$ has very different
properties than its 3d version. Thus, in 1d-model calculations the soft-
Coulomb potential
$w(x_{1},x_{2})=\frac{1}{\sqrt{(x_{1}-x_{2})^{2}+\epsilon^{2}}}$ with
$\epsilon>0$ is usually considered (see also Sec. 2.5). that we add to
Hamiltonian (2.14), i.e.,

$\displaystyle\hat{H}=\sum_{i=1}^{2}\hat{t}(x_{i})+\hat{v}(x_{i})+\hat{w}(x_{1},x_{2}).$
(2.17)

We start with the orbital set that we need to construct the Slater
determinants. A reasonable choice could be the eigenstates of
$\hat{t}+\hat{v}$ that we know to calculate numerically. However, Slater
determinants of these orbitals cannot be eigenstates of the interacting
problem and thus, we have to consider the general ansatz

$\displaystyle\Psi(x_{1},x_{2})=$
$\displaystyle\tfrac{1}{\sqrt{2}}\sum_{i,j}c_{ij}\left(\psi_{i}(x_{1})\psi_{j}(x_{2})-\psi_{j}(x_{1})\psi_{i}(x_{2})\right)$
(cf. 2.11) $\displaystyle=$ $\displaystyle\sum_{i,j}c_{ij}\Psi_{ij}.$

If we apply the full Hamiltonian to one of these Slater determinants
$\Psi_{ij}$,

$\displaystyle\hat{H}\Psi_{ij}=$
$\displaystyle\left[\sum_{i=1}^{2}\hat{t}(\mathbf{r}_{i})+\hat{v}(\mathbf{r}_{i})+\hat{w}(\mathbf{r}_{1},\mathbf{r}_{2})\right]\Psi_{ij}$
$\displaystyle=$
$\displaystyle\underbrace{\left[\sum_{i=1}^{2}\hat{t}(\mathbf{r}_{i})+\hat{v}(\mathbf{r}_{i})\right]\Psi_{ij}}_{2(e_{1}+e_{2})\Psi_{ij}}+\hat{w}(\mathbf{r}_{1},\mathbf{r}_{2})\tfrac{1}{\sqrt{2}}\left(\psi_{i}(x_{1})\psi_{j}(x_{2})-\psi_{j}(x_{1})\psi_{i}(x_{2})\right),$

we see that the interaction term _couples_ the two orbitals of $\Psi_{ij}$ for
every $i,j$. This means that (without further information on the systems like
symmetries etc.) we cannot simplify the general ansatz (2.11) for an
interacting problem. Even if we chose a highly optimized orbital basis, we
have to perform the minimization

$\displaystyle E[c_{ij}]=\braket{\Psi[c_{ij}]}{\hat{H}\Psi[c_{ij}]}.$ (cf.
2.13)

with respect to in principle _all_ entries of the coefficient matrix $c_{ij}$.

#####  2.1.2 Solving a generic many-electron problem and the exponential wall

Let us now transfer the concepts and tools from the two-electron example to
the general case of $N$ electrons in three spatial dimensions, including spin.
Thus, we consider spin-spatial coordinates $\mathbf{x}=(\mathbf{r},\sigma)$,
where $\mathbf{r}\in\mathbb{R}^{3}$ describes a position and
$\sigma\in\\{\uparrow,\downarrow\\}$ is the spin-quantum number.

###### The antisymmetric $N$-electron space

We start with the definition of the antisymmetric $N$-particle space, i.e.,
the generalization of Eq. (2.12). Therefore, we consider again the single-
particle Hilbert space $\mathfrak{h}^{1}$ that is spanned by some spin-spatial
orbital basis $\\{\bar{\psi}_{i}(\mathbf{x})\\}$ (where we assume basis states
to be orthonormalized, i.e.,
$\int\mathrm{d}\mathbf{r}\bar{\psi}_{i}^{*}(\mathbf{r})\bar{\psi}_{j}(\mathbf{r})=\delta_{ij}$),
such that every single-particle state can be expressed as

$\displaystyle\Psi(\mathbf{x})=\sum_{i=1}^{\infty}c_{i}\bar{\psi}_{i}(\mathbf{x}),$
(2.18)

where the $c_{i}\in\mathbb{C}$ are the expansion coefficients (with
corresponding sum rule $\sum_{i,j}|c_{ij}|^{2}=1$ for normalization). From
this starting point, we can now successively build Hilbert spaces
$\mathfrak{h}^{2},\mathfrak{h}^{3}$ for more particles as antisymmetric
product spaces, i.e.,

$\displaystyle\mathfrak{h}_{A}^{N}=\bigwedge\limits_{l=1}^{N}\mathfrak{h}^{1}.$
(2.19)

The basis states of $\mathfrak{h}^{N}$ are the $N$-body Slater determinants

$\displaystyle\Psi_{I}(\mathbf{x}_{1},\dotsc,\mathbf{x}_{N})=$
$\displaystyle\frac{1}{\sqrt{N!}}\sum_{\pi_{j}\in
P_{N}}(-1)^{j}{\bar{\psi}}_{\pi_{j}(I_{1})}(\mathbf{x}_{1})\dotsb{\bar{\psi}}_{\pi_{j}(I_{N})}(\mathbf{x}_{N})$
$\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt{N!}}|\bar{\psi}_{I_{1}},\dotsc,\bar{\psi}_{I_{N}}|_{-},$
(2.20)

where $P_{N}$ denotes the permutation group on $N$ elements and the index $j$
is chosen such that it is even (odd) for an even (odd) permutation $\pi_{j}\in
S_{N}$. Thus, every wave function
$\Psi(\mathbf{x}_{1},\dotsc,\mathbf{x}_{N})\in\mathfrak{h}_{A}^{N}$ can be
parametrized as

$\displaystyle\begin{split}\Psi(\mathbf{x}_{1},\dotsc,\mathbf{x}_{N})=&\sum_{I}c^{I}\Psi_{I}(\mathbf{x}_{1},\dotsc,\mathbf{x}_{N})\\\
=&\Psi[c^{I}],\end{split}$ (2.21)

where $c^{I}=c^{I}_{I_{1},\dotsc,I_{N}}\in\mathbb{C}$ is the coefficient of
the $I$-th Slater determinant including the orbitals with indices
$I_{1},\dots,I_{N}$. The according sum rule is $\sum_{I}|c^{I}|=1$. Note that
the introduction of Slater determinants is one of the few known ways to
guarantee the antisymmetry of a wave function and it is thus a fundamental
tool of most electronic-structure theories.

###### The closed-shell setting

In the following, we will usually consider the special case of _closed-shell_
systems, which is very common in electronic-structure theory [132] and covers
a broad range of in practice relevant systems.161616For instance, this
excludes systems where certain magnetic effects or spin-orbit coupling play a
role [137]. The generalization to odd $N$ and spin-dependent Hamiltonians is
in principle straightforward. However, approximate methods for the spin-
restricted case are generally better understood. For that, we assume that $N$
is even and define

$\displaystyle\begin{split}\bar{\psi}_{2i}(\mathbf{r},\sigma)=&\psi_{i}(\mathbf{r})\alpha(\sigma)\\\
\bar{\psi}_{2i+1}(\mathbf{r},\sigma)=&\psi_{i}(\mathbf{r})\beta(\sigma),\end{split}$
(2.22)

where $\\{\alpha,\beta\\}\equiv\\{\uparrow,\downarrow\\}$ denote the two spin
functions.171717We choose this nomenclature, because it is very common in the
quantum chemistry literature.

This means that each spatial orbital is occupied twice and the system is spin-
saturated. We can thus remove the spin-coordinate from our description (by a
trivial integration).

###### The ground state of a many-electron system

With the ansatz (2.21), we can reformulate the original minimization problem
(2.2) in the following way: First, we choose a specific _basis_
$\\{\psi_{j}(x)\\}$, expand $\Psi=\Psi[c^{I}]$ in Slater determinants of this
basis according to Eq. (2.21) and calculate the general energy expression as a
functional of the expansion coefficients $c^{I}$, i.e.,

$\displaystyle E[c^{I}]=\braket{\Psi[c^{I}]}{\hat{H}\Psi[c^{I}]}.$ (2.23)

To find the ground state, we then minimize $E[c^{I}]$ with respect to the
$c^{I}$. The solution $c^{I}_{0}$ defines $\Psi_{0}=\Psi[c^{I}_{0}]$.

###### The non-interacting case

Importantly, the just described strategy strongly simplifies for every systems
that is described by operators of order one, i.e., operators that can be
decomposed as

$\displaystyle\hat{O}(\mathbf{r}_{1},...,\mathbf{r}_{N})=\sum_{i=1}^{N}\hat{o}(\mathbf{r}_{i}).$
(2.24)

An important example for this are non-interacting systems with the Hamiltonian

$\displaystyle\hat{H}=\sum_{i}^{N}\left[\frac{1}{2}\nabla_{\mathbf{r}_{i}}^{2}+v(\mathbf{r}_{i})\right].$
(2.25)

To find the ground state of $\hat{H}$, we simply have to calculate the $N/2$
(assuming that $N$ is even) lowest eigenstates of
$\frac{1}{2}\nabla_{\mathbf{r}}^{2}+v(\mathbf{r})$ and construct a Slater
determinant out of those, using Eq. (2.22). This construction of a many-body
state by successively “filling up” one-body states is often called _aufbau
principle_ (see Sec. 2.1.3).

Thus, the non-interacting $N$-electron problem (and any other problem of order
one) _reduces_ effectively to $N$ one-electron or orbital problems. With
state-of-the-art algorithms and higher-performance clusters, such (and
nonlinear versions of such) eigenvalue problems for orbitals can be very
efficiently solved even for spatially very extended systems (that in turn
require large orbital bases). We call a wave function consisting only of one
Slater determinant a _single-reference_ wave function.

###### Interaction and the exponential wall

As in the two-particle case, any form of interaction spoils this simple
solution strategy. For instance, the electronic-structure Hamiltonian (Eq.
(2.1)) includes the Coulomb interaction operator $\frac{1}{2}\sum_{i\neq
j}^{N}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}$ that has order two. This
operator in principle couples all Slater determinants of the general (or
_multi-reference_) ansatz

$\displaystyle\Psi(\mathbf{x}_{1},\dotsc,\mathbf{x}_{N})=\sum_{I}c^{I}\Psi_{N}^{I}(\mathbf{x}_{1},\dotsc,\mathbf{x}_{N}),$
(cf. 2.21)

with the $N$-dimensional coefficient tensor

$\displaystyle c^{I}=c^{I}_{j_{1},\dotsc,j_{N}}.$ (2.26)

Thus, we cannot simplify the expansion as we have done for the non-interacting
case and need to consider the full coefficient tensor. If we consider $M$
basis states, this tensor has a dimension of

$\displaystyle D=M^{N},$ (2.27)

that grows _exponentially_ with the particle number. If we wanted to compute
the wave function for say a Benzene molecule ($C_{6}H_{6}$) that has
$N=6^{*}6+6=42$ electrons with a small basis set of say $M=10$ orbitals, the
coefficient tensor would have $D=10^{42}$ entries, which could hardly be
determined by any imaginable computation system. For a system with $N=80$
electrons, we approach the approximated number of particles of the
universe.181818Source: https://www.universetoday.com/36302/atoms-in-the-
universe/amp/, accessed 15.07.2020. Clearly, we can exclude some of this
entries by smart basis choices and symmetry considerations. State-of-the-art
methods have been employed to calculate the (real-space) many-body wave
function (for certain systems) of about 50 electrons [12].191919In certain
(usually effectively one-dimensional) model descriptions, the many-body wave
function has been calculated even for larger $N$. See for example Motta et al.
[139]. However, the _scaling_ of Eq. (2.27) remains and represents an
“exponential wall” as Nobel laureate Walter Kohn called it [13].

The problem of the exponential wall is so severe that the concept of the many-
body wave function for larger systems has basically no value beyond its use
for certain theoretical considerations. Walter Kohn went so far as to suggest
that

> “in general the many-electron wavefunction
> $\Psi(\mathbf{r}_{1},...,\mathbf{r}_{N})$ for a system of $N$ electrons is
> not a legitimate concept, when $N>N_{0}$, where $N_{0}\approx 1000$.”
> (Kohn,1999 [13])

Regarding the above-mentioned current state-of-the-art, Kohn’s guess of
$N_{0}$ seems even understated. In practice, this means that most of our
knowledge on electronic-structure relies on alternative descriptions that
carry less information than the full many-body wave function.

#####  2.1.3 Excursion: exchange symmetry in quantum mechanics

Since exchange symmetry will play a central role in the course of this thesis,
we want to make a little detour recapitulating the concept of exchange
symmetry. This more historically excursion supplements the chapter and is not
necessary to follow the course of this text. The fast reader can directly
continue with Sec. 2.2.

###### The origin of quantum indistinguishability

To the best of our knowledge, all electrons have identical physical
properties, i.e., the same mass, charge and spin and thus, a valid definition
of a many-body Hamiltonian like Eq. (2.6) or Eq. (2.1) must be _symmetric_
with respect to the exchange of particles coordinates. From the symmetry of
the Hamiltonian we can deduce the exchange symmetry or “special circumstances”
of the wave function as Wolfgang Pauli called it in one of his famous articles
on quantum mechanics:

> Wenn wir es mit mehreren gleichartigen Teilchen zu tun haben, treten
> besondere Verhältnisse ein, die daher rühren, daß der Hamiltonoperator stets
> invariant ist bei irgendwelchen Vertauschungen der Teilchen.
> _If we deal with many similar [indistinguishable] particles, special
> circumstances occur that arise from the Hamilton operator being invariant
> under any permutations of particles._ (Pauli, 1933 [140])

It is clear that also a classical system consisting of many identical
particles would be described by a symmetric Hamiltonian. However, such
classical particles (imagine billiard balls) are distinguishable, even if they
have identical physical properties. The reason is that we can in principle
locate their position at any instance of time and collect them in a
trajectory. This allows to define a label for every trajectory, which in turn
distinguishes the balls. The situation changes, if we employ, e.g., a
_statistical_ description. This means that instead of describing explicitly
the trajectory of each billiard ball, we employ a probability distribution
that describes the probability to find _a_ billiard ball at a certain position
at a given time instance. In such a description, we obviously lose the
distinguishability.

The most common formulations202020For a discussion on possible alternative
approaches, see for example Ref. [141]. of quantum mechanics are also
probabilistic and thus they are (at least in practice) similar to the latter
case.212121However, ontologically there are huge differences between classical
mechanics and quantum mechanics, were the details depend on the specific
interpretation of quantum mechanics. We refer the interested reader to the
well-documented Wikipedia article on the topic and the reference therein:
https://en.wikipedia.org/wiki/Interpretations_of_quantum_mechanics.A typical
argument why quantum particles must be indistinguishable extends on the
trajectory example above: the knowledge of a trajectory in time and space
contradicts Heisenberg’s uncertainty principle [142].222222From this argument
follows that quantum particles that are very far away from each other (such
that they can be described by separate wave functions that do not overlap) are
indeed distinguishable. Dicke used this argument, when he introduced his model
of $N$ two-level system coupled to one photon mode [72]. A less “hand-wavy”
argument by means of the counting of states in a fictive experiment was given
by Dirac already in the very beginnings of the development of quantum
mechanics [143]. For a very detailed discussion of this argument, the reader
is referred to the book by Stefanucci and Van Leeuwen [144, part 1].

###### The consequence of indistinguishability: particle statistics

To the best of our knowledge, the indistinguishability of quantum particles is
unquestioned until nowadays and usually regarded as a fact [142, 137, 144].
This means in practice that although we formally need to introduce (labelled)
coordinates to describe $N$ electrons, we have to make sure that the
predictions of the theory cannot be used to differentiate the electrons. Since
only _expectation values_ of observable operators correspond to possible
measurements in experiments, we must guarantee indistinguishability on this
level. For instance, to obtain the expectation value of the energy of a two-
particle system, we need to calculate the integral

$\displaystyle
E=\int\mathrm{d}\mathbf{x}_{1}\mathrm{d}\mathbf{x}_{2}\Psi^{*}(\mathbf{x}_{1},\mathbf{x}_{2})\hat{H}(\mathbf{x}_{1},\mathbf{x}_{2})\Psi(\mathbf{x}_{1},\mathbf{x}_{2}).$
(2.28)

Since the wave function occurs quadratically in such expressions, we have more
freedom to do such a symmetrization on the level of $\Psi$ than on the
operator level. Dirac [143] concluded already in 1926 that as a consequence,
many-body wave functions need to be either symmetric or
antisymmetric232323Note that there are effective Hamiltonians that have even
more complicated symmetries, e.g., the effective Hamiltonian of certain
topologically nontrivial systems. See Ref. [137, part 9] for a good overview
about such scenarios. under a pair-wise permutation
$\mathcal{P}_{\mathbf{x},\mathbf{y}}$ of two coordinates
$\mathbf{x},\mathbf{y}$:

$\displaystyle\mathcal{P}_{\mathbf{x}_{1},\mathbf{x}_{2}}\Psi(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle=\Psi(\mathbf{x}_{2},\mathbf{x}_{1})\overset{!}{=}\Psi(\mathbf{x}_{1},\mathbf{x}_{2})\quad$
$\displaystyle\longleftrightarrow\quad\text{symmetry}$ (2.29a)
$\displaystyle\mathcal{P}_{\mathbf{x}_{1},\mathbf{x}_{2}}\Psi(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle=\Psi(\mathbf{x}_{2},\mathbf{x}_{1})\overset{!}{=}-\Psi(\mathbf{x}_{1},\mathbf{x}_{2})\quad$
$\displaystyle\longleftrightarrow\quad\text{antisymmetry}.$ (2.29b)

These symmetry properties _define_ the two fundamental particle classes, that
we consider until today:242424Note that particles with different symmetry
properties can be defined (See also the previous footnote). However, such
definitions require some kind of effective theory and thus the according
particles are usually not considered as fundamental. This holds also for the
polaritons that we will introduce in part II. They are electron-photon hybrid
particles. we call particles that are described by a symmetric or
antisymmetric many-particle wave function _bosons_ or _fermions_ ,
respectively. In Ref. [143], Dirac further connected the fundamental exchange
symmetry of the particle to its statistical properties: symmetric particles
follow Bose-Einstein statistics [145, 146] and antisymmetric particles follow
Fermi-Dirac [147] statistics. The different particle statistics are
undoubtedly one of the most important elements of quantum theory and alone,
i.e., _without_ taking into account many other elements of the theory, they
have an impressive explanatory power. Einstein showed for example already in
1925 [146], i.e., _before_ Heisenberg and Schrödinger formulated the first
well-defined theories of quantum mechanics [148, 149], that Bose-Einstein
gases exhibit a phase state, the _condensate_ , which is characterised by
extreme coherence between all the gas particles, letting them act as one large
object.252525Already Einstein predicted this peculiar phase state in the cited
paper but it took until 1995 for its first experimental realization [150].
Since then, many more Bose-Einstein condensates have been realized
experimentally and one can find an considerable amount of literature about the
topic. A good introduction to the topic can be found in the textbook by
Altland and Simons [137, part 6.3].

###### Pauli’s principle

Regarding the explanatory potential of fermion statistics, many phenomena
could be named, especially in the realm of electronic-structure theory, where
electrons are described as fermions. Most importantly, one can show that
fermions obey _Pauli’s exclusion principle_ , which is a crucial element of
the modern model of atoms. Pauli defined his principle in the following way:

> Es kann niemals zwei oder mehrere äquivalente Elektronen im Atom geben, für
> welche […] die Werte aller Quantenzahlen […] übereinstimmen. Ist ein
> Elektron im Atom vorhanden, für das diese Quantenzahlen bestimmte Werte
> haben, so ist dieser Zustand “besetzt”.
> There are never two or more equivalent electrons in the atom which have the
> same quantum number. If there is an electron in the atom with certain
> quantum numbers, this state is “occupied.” (Pauli, 1925 [151])

A remarkable deduction from Pauli’s principle is for example the explanation
of the _stability of matter_. Lieb [152] showed under very weak assumptions
employing the atom model of electrons and static nuclei that Pauli’s principle
is necessary and sufficient to guarantee that _any_ matter system built by
such constituents is stable, i.e., it has a well-defined ground state.

###### A historical example of the importance of Pauli’s principle

To illustrate the importance of Pauli’s principle, let us make a very brief
detour to the 1930ies that is based on the first chapters of the book by
Gavroglu and Simões [153]. At that time, chemistry and physics were two
(almost completely) separate research fields. But with the formulation of
quantum mechanics and the corresponding model of the atom, physicists paved
the way for establishing _quantum chemistry_. In this since then ever growing
research field, researchers try to understand physical models and experiments
of chemical systems by connecting quantum theory with the theoretical
knowledge of chemistry. However, the development of quantum chemistry was
initially very slow and it took a lot of time and especially modern computers
until the real success story of the field started. Physicists like Dirac
started to promote the point of view that chemistry was merely an application
of physical laws, but they had to realize very fast that this application is
much more difficult than anybody could have ever imagined. Walter Kohn
summarized this in his Nobel lecture in the following way:

> “There is an oral tradition that […] Dirac declared that chemistry had come
> to an end - its content was entirely contained in that powerful [i.e., the
> Schrödinger] equation. Too bad, he is said to have added, that in almost all
> cases, this equation was far too complex to allow solution.”  (Kohn,1999
> [13])

On the other hand, Chemistry always relied (and is still relying) on heuristic
rules and tables, like Mendeleev’s perodic table [154], Hund’s rule [155] or
Mulliken’s “correlation diagram” [156]. And most chemists doubted that the
mathematically very challenging quantum theory could be ever useful for
understanding such kind of rules. The only exception is Pauli’s principle,
which is somewhat special in comparison to other concepts like the wave
function or the Hamilton operator. In fact, Pauli’s principle is a rule: Only
by utilizing the principle and a simple version of the physicist’s atom model,
chemists could generalize and combine several of their own rules and concepts.
Thus, from all the complex apparatus of the newly developed quantum theory, it
was the conceptually most simple part, Pauli’s principle, that played a key
role in establishing the fruitful connection between physics and chemistry. It
is needless to say that this connection played a very important role in the
further development of quantum mechanics itself. Nowadays, we know that
solving Schrödinger’s equation for many-particle systems like larger molecules
is de facto impossible and a first-principles description of such systems
requires methods that have been developed by mathematicians, physicists and
chemists, often in close collaboration. We will see in the following that
Pauli’s principle plays an important role in all these methods.

###### Bohr’s aufbau principle

Another consequence of fermion statistics that plays an important role in
electronic-structure theory is _Bohr’s aufbau principle_ (building-up
principle), which Pauli defined in the following way:

> Dieses von Bohr aufgestellte Prinzip besagt, daß bei Anlagerang eines
> weiteren Elektrons an ein Atom die Quantenzahlen der schon gebundenen
> Elektronen dieselben Werte behalten, die ihnen im zugehörigen stationären
> Zustand des freien Atomrestes zukommen.
> This by Bohr established principle states that when a further electron is
> absorbed, the quantum numbers of an atom’s bound electron keep the values,
> that they would have in the corresponding stationary state of the free atom
> rest. (Pauli, 1925 [151])

The strength of the aufbau principle is especially visible in effective
single-particle methods such as HF theory (see Sec. 2.2) to describe many-body
systems. In such methods, we solve (nonlinear) eigenvalue equations for
single-particle states (orbitals). The aufbau principle then tells us how to
“occupy” these orbitals with electrons to obtain the many-body state: if our
systems consists of $N$ electrons, we will occupy the $N$ energetically lowest
orbitals. In this sense, the aufbau principle combines the exclusion principle
(we only populate every state by one electron) with the variational principle
(we try to obtain the lowest energy). However, other aufbau principles
inspired by Bohr’s principle have been developed, which in one or the other
way fill up some levels.262626See for instance Refs. [157, 158, 159]. This
shows the strength of the very concept of an aufbau principle in electronic-
structure theory.

####  2.2 Hartree-Fock theory

In the last section, we have discussed the problems connected to the
$N$-electron wave function that for a general interacting system is a
superposition of arbitrary Slater determinants, i.e.,

$\displaystyle\Psi(\mathbf{x}_{1},\dotsc,\mathbf{x}_{N})=\sum_{I}c^{I}\Psi_{N}^{I}(\mathbf{x}_{1},\dotsc,\mathbf{x}_{N}).$
(cf. 2.21)

There is however the special setting of non-interacting electrons, for which
the $N$-body problem basically reduces to a one-body problem. The
corresponding many-body state is given by one Slater determinant. In HF
theory, we try to find the “best fit” between the non-interacting and
interacting problem by simply truncating the expansion (2.21) after the first
element, i.e., we consider a single-reference wave function ansatz

$\displaystyle\Psi(\mathbf{x}_{1},\dotsc,\mathbf{x}_{N})\approx\Psi^{HF}(\mathbf{x}_{1},\dotsc,\mathbf{x}_{N})=\tfrac{1}{\sqrt{N!}}|\bar{\psi}_{1}\dotsc,\bar{\psi}_{N}|_{-}.$
(2.30)

to describe the interacting problem.

###### The HF minimization problem

For our derivations, we assume _spin-restriction_ , i.e., we consider systems
with even particle number and define

$\displaystyle\begin{split}\bar{\psi}_{2i}(\mathbf{r},\sigma)=&\psi_{i}(\mathbf{r})\alpha(\sigma)\\\
\bar{\psi}_{2i+1}(\mathbf{r},\sigma)=&\psi_{i}(\mathbf{r})\beta(\sigma),\end{split}$
(cf. 2.22)

where $\\{\alpha,\beta\\}$ denote the two spin functions (see the discussion
around Eq. (2.22)). This _restricted HF_ is the most common version of HF
theory. We saw already that wave functions like $\Psi^{HF}$ are the solutions
of non-interacting problems, i.e., $\Psi^{HF}$ is the ground state of any
Hamiltonian of the form

$\displaystyle\hat{H}_{ni}=\sum_{i=1}^{N}\hat{t}(\mathbf{r}_{i})+\hat{v}(\mathbf{r}_{i}),$
(2.31)

if we identify $\psi_{1},..,\psi_{N/2}$ according to the aufbau principle with
the lowest $N/2$ eigenstates of $\hat{t}+\hat{v}$. One can thus say that in
single-reference methods, we approximate the interacting with a non-
interacting wave function.272727This statement is independent of the assumed
spin-restriction of our example. However, there are many possibilities for
such an approximation and we have to define a qualifier that determines one
out of these. In wave-function methods such as HF theory, this qualifier is
always the energy (we will employ a different qualifier in Sec. 2.3). We thus
can define the HF ground state by the variational principle, i.e.,

$\displaystyle
E_{0}^{HF}=\braket{\Psi^{HF}_{0}}{\hat{H}\Psi^{HF}_{0}}=\inf_{\Psi^{HF}}\braket{\Psi^{HF}}{\hat{H}\Psi^{HF}},$
(2.32)

with the $N$-electron electronic-structure Hamiltonian that we defined in the
previous section

$\displaystyle\hat{H}=\sum_{i=1}^{N}\hat{t}(\mathbf{r}_{i})+\hat{v}(\mathbf{r}_{i})+\tfrac{1}{2}\sum_{i,j=1}^{N}\hat{w}(\mathbf{r}_{i},\mathbf{r}_{j}).$
(cf. 2.1)

###### The features of HF theory: the upper energetic bound and the definition
of correlation

It is clear that the ansatz $\Psi\approx\Psi^{HF}$ constrains the full
configuration space that we have to consider to find the exact ground state of
$\hat{H}$. Since due the variational principle, $\Psi$ must be the wave
function with the lowest possible energy expectation value, we can deduce that
the HF energy

$\displaystyle E_{0}^{HF}\geq E_{0}$ (2.33)

is an upper bound of the exact many-body energy $E_{0}$. We call methods with
this property _variational_. The difference between exact and HF energy, the
so-called correlation energy

$\displaystyle E_{c}=E_{0}-E_{0}^{HF}\leq 0,$ (2.34)

thus must be zero or negative. Consequently, HF defines a kind of “baseline”
for all electronic-structure methods that from this perspective “only” aim to
accurately describe $E_{c}$. Because the HF wave function corresponds to a
non-interacting system that as we have seen in the last section can be reduced
to an effective one-body problem, the correlation energy is often called
_many-body energy_.

###### The HF energy: the “quantum” contribution and exchange symmetry

Let us now calculate this energy expression explicitly. We have

$\displaystyle E_{HF}=$ $\displaystyle\braket{\Psi^{HF}}{\hat{H}\Psi^{HF}}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\braket{\Psi^{HF}}{\left[\hat{t}(\mathbf{r}_{i})+\hat{v}(\mathbf{r}_{i})\right]\Psi^{HF}}+\tfrac{1}{2}\sum_{i,j=1}^{N}\braket{\Psi^{HF}}{\hat{w}(\mathbf{r}_{i},\mathbf{r}_{j})\Psi^{HF}},$

where the first, i.e., the one-body part (after some algebra, see Ref. [132,
part 5.4]) reads

$\displaystyle E^{(1)}[\\{\psi_{i}\\}]=$ $\displaystyle
2\sum_{i=1}^{N/2}\int\mathrm{d}\mathbf{r}\psi_{i}^{*}(\mathbf{r})\left[\hat{t}(\mathbf{r})+\hat{v}(\mathbf{r})\right]\psi_{i}(\mathbf{r})$
$\displaystyle=$ $\displaystyle
2\sum_{i=1}^{N/2}\braket{\psi_{i}}{\hat{t}+\hat{v}}{\psi_{i}},$

where we have defined the one-body integrals
$\braket{\psi_{i}}{\hat{o}^{(1)}}{\psi_{i}}=\int\mathrm{d}\mathbf{r}\psi_{i}^{*}(\mathbf{r})\left[\hat{o}(\mathbf{r})+\hat{v}(\mathbf{r})\right]\psi_{i}(\mathbf{r})$
for a given one-body $\hat{o}^{(1)}$. Note that the factor 2 is due to the
spin-restriction. The second term involves the two-body operator
$\hat{w}(\mathbf{r},\mathbf{r}^{\prime})$ and reads [132, part 5.4]

$\displaystyle E^{(2)}[\\{\psi_{i}\\}]=$
$\displaystyle\sum_{i,j=1}^{N/2}\left[2\int\mathrm{d}\mathbf{r}\mathrm{d}\mathbf{r}^{\prime}\psi_{i}^{*}(\mathbf{r})\psi_{j}^{*}(\mathbf{r}^{\prime})\hat{w}(\mathbf{r},\mathbf{r}^{\prime})\psi_{i}(\mathbf{r})\psi_{j}(\mathbf{r}^{\prime})\right.$
$\displaystyle\left.-\int\mathrm{d}\mathbf{r}\mathrm{d}\mathbf{r}^{\prime}\psi_{i}^{*}(\mathbf{r})\psi_{j}^{*}(\mathbf{r}^{\prime})\hat{w}(\mathbf{r},\mathbf{r}^{\prime})\psi_{j}(\mathbf{r})\psi_{i}(\mathbf{r}^{\prime})\right]$
$\displaystyle\equiv$
$\displaystyle\sum_{i,j=1}^{N/2}\left[2\braket{\psi_{i}\psi_{j}}{\hat{w}}{\psi_{i}\psi_{j}}-\braket{\psi_{i}\psi_{j}}{\hat{w}}{\psi_{j}\psi_{i}}\right].$

Here, we have introduced the notation
$\braket{\psi_{i}\psi_{j}}{\hat{o}^{(2)}}{\psi_{k}\psi_{l}}=\int\mathrm{d}\mathbf{r}\mathrm{d}\mathbf{r}^{\prime}\psi_{i}^{*}(\mathbf{r})\psi_{j}^{*}(\mathbf{r}^{\prime})\hat{o}^{(2)}(\mathbf{r},\mathbf{r}^{\prime})\psi_{i}(\mathbf{r})\psi_{j}(\mathbf{r}^{\prime})$
to abbreviate the integrals with respect to a two-body operator
$\hat{o}^{(2)}$. We call the first term of $E^{(2)}$ the _Hartree_ or mean-
field energy contribution,

$\displaystyle
E_{H}=2\sum_{i,j=1}^{N/2}\braket{\psi_{i}\psi_{j}}{\hat{w}}{\psi_{i}\psi_{j}},$
(2.35)

that can be evaluated in two steps by first calculating the Hartree potential

$\displaystyle\hat{v}_{H}(\mathbf{r})=\sum_{j=1}^{N/2}\int\mathrm{d}\mathbf{r}^{\prime}\psi_{j}^{*}(\mathbf{r}^{\prime})\hat{w}(\mathbf{r},\mathbf{r}^{\prime})\psi_{j}(\mathbf{r}^{\prime})$
(2.36)

and then the energy expression
$E_{H}=2\sum_{i=1}^{N/2}\braket{\psi_{i}}{\hat{v}_{H}}{\psi_{i}}$. The name
mean-field stems from the usual interpretation of $\braket{\hat{v}_{H}}$ as
the mean Coulomb potential of all the electrons. Importantly,
$\braket{\hat{v}_{H}}$ is equivalent to the Coulomb potential of a classical
charge distribution. The second term of $E^{(2)}$ instead is called the
_exchange_ (or Fock) contribution,

$\displaystyle
E_{X}=-\sum_{i,j=1}^{N/2}\braket{\psi_{i}\psi_{j}}{\hat{w}}{\psi_{j}\psi_{i}}$
(2.37)

that arises because of the negative terms of the Slater determinant and thus
has no classical counter-part. The name stems from the fact the two body
integrals occurring in Eq. (2.37) can be obtained from the integrals in Eq.
(2.35) by exchanging the indices on one of the two sides. The exchange
contribution represents thus a quantum-mechanical correction to the classical
or mean-field energy. Although this correction is usually much smaller than
the other energy terms in $E^{HF}$, it makes the HF description considerably
more accurate than the older Hartree theory that only takes the mean-field
into account [132].

Collecting all the terms, we get the following expression for the HF energy

$\displaystyle
E^{HF}[\\{\psi_{i}\\}]=2\sum_{i=1}^{N/2}\braket{\psi_{i}}{\hat{t}+\hat{v}}{\psi_{i}}+\sum_{i,j=1}^{N/2}\left[2\braket{\psi_{i}\psi_{j}}{\hat{w}}{\psi_{i}\psi_{j}}-\braket{\psi_{i}\psi_{j}}{\hat{w}}{\psi_{j}\psi_{i}}\right]$
(2.38)

We see that the HF energy is a functional of the $N/2$ orbitals
$\\{\psi_{1},...,\psi_{N/2}\\}$.

###### The HF minimization problem and the Fock operator

To minimize $E^{HF}$ with respect to the orbitals, we need to guarantee their
orthonormality, i.e.,

$\displaystyle\braket{\psi_{i}}{\psi_{j}}=\delta_{ij}\quad\forall\,i,j.$
(2.39)

The standard way to do so, is to constrain the minimization of the functional
by introducing a Lagrangian multiplier, $\tilde{\epsilon}_{ij}$, for every
condition. Instead of the original minimization problem of Eq. (2.32), we
minimize the Lagrangian

$\displaystyle
L^{HF}[\\{\psi_{i}\\},(\lambda_{ij})]=E^{HF}[\\{\psi_{i}\\}]-2\sum_{i,j=1}^{N/2}\tilde{\epsilon}_{ij}(\braket{\psi_{i}}{\psi_{j}}-\delta_{ij})$
(2.40)

with respect to the orbitals _and_ the Lagrange multipliers.282828Note that we
introduced the factor 2 in front of the Langrange-multiplier term for later
convenience. As it is customary, we treat the orbitals $\psi_{i}$ and their
complex conjugates $\psi^{*}_{i}$ as independent and consider
$L^{HF}[\\{\psi_{i},\psi_{i}^{*}\\},(\lambda_{ij})]$.292929It is important to
realize that our goal is to perform a free minimization over all possible
functions $f(\mathbf{r})$ that depend on one variable. In this space, there
are in principle functions with $f(\mathbf{r})\neq f^{*}(\mathbf{r})$.
However, our formalism is constructed such that the solution $f_{0}$ has the
correct property $f_{0}=f_{0}^{*}$. A necessary condition for a minimum of
$L^{HF}$ is _stationarity_ with respect to all variables

$\displaystyle 0=\delta
L^{HF}[\\{\psi_{i}\\},(\lambda_{ij})]=\sum_{i=1}^{N/2}\left[\int\mathrm{d}\mathbf{r}\frac{\delta
L^{HF}}{\delta\psi_{i}(\mathbf{r})}\delta\psi_{i}(\mathbf{r})+\int\mathrm{d}\mathbf{r}\frac{\delta
L^{HF}}{\delta\psi_{i}^{*}(\mathbf{r})}\delta\psi_{i}^{*}(\mathbf{r})\right]+\sum_{i,j=1}^{N/2}\frac{\partial
L^{HF}}{\partial\lambda_{ij}}\mathrm{d}\lambda_{ij}.$ (2.41)

This leads to three independent sets of conditions, where the latter just
gives back the orthonormality conditions of Eq. (2.39) and the former two are
equivalent to each other. For notational convenience, we utilize the second
set of conditions (with respect to $\psi_{i}^{*}$) and obtain

$\displaystyle\sum_{j=1}^{N/2}\tilde{\epsilon}_{ij}\psi_{j}(\mathbf{r})=\left[\hat{t}+\hat{v}\right]\psi_{i}(\mathbf{r})+\sum_{j=1}^{N/2}\left[2\int\mathrm{d}\mathbf{r}\psi_{j}^{*}(\mathbf{r}^{\prime})\hat{w}(\mathbf{r},\mathbf{r}^{\prime})\psi_{j}(\mathbf{r}^{\prime})\psi_{i}(\mathbf{r})-\int\mathrm{d}\mathbf{r}\psi_{j}^{*}(\mathbf{r}^{\prime})\hat{w}(\mathbf{r},\mathbf{r}^{\prime})\psi_{i}(\mathbf{r}^{\prime})\psi_{j}(\mathbf{r})\right].$
(2.42)

To simplify this expression, we can define the Coulomb-operator
$\hat{J}_{j}$303030Note that the Coulomb-operator is connected to the Hartree
potential by $\hat{v}_{H}=\sum_{j}\hat{J}_{j}$. that acts as

$\displaystyle\hat{J}_{j}\psi_{i}(\mathbf{r})$
$\displaystyle=\int\mathrm{d}\mathbf{z}^{\prime}\psi_{j}^{*}(\mathbf{r}^{\prime})w(\mathbf{r},\mathbf{r}^{\prime})\psi_{j}(\mathbf{r}^{\prime})\psi_{i}(\mathbf{r})$
(2.43) and the Exchange-operator $\hat{K}_{j}$ that acts as [160]
$\displaystyle\hat{K}_{j}\psi_{i}(\mathbf{r})$
$\displaystyle=\int\mathrm{d}\mathbf{z}^{\prime}\psi_{j}^{*}(\mathbf{r}^{\prime})w(\mathbf{r},\mathbf{r}^{\prime})\psi_{i}(\mathbf{r}^{\prime})\psi_{j}(\mathbf{r}).$
(2.44)

The right-hand side of Eq. (2.42) can then be expressed as

$\displaystyle\hat{H}^{1}\psi_{i}(\mathbf{r})=\left[\hat{t}+\hat{v}\right]\psi_{i}(\mathbf{r})+\sum_{j=1}^{N/2}\left[2\hat{J}_{j}-\hat{K}_{j}\right]\psi_{i}(\mathbf{r}),$
(2.45)

where we introduced the _Fock-operator_ $\hat{H}^{1}$, that acts on one
orbital and thus is can be seen as an effective _one-body Hamiltonian_.

###### Solving the HF equations: a nonlinear orbital eigenvalue problem

Importantly, $\hat{H}^{1}$ is hermitian and thus, we can apply a unitary
transformation to Eq. (2.42) such that
$\tilde{\epsilon}_{ij}=\epsilon_{i}\delta_{ij}$ becomes diagonal. In this
canonical formulation, Eq. (2.42) obtains the form of an _eigenvalue equation_

$\displaystyle\hat{H}^{1}\psi_{i}(\mathbf{r})=\epsilon_{i}\psi_{i}(\mathbf{r}).$
(2.46)

for HF orbitals $\psi_{i}$ and the orbital energies

$\displaystyle\epsilon_{i}=$
$\displaystyle\braket{\psi_{i}}{\hat{H}^{1}\psi_{i}}$ $\displaystyle=$
$\displaystyle\braket{\psi_{i}}{\left[\hat{t}+\hat{v}\right]\psi_{i}}+\braket{\psi_{i}}{\sum_{j=1}^{N/2}\left[2\hat{J}_{j}-\hat{K}_{j}\psi_{j}(\mathbf{r})\right]\psi_{i}}$
$\displaystyle=$
$\displaystyle\braket{\psi_{i}}{\left[\hat{t}+\hat{v}\right]\psi_{i}}+\sum_{j=1}^{N/2}\left[2\braket{\psi_{i}\psi_{j}}{\hat{w}}{\psi_{i}\psi_{j}}-\braket{\psi_{i}\psi_{j}}{\hat{w}}{\psi_{j}\psi_{i}}\right].$
(2.47)

The set of these equations for all $N/2$ orbitals are called the _Hartree-Fock
equations_ , which have to be solved (numerically) to obtain the Hartee-Fock
ground state of the Hamiltonian (2.1).

Importantly, Eq. (2.46) is not a usual, i.e., linear eigenvalue equation. The
Fock-operator $\hat{H}^{1}=\hat{H}^{1}[\\{\psi_{i}\\}]$ _depends_ on its own
eigenstates and thus the eigenvalue problem of Eq. (2.46) is in fact
_nonlinear_. We call $\hat{H}^{1}$ a nonlinear one-body operator, which has
very different properties from linear operators. We can utilize $\hat{H}^{1}$
to find the HF ground state or to determine approximate ionization potentials
and electron affinities (Koopman’s theorem [132, Sec. 10.5]). However, the
higher lying “eigenstates” of $\hat{H}^{1}$, the so-called unoccupied
orbitals, have only limited physical meaning. In fact, the validity Koopman’s
theorem is a very special feature of the Fock-matrix and similar nonlinear
one-body operators (see Sec. 2.4) in other electronic-structure theories do
not have the same level of explanatory power. In comparison, the
eigenfunctions and eigenvalues of the many-body Hamiltonian have a clear
physical interpretation: they simply denote the possible states with
corresponding energy levels that the physical system can adopt.

Another practical consequence of the nonlinearity of $\hat{H}^{1}$ is that
solving the HF equations is numerically a very challenging task that is not
comparable with an ordinary eigenvalue equation. For instance, there are
simple many-body systems, for which the Schrödinger equation is analytically
solvable.313131For example the _Harmonium_ model that consists of two
interacting electrons in a Harmonic potential. But there is _no known analytic
solution_ of the HF equations. The only known way to solve Eq. (2.46) is a so-
called _self-consistent field_ (SCF) procedure. For that, we linearize
$\hat{H}^{1}[\\{\psi_{i}\\}]\approx\hat{H}^{1}[\\{\psi_{i}^{0}\\}]$ locally
around some starting guess of orbitals $\\{\psi_{i}^{0}\\}$, solve Eq. (2.46)
for this $\hat{H}^{1}[\\{\psi_{i}^{0}\\}]$ to obtain new orbitals
$\\{\psi_{i}^{1}\\}$ and update $\hat{H}^{1}[\\{\psi_{i}^{1}\\}]$. We proceed
with this iterative procedure until the orbitals
$\\{\psi_{i}^{n+1}\\}\approx\\{\psi_{i}^{n}\\}$ do not change significantly,
i.e., until self consistence of Eq. (2.46).323232In part III, we will present
some algorithms that are capable to perform HF minimizations.

Although we know from experience that this procedure converges in most cases,
we can neither guarantee convergence, nor can we be sure that the converged
result is the global minimum of $E_{HF}$ that we are searching for. The reason
is again the nonlinearity of the problem that makes a mathematical analysis
very difficult. There are certain hints that $E_{HF}$ is not convex, which
would suggest that solutions of the HF equations could correspond to local
minima, but this has not yet been proven.333333See for example the discussion
in Ref. [161] and the references therein.

###### The cornerstone of (molecular) electronic-structure theory

Let us finish our brief discussion of HF theory with a quotation of one of the
most important textbooks on molecular electronic-structure theory:

> “The Hartree-Fock wave function is the cornerstone of ab initio electronic-
> structure theory. [… It] yields total electronic energies that are in error
> by less than 1% and a wide range of important molecular properties such as
> dipole moments, electric polarizabilities, electronic excitation energies,
> magnetizabilities, force constants and nuclear magnetic shieldings are
> usually reproduced to within 5-10 % accuracy. Molecular geometries are
> particularly well reproduced and are mostly within a few picometres of the
> true equilibrium structure.
>
> The Hartree-Fock wave function is often used in qualitative studies of
> molecular systems, particularly larger systems. Indeed, the Hartree-Fock
> wave function is still the only wave function that can be applied routinely
> to large systems, and systems containing several hundred atoms have been
> studied at this level of approximation. For accurate, quantitative studies
> of molecular systems, the Hartree-Fock wave function is by itself not useful
> but it constitutes the starting point for more accurate treatments.”
>
> (Helgaker et al., 2000 [132, p. 433])

This summarizes the performance of HF theory and its importance in the big
field of electronic-structure methods.343434Note that HF plays a less
pronounced role for extend systems, such as solids. One reason is simply that
performing HF-calculations in $k$-space is in most cases numerically very
expensive. Here, the local-density approximation of KS-DFT has a comparable
role to HF for molecular systems (see next section). Being one of the
conceptually and computationally simplest electronic-structure methods, we
will repeatedly consider HF theory in the following chapters, when we discuss
the coupled electron-photon problem.

However, as Helgaker et al. have pointed out, in many cases, HF is not
sufficient for an accurate description of many-body systems and we need more
accurate methods . In the next two sections, we will thus discuss some
approaches to go beyond HF theory and try to include the aforementioned
correlation energy in the description.

####  2.3 Density functional theory

In the last section, we have discussed how we can describe an interacting
many-body system by a single Slater determinant, i.e., an effective non-
interacting wave function. This is an approximate description that neglects
correlation, i.e., all effects that intrinsically require a multi-reference
description. However, if we give up the concept of the many-body wave function
(and with this the standard description of quantum states), we can in
principle describe any many-particle system by such a single-reference wave
function. This is called the _Kohn-Sham_ (KS) construction of DFT and it
represents for a huge class of systems the best trade-off between efficiency
and accuracy in comparison to all alternative descriptions. It is not
overstated to say that KS-DFT “is indispensable for modern quantum-chemical
modeling of materials and molecules” [162].

Because of the single-reference ansatz, KS-DFT is in practice very similar to
HF. A KS-DFT calculation means to solve the KS equations, which constitute a
nonlinear eigenvalue equation. In many cases, this is computationally even
less demanding than solving the HF equations. Conceptually however, DFT is
very distinguished from any type of wave-function method. In fact, one can see
DFT as an alternative formulation of quantum mechanics.353535For instance,
Tokatly argues that DFT can be interpreted as a _hydrodynamic formulation of
quantum mechanics_ [163, 164].

###### The density as basic variable

As the name suggests, DFT is based on the (single-particle) electron density,
that is defined as

$\displaystyle\rho(\mathbf{r})=N\sum_{\sigma_{1},...,\sigma_{N}}\int\mathrm{d}\mathbf{r}_{2}\cdots\mathrm{d}\mathbf{r}_{N}\Psi_{0}^{*}(\mathbf{x},\mathbf{x}_{2},...,\mathbf{x}_{N})\Psi_{0}(\mathbf{x},\mathbf{x}_{2},...,\mathbf{x}_{N})$
(2.48)

for an $N$-electron system with ground state wave function
$\Psi_{0}(\mathbf{x}_{1},...,\mathbf{x}_{N})$ that is the lowest eigenstate of
the electronic-structure Hamiltonian, cf. (2.1),

$\displaystyle\hat{H}=$
$\displaystyle\sum_{i=1}^{N}\left[\hat{t}(\mathbf{r}_{i})+v(\mathbf{r}_{i})\right]+\tfrac{1}{2}\sum_{i,j=1}^{N}w(\mathbf{r}_{i},\mathbf{r}_{j})$
(2.49) $\displaystyle=$ $\displaystyle\hat{T}+\hat{V}+\hat{W}.$

We renamed the three occurring terms,
$\hat{T}=\sum_{i=1}^{N}\hat{t}(\mathbf{r}_{i})$, etc. for later convenience.
Already in 1964, Walter Kohn and Pierre Hohenberg have proven that
$\rho(\mathbf{r})$ _entirely_ determines the ground state of $\hat{H}$ [165].
On a first glance, it might sound counter-intuitive that all the relevant
information of the $4N$-dimensional wave function $\Psi_{0}$ is contained in a
merely 3-dimensional quantity. However, the wave function is simply an
auxiliary quantity that allows us to formulate the standard theory of quantum
mechanics, based on linear operators. The information that determines the
physical setting in the electronic-structure Hamiltonian is actually the local
potential $v(\mathbf{r})$ that is also merely three-dimensional.

For instance, if we want to describe a diatomic molecule with nuclei at the
positions $\mathbf{R}_{1},\mathbf{R}_{2}$ and charges $Z_{1},Z_{2}$, we would
model this by considering the local potential
$v(\mathbf{r})=\frac{Z_{1}}{|\mathbf{r}-\mathbf{R}_{1}|}+\frac{Z_{2}}{|\mathbf{r}-\mathbf{R}_{2}|}$.
We can arbitrarily enlarge this system with further nuclei by adding Coulomb
terms to the local potential. But for all these systems, the kinetic and the
interaction operators are the same. These operators are thus said to be
_universal_. Any many-electron system in the non-relativistic setting is thus
completely defined by the 3-dimensional local potential and the number of
electrons $N$ and judged from this perspective, the concept of a theory based
on the electronic density might sound more reasonable.

In the realm of DFT, $v(\mathbf{r})$ is called the _external_ variable that
determines the corresponding _internal_ variable $\rho(\mathbf{r})$. This
variable pair $(\rho,v)$ is defined by the structure of the Hamiltonian
(2.49). To see this, let us calculate the expectation value of the potential
term

$\displaystyle\braket{\hat{V}}=$
$\displaystyle\braket{\Psi}{\sum_{i=1}^{N}v(\mathbf{r}_{i})\Psi}$
$\displaystyle=$
$\displaystyle\sum_{\sigma_{1},...,\sigma_{N}}\int\mathrm{d}\mathbf{r}_{1}\mathbf{r}_{2}\cdots\mathrm{d}\mathbf{r}_{N}\Psi^{*}(\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N})\left[\sum_{i=1}^{N}v(\mathbf{r}_{i})\right]\Psi(\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N})$
$\displaystyle=$
$\displaystyle\int\mathrm{d}\mathbf{r}v(\mathbf{r})N\sum_{\sigma_{1},...,\sigma_{N}}\int\mathbf{r}_{2}\cdots\mathrm{d}\mathbf{r}_{N}\Psi^{*}(\mathbf{x},\mathbf{x}_{2},...,\mathbf{x}_{N})\Psi(\mathbf{x},\mathbf{x}_{2},...,\mathbf{x}_{N})$
$\displaystyle=$
$\displaystyle\int\mathrm{d}\mathbf{r}v(\mathbf{r})\rho(\mathbf{r}),$

where we used the exchange symmetry of $\Psi$ from the second to the third
line. This term _defines_ the pair of external and internal variable and it is
at the center of the Hohenberg-Kohn proof [166] of the existence of the one-
to-one mapping

$\displaystyle
v(\mathbf{r})\overset{1:1}{\longleftrightarrow}\rho(\mathbf{r}).$ (2.50)

The simplicity of this argument illustrates the fact that there are many
possible functional theories, depending on the structure of the Hamiltonian. A
famous example is _current-density functional theory_ [167] that can describe
systems with external magnetic field or _reduced density matrix functional
theory_ that we will introduce in Sec. 2.4.2. All these functional theories
are based on the same type of proof,363636For very general proof of this
statement, see Penz et al. [166]. which we do not recapitulate at this point.
We will instead present (one version of) the proof in the realm of coupled
electron-photon systems in Sec. 3.3.2.

###### The universal density functional

The Hohenberg-Kohn theorem establishes DFT with the density $\rho$ as the
fundamental variable to describe many-electron systems. We thus can replace
the wave function in the variational principle, cf. Eq. (2.2), by the density
and write

$\displaystyle E_{0}=\inf_{\rho}E[\rho],$ (2.51)

where $E[\rho]$ describes the energy of the system with Hamiltonian (2.49) as
a functional of the density $\rho$. We can further specify this expression,
because also the ground state wave function $\Psi=\Psi[\rho]$ is uniquely
connected to $\rho$ and thus we have

$\displaystyle\begin{split}E[\rho]=&\braket{\Psi[\rho]}{\hat{H}\Psi[\rho]}\\\
=&\braket{\Psi[\rho]}{\hat{T}\Psi[\rho]}+\braket{\Psi[\rho]}{\hat{V}\Psi[\rho]}+\braket{\Psi[\rho]}{\hat{W}\Psi[\rho]}\\\
=&T[\rho]+\int\mathrm{d}\mathbf{r}v(\mathbf{r})\rho(\mathbf{r})+W[\rho].\end{split}$
(2.52)

We used here that the potential energy
$\braket{\hat{V}}=\int\mathrm{d}\mathbf{r}v(\mathbf{r})\rho(\mathbf{r})$ can
be explicitly written in terms of $\rho$. The remaining part of $E[\rho]$ is
called the _universal functional_ 373737For details on mathematical issues of
the universal functional as defined by Hohenberg and Kohn and an alternative
more general definition, the reader is referred to the publications by Levy
[168] and Lieb [169].

$\displaystyle F[\rho]=\hat{T}[\rho]+\hat{W}[\rho].$ (2.53)

This reformulation has however a _fundamental problem_ : Although the
Hohenberg-Kohn theorem proves the existence of $F$, it does _not provide an
explicit expression_ for it, because the proof is not constructive (see Sec.
3.3.2). As a matter of fact, most of the research in the field of DFT aims on
understanding how this functional looks like. This is a very difficult task
and only many years of research lead to the remarkable accuracy of state-of-
the-art DFT.

###### The Kohn-Sham construction

The first crucial step toward a useful density-functional is the KS
construction [170] that makes use of the increased flexibility of the density
description. The basic idea is to employ as in HF, a non-interacting wave
function, the KS wave function $\Psi^{s}$ (we denote all quantities in the
auxiliary system with an s). We have seen in HF theory, that this is a very
good starting point to describe the interacting problem and at the same time,
a single Slater determinant can be calculated relatively cheaply by effective
orbital equations. However, in contrast to HF theory, the fundamental variable
in DFT is the density and not the wave function. Thus, the KS wave function is
a purely auxiliary quantity that in principle is _not related_ to the many-
body wave function of the system. The connection between the non-interacting
_Kohn-Sham system_ described by the Hamiltonian

$\displaystyle\hat{H}^{s}=\sum_{i=1}^{N}\hat{t}(\mathbf{r}_{i})+v^{s}(\mathbf{r}_{i})$
(2.54)

and the physical system described by Hamiltonian (2.49) is given by their
densities, i.e., we demand

$\displaystyle\rho^{s}(\mathbf{r})=\rho(\mathbf{r}).$ (2.55)

Since the connection between density and potential is unique, i.e., there is
one and only one interacting system with the density $\rho$, there is also one
and only one non-interacting system with the density $\rho^{s}=\rho$. Thus,
the Hohenberg-Kohn theorem also proves that the KS potential is unique.

###### Connecting the physical and the KS system

The original construction due to Kohn and Sham [170] to connect the KS and the
physical system is based on the energy. For that we define the Kohn-Sham
energy (using Eq. (2.55))

$\displaystyle
E^{s}[\rho]=T^{s}[\rho]+\int\mathrm{d}\mathbf{r}v^{s}(\mathbf{r})\rho(\mathbf{r})$
(2.56)

and then re-express the physical (many-body) energy as

$\displaystyle\begin{split}E[\rho]=&\int\mathrm{d}\mathbf{r}v(\mathbf{r})\rho(\mathbf{r})+F[\rho]\\\
=&T^{s}[\rho]+\int\mathrm{d}\mathbf{r}v(\mathbf{r})\rho(\mathbf{r})+F[\rho]-T^{s}[\rho]\\\
\equiv&T^{s}[\rho]+\int\mathrm{d}\mathbf{r}v(\mathbf{r})\rho(\mathbf{r})+E_{Hxc}[\rho].\end{split}$
(2.57)

This defines the Hartree-exchange-correlation energy

$\displaystyle E_{Hxc}[\rho]=$ $\displaystyle T[\rho]-T^{s}[\rho]+W[\rho]$
$\displaystyle=$ $\displaystyle E_{H}[\rho]+E_{xc}[\rho],$ (2.58)

that usually is further divided in the mean-field or Hartree part
$E_{H}[\rho]=\frac{1}{2}\int\mathrm{d}\mathbf{r}\mathbf{r}^{\prime}\rho(\mathbf{r})w(\mathbf{r},\mathbf{r}^{\prime})\rho(\mathbf{r}^{\prime})$,
that is an explicit functional of $\rho$, and the remaining unknown part
$E_{xc}$, that is called _exchange-correlation energy_ (in analogy to the HF
energy definitions, see Sec. 2.2). If we assume that this functional is
differentiable with respect to $\rho$-variations,383838In fact, it can be
shown that this functional is not differentiable. There are ways to circumvent
this problem (see, e.g., [166]), but from a mathematical point of view this
issue is not yet resolved. In practice, this seems to play only a secondary
role. However, one should not underrate such mathematical issues. Many
advances in DFT were in fact triggered by mathematical research such as the
Levy-Lieb constrained search formulation [168, 171, 169]. For a good
introduction in this topic from a physicist’s perspective, the reader is
referred to the respective section in the textbook by Dreizler and Gross
[172]. we can define the ground state by the stationarity condition

$\displaystyle 0$ $\displaystyle=\frac{\delta E[\rho]}{\delta\rho}$
$\displaystyle=\frac{\delta
T^{s}[\rho]}{\delta\rho}+v(\mathbf{r})+\frac{\delta
E_{H}[\rho]}{\delta\rho}+\frac{\delta E_{xc}[\rho]}{\delta\rho}.$

The first term

$\displaystyle\frac{\delta T^{s}[\rho]}{\delta\rho}=-v^{s}[\rho]$

can be calculated by Eq. (2.56) (see [172, part 4.1] for details), the third
term is just the Hartree-potential identical to the definition of HF theory,

$\displaystyle
v_{H}[\rho]=\int\mathrm{d}\mathbf{r}^{\prime}\rho(\mathbf{r}^{\prime})w(\mathbf{r},\mathbf{r}^{\prime}),$
(cf. 2.36)

and the last term defines the unknown _exchange-correlation potential_

$\displaystyle v_{xc}[\rho]=\frac{\delta E_{Hxc}[\rho]}{\delta\rho}.$ (2.59)

It follows that the potential of the auxiliary system is

$\displaystyle
v^{s}[\rho](\mathbf{r})=v(\mathbf{r})+v_{H}[\rho](\mathbf{r})+v_{xc}[\rho](\mathbf{r}).$
(2.60)

Under our assumption that $E[\rho]$ is differentiable, this relation uniquely
defines the KS system. We want to remark at this point that there are
alternative constructions to define the KS system avoid certain mathematical
subtleties. For instance via the force-balance equations [173] or by
considering more general definition spaces [166].

###### The KS equations

Since $\hat{H}^{s}$ does not contain interaction terms, its ground state is
(assuming no degeneracy393939This assumption makes the derivations in the
following simpler, but it is not necessary for the KS construction. For
details, see Ref. [172].) described by one Slater determinant

$\displaystyle\Psi^{s}(\mathbf{x}_{1},...,\mathbf{x}_{N})=|\psi_{1}\cdots\psi_{N}|_{-},$
(2.61)

which can be determined by a set of orbital equations

$\displaystyle\left[\hat{t}(\mathbf{r})+v^{s}(\mathbf{r})\right]\psi_{i}(\mathbf{x})=\epsilon_{i}\psi_{i}(\mathbf{x}).$
(2.62)

These are called the _KS equations_. We see that in contrast to HF theory, the
KS equations do not include a nonlocal term like the Exchange operator (cf.
Eq.(2.44)) and are thus formally simpler. However, approximations to the
unknown exchange correlation potential usually have strong nonlinear
dependencies on the density, or derivatives of the density. A very
sophisticated class of approximate functionals, so-called _hybrids_ , even
lift the assumption of a local $v^{s}$ completely and include the HF exchange
energy.404040See for example the review by Perdew and Kurth [174] for an
overview on the main classes of known functionals. Nevertheless, the KS
equations define a nonlinear eigenvalue problem for the effective one-body
Hamiltonian $\hat{H}^{s}$. This is very similar to HF theory and accordingly
the solution strategies are basically the same.414141Although the locality of
most exchange-correlation potentials makes a large difference in practical
calculations. Especially, KS-DFT can be applied to condensed matter systems
considerably more easily than HF. For example KS-DFT calculations also require
an SCF procedure (see Sec. 2.2) that is called the _KS scheme_.

###### The quest for the universal functional: approximations to the exchange-
correlation potential

Although the KS scheme simplifies the functional construction, we have not
gained much so far with respect to HF theory, because we do not know the form
of the exchange-correlation functional. The formalism itself does not give us
any hint, how $v_{xc}$ could look like and thus, there is literally nothing
else than “intuition” or” experience” to construct functionals. The
historically first functional has been proposed by Kohn and Sham in Ref.
[170], where they introduced the KS scheme. This _local-density approximation_
(LDA) is based on one of the few many-electron systems that can be solved
analytically, i.e., the homogeneous electron gas (HEG). Importantly, the
electron-density of the HEG $\rho^{HEG}$ is constant and thus, we can
explicitly parametrize the total energy $E^{HEG}(\rho^{HEG})$ as a function of
the density.424242Note that $E^{HEG}(\rho^{HEG})$ is really a function and not
a functional, because $\rho^{HEG}$ is just a number. Utilizing definition
(2.57), we can then (analytically) derive the exchange-correlation energy
expression [175]

$\displaystyle E_{xc}^{HEG}(\rho)\equiv
E_{X}^{LDA}[\rho]=-\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\int\mathrm{d}\mathbf{r}\rho^{4/3}(\mathbf{r}),$
(2.63)

where we formally reintroduced the spatial dependence of
$\rho(\mathbf{r})=\rho^{HEG}=const.$ Additionally, we removed the index “C,”
because a non-interacting system as the HEG per definition does not have
correlation. This procedure can be generalized to the interacting homogeneous
electron gas, yielding a functional $E_{xc}^{LDA}$ that includes correlation
effects.434343This model cannot be solved analytically anymore, but it still
can be reduced to the numerical calculation of one integral in the many-body
coordinate space. This can be done very efficiently by Monte-Carlo
integration. For an exhaustive discussion on the homogeneous electron gas, the
reader is referred to Ref. [176]. The local-density approximation is then
simply to apply $E_{xc}^{LDA}$ to systems that is not homogeneous, i.e., to
consider $\rho=\rho(\mathbf{r})$ as we have already anticipated in Eq. (2.63).
Despite its simplicity, the LDA is impressively accurate for a large range of
systems, especially in the realm of condensed matter, where the LDA has a
similar significance as HF for molecular systems. But also for certain classes
of chemical systems, the LDA is very accurate and often also better than HF
[174]. However, the LDA is by far not sufficient to reach “chemical accuracy,”
i.e., the accuracy required in electronic-structure calculations to make
realistic predictions.444444See for example the Nobel-lecture of Pople [177].
This has only been achieved with the historically second class of functionals,
which are the _generalized gradient approximations_ (GGAs). The GGA
functionals are conceptually a straightforward generalization of the LDA,
because they consider besides $\rho$ also its gradient, i.e.,

$\displaystyle E_{xc}^{GGA}=E_{xc}^{GGA}[\rho,\nabla\rho].$ (2.64)

However, from the definition of the LDA in 1966 until 1985, when Perdew [178]
constructed the first GGA functional that provided significantly better
results, passed _almost 20 years_. This illustrates how difficult functional
construction in DFT really is. However, once an accurate functional is
constructed, its potential is enormous because the the KS equations that have
to be solved remain basically the same. Thus it is probably no wonder that the
publication of Perdew marks a turning point in the history of DFT. Since 1985,
the number of proposed exchange-correlation functionals and publications has
grown exponentially [179] and started the “incredible success story” [180] of
DFT. Modern functionals include, e.g., even higher derivatives of the density
(for example the _meta-GGAs_) and the Fock or exchange contribution of HF
(_hybrid functionals_) [174, 172]. Such theoretical ideas define the principal
form of the functionals, but the precise contribution of the different
theoretical levels usually cannot be predicted. Many successful KS-DFT
functionals have thus been constructed by fitting to experimental
data.454545The interested reader is referred to the recent report by Medvedev
et al. [162] that highlights the history of functional construction and
reflects critically on the strong focus on fitting that has become common
practice nowadays.

Despite these huge efforts and certain systematic improvements, it is
difficult to compare and judge the performance of the many existing
functionals. Theoretically, there is one universal functional and one could
expect that the known functionals somehow “converge” into one direction, i.e.,
they approximate increasingly better the exact functional. However, the
opposite seems to be true: there is basically no functional that is accurate
for all material classes, but instead, there are many specific problems that
require specific functionals.464646See for example the extensive assessment of
density functionals in quantum chemistry by Mardirossian and Head-Gordon [84].
Consequently, applying KS-DFT in practice is everything but easy and requires
experience and a careful literature research. Discussing functional
construction in more depth is beyond the scope of this text and thus, we refer
the reader for further details to the already cited literature. To name a few
further examples, the review of Jones [179] provides a very good introduction
to the topic of density functionals including its very interesting history and
the review by Burke [180] provides a perspective on the recent state of DFT.

As a last example for an exchange correlation functional, we want to mention
the _exact exchange_ (EXX) approximation which will play an important role in
the generalization of DFT to coupled electron-photon systems (see Sec. 3.2).
The EXX approximation is exceptional in comparison to most other known KS-DFT
functionals, because it does not rely on empirical (or in the special case of
the LDA analytical) knowledge. The idea is to approximate the interacting
system with a single-reference wave function exactly as in HF. However, the
EXX functional is _local_ , which requires to solve an auxiliary equation
inside the SCF routine. This is a non-unique problem and thus there are
several possibilities to define this auxiliary equation. The most common way
is the _optimized effective potential_ (OEP) approach together with the
Krieger-Li-Iafrate (KLI) approximation [172], but there are also other
possible constructions [181]. Importantly, the EXX functional is usually more
accurate than HF but at the same time as generic. This makes it an important
tool for developing first-principles methods in settings that are not covered
by the electronic-structure Hamiltonian.

We conclude our brief survey on KS-DFT functionals with a small example that
shows the significance of density functionals in modern material science and
chemistry. Despite all the problems and difficulties with KS-DFT, accurate
functionals are known nowadays and they are basically the only tool that
exists to calculate the properties of many realistic matter systems. This is
well-reflected in the extraordinary number of citations of the corresponding
publications. The most famous example is the paper [182] in which Perdew,
Burke and Enzerhof propose their PBE functional, which has become within
merely two decades the _most cited paper of all physics_. According to Google
scholar, there are more than 110000 publications that refer to this
publication.

###### A remark on the KS scheme and alternatives

With the KS scheme, we have reduced the problem of calculating the
$4N$-dimensional wave function on finding the $N$ lowest eigenstates of the
nonlinear one-particle operator $\hat{H}^{s}$. This is on obvious
computational advantage over the exponential wall of the many-body problem.
Still, it is much more expensive than our starting point, i.e., a direct
minimization of $E[\rho]$ (Eq. (2.51)). In fact, people have tried to
construct theories that are directly based on $\rho$ (nowadays known under the
term _orbital-free_ DFT) even before DFT was developed. A famous example is
_Thomas-Fermi_ theory,474747See for example the review by Jones [179] and the
references therein. developed already in 1927, which provides however a very
poor description of quantum systems. For example, matter is not stable in
Thomas-Fermi theory [183]. There have been advances in the field, but still
the accuracy of known explicit density functionals is very low.484848See Ref.
[184] for a recent review on the topic. The crucial point of the KS
construction is that the contribution of the unknown $v_{xc}[\rho]$ to the
total energy is _considerably smaller_ than the contribution of the universal
functional $F[\rho]$ (cf. Eq. (2.52)). The reason is that the Slater-
determinant includes already a very large part of the quantum mechanical
problem into our description (as we know from HF theory). Thus, DFT is almost
exclusively utilized in the KS picture and often DFT and KS-DFT are used
interchangeably.

###### When DFT fails: the strong-correlation regime

In this last paragraph, we want to discuss the (actual) limitations of KS-DFT
(and other single-reference methods). Although density functionals can in
principle exactly describe every physical setting that is included in the
electronic-structure Hamiltonian (Eq. (2.1)), there are scenarios, where all
known KS-DFT functionals are severely inaccurate. The research on functional
construction has shown that in such settings, the form of the exact functional
is very intricate. For instance, the exact functional of the non-interacting
HEG is given by a formula as simple as Eq. (2.63). Systems that do not deviate
too much from the HEG like typical conductors are usually very well described
by KS-DFT. However, for other seemingly simple systems, such as two separate
Hydrogen atoms, it turns our that the exact functional has a very complicated
form.

The principal reason for this difference is the single Slater determinant
ansatz in the KS system that is a very good starting point for the HEG, but a
very bad starting point for two separate atoms. The (spatial part of the) KS
wave function for the ground state of two electrons reads simply

$\displaystyle\Psi^{s}(\mathbf{r}_{1},\mathbf{r}_{2})=\psi(\mathbf{r}_{1})\psi(\mathbf{r}_{2}),$
(2.65)

where $\psi$ is the doubly-occupied spatial orbital of a spin-
singlet.494949Ground-states are normally singlets, i.e., spin-saturated. The
exact ground state of two separate Hydrogen atoms after a proper
symmetrization reads instead

$\displaystyle\Psi(\mathbf{r}_{1},\mathbf{r}_{2})=$
$\displaystyle\tfrac{1}{\sqrt{2}}\left(\psi_{a}(\mathbf{r}_{1})\psi_{b}(\mathbf{r}_{2})+\psi_{a}(\mathbf{r}_{2})\psi_{b}(\mathbf{r}_{1})\right),$
(2.66)

where $\psi_{a(b)}=\exp(-(r-R_{a(b)}))/\sqrt{\pi}$ is the ground state of the
Hydrogen located at $R_{a}(R_{b})$, where $|R_{a}-R_{b}|>>1$ such that the
orbitals do not overlap. Also $\Psi$ describes a spin-singlet,505050Note that
the antisymmetry is contained in the spin-function:
$\tfrac{1}{\sqrt{2}}(\alpha(\sigma_{1})\beta(\sigma_{2})-\alpha(\sigma_{2})\beta(\sigma_{1}))$.
but it requires two _different_ spatial orbitals, which are degenerate, i.e.,
they have the same energy eigenvalue. From the wave function perspective, this
special situation requires thus a multi-reference ansatz, i.e., a wave