3.1.1 ($\tilde{\Phi}$ in Eq. (3.10)). Thus, we consider an $N_{b}$ boson state
in the representation $\psi_{b}(\alpha_{1},...,\alpha_{N_{b}})$ and define the
corresponding bosonic $q$RDM

$\displaystyle\Gamma_{b}^{(q)}(\alpha_{1},...,\alpha_{q};$
$\displaystyle\alpha_{1}^{\prime},...\alpha_{q}^{\prime})=\tfrac{N_{b}!}{(N_{b}-q)!}\sum_{\alpha_{q+1},...,\alpha_{N_{b}}=1}^{M}$
$\displaystyle\psi_{b}^{*}(\alpha_{1}^{\prime},...,\alpha_{q}^{\prime},\alpha_{q+1},...,\alpha_{N_{b}})\psi_{b}(\alpha_{1},...,\alpha_{q},\alpha_{q+1},...,\alpha_{N_{b}}).$
(3.41)

Equivalently to the electronic case, we denote the 1RDM by
$\gamma_{b}(\alpha,\beta)=\Gamma_{b}^{(1)}(\alpha;\beta)$. But this definition
is not sufficient to describe expectation values of the Hamiltonian (1.37).
Because of the form of the electron-photon interaction, the number of photons
is undetermined and we need to work with Fock-space wave functions
$\ket{\Phi}$. For the 1RDM, we generalize the definition (3.3.1) as

$\displaystyle\gamma_{b}(\alpha,\beta)=\braket{\Phi}{\hat{a}^{\dagger}_{\beta}\hat{a}_{\alpha}\Phi}=\sum_{N_{b}=0}^{\infty}N_{b}\left(\sum_{\alpha_{2},...,\alpha_{N_{b}}=1}^{M}\psi_{b}^{*}(\beta,\alpha_{2},...,\alpha_{N_{b}})\psi_{b}(\alpha,\alpha_{2},...,\alpha_{N_{b}})\right).$
(3.42)

We see here why the coordinate-representation is normally not used for photon
states. Since their particle-number is usually not fixed, the description
becomes highly cumbersome. The abstract definition of $\ket{\Phi}$ together
with the annihilation (creation) operators $\hat{a}_{\alpha}^{(\dagger)}$ are
much simpler. Thus, we consider in the following the general bosonic Fock-
space $q$RDM directly via

$\displaystyle\Gamma^{(q)}_{b}(\alpha_{1},...,\alpha_{q};\alpha_{1}^{\prime},...,\alpha_{q}^{\prime})=\braket{\Phi}{\hat{a}^{\dagger}_{\alpha_{1}^{\prime}}\cdots\hat{a}^{\dagger}_{\alpha_{q}^{\prime}}\hat{a}^{\vphantom{\dagger}}_{\alpha_{q}}\cdots\hat{a}^{\vphantom{\dagger}}_{\alpha_{1}}\Phi}.$
(3.43)

The fermionic and bosonic RDMs can be extended to the coupled fermion-boson
case straightforwardly by just integrating/summing out the other degrees of
freedom. That is, if we have a general electron-boson state of the form of Eq.
(3.2) we can accordingly define
$\Gamma_{e}^{(q)}\equiv\frac{N!}{(N-q)!}\sum_{\sigma_{1},...,\sigma_{N}}\int\mathrm{d}^{3(N-q)}r\;\mathrm{d}^{M}p\;\Psi^{*}\Psi$
as well as
$\Gamma^{(q)}_{b}\equiv\braket{\Psi}{\hat{a}^{\dagger}_{\alpha_{1}^{\prime}}\cdots\hat{a}^{\dagger}_{\alpha_{q}^{\prime}}\hat{a}_{\alpha_{q}}\cdots\hat{a}_{\alpha_{1}}\Psi}$.191919Note
that to actually calculate the latter case, we have to employ the connection
(3.9).

In a next step, we see whether these standard ingredients of RDM theories are
sufficient to express the energy expectation value of the Hamiltonian of Eq.
(3.1). For the purely electronic part, the different contributions can be
expressed either explicitly by the electronic 1RDM or by the electronic 2RDM.
As the general prescription of Def. (2.77) tells us, all expectation values of
the single-particle operators of $\hat{H}$ are given in terms of the 1RDM.
These are the standard electronic operators $\hat{T}$ and $\hat{V}$ but also
the single-particle part of the dipole self-energy $\hat{H}_{self,1}$. We can
thus write

$\displaystyle T[\gamma_{e}]$
$\displaystyle=\int\mathrm{d}^{3}r\left[-\tfrac{1}{2}\nabla_{\mathbf{r}}^{2}\right]\gamma_{e}(\mathbf{r};\mathbf{r}^{\prime})|_{\mathbf{r}^{\prime}=\mathbf{r}},$
(3.44a) $\displaystyle V[\gamma_{e}]$
$\displaystyle=\int\mathrm{d}^{3}r\left[\vphantom{\sum}v(\mathbf{r})\right]\gamma_{e}(\mathbf{r};\mathbf{r}),$
(3.44b) $\displaystyle H_{self,1}[\gamma_{e}]$
$\displaystyle=\int\mathrm{d}^{3}r\left[\sum_{\alpha=1}^{M}\tfrac{1}{2}(\bm{\lambda}_{\alpha}\cdot\mathbf{r})^{2}\right]\gamma_{e}(\mathbf{r};\mathbf{r})$
(3.44c)

as functionals of $\gamma_{e}$. The latter two energies are actually
functionals of merely the electronic density
$\rho(\mathbf{r})=\gamma_{e}(\mathbf{r};\mathbf{r})$. As before, the subscript
$|_{\mathbf{r}^{\prime}=\mathbf{r}}$ indicates that $\mathbf{r}^{\prime}$ is
set to $\mathbf{r}$ after the application of the semi-local single-particle
operator $-\tfrac{1}{2}\nabla_{\mathbf{r}}^{2}$. The expectation value of the
electronic interaction energy $\hat{W}$ and the two-body part of the dipole
self-energy $\hat{H}_{d}^{(2)}$ are given in terms of the (diagonal) of the
2RDM by

$\displaystyle W[\Gamma^{(2)}_{e}]=$
$\displaystyle\tfrac{1}{2}\\!\int\mathrm{d}^{3}r\mathrm{d}^{3}r^{\prime}w(\mathbf{r},\mathbf{r}^{\prime})\Gamma_{e}^{(2)}(\mathbf{r},\mathbf{r}^{\prime};\mathbf{r},\mathbf{r}^{\prime}),$
(3.45a) $\displaystyle H_{d}^{(2)}[\Gamma^{(2)}_{e}]=$
$\displaystyle\tfrac{1}{2}\\!\int\mathrm{d}^{3}r\mathrm{d}^{3}r^{\prime}\left[\sum_{\alpha=1}^{M}\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}\right)\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}^{\prime}\right)\right]\Gamma_{e}^{(2)}(\mathbf{r},\mathbf{r}^{\prime};\mathbf{r},\mathbf{r}^{\prime}).$
(3.45b)

Hence, for the electronic operator expectation values little changes in
comparison to a purely fermionic problem, except that we have a coupled
electron-boson wave function and the extra contributions of the dipole self-
energy. For the purely bosonic part of the coupled Hamiltonian we can do the
same and find (see App. A.1 for further details about the Hamiltonian for
different representations)

$\displaystyle H_{ph}[\gamma_{b}]$
$\displaystyle=\braket{\Psi}{\left\\{\sum_{\alpha=1}^{M}\left[-\tfrac{1}{2}\tfrac{\partial^{2}}{\partial
p_{\alpha}^{2}}+\tfrac{\omega^{2}_{\alpha}}{2}p_{\alpha}^{2}\right]\right\\}\Psi}=\sum_{\alpha=1}^{M}\left(\omega_{\alpha}+\frac{1}{2}\right)\gamma_{b}(\alpha,\alpha).$
(3.46)

We see that the photon-energy functional resembles structurally exactly the
electronic one-body expressions (3.44a). This is of course due to the
construction of RDMs, but nevertheless it is remarkable, since we deal with
two different particle classes. The underlying exchange symmetry of an RDM is
instead encoded in what we introduced in Sec.2.4.1 as $N$-representability
conditions. For the 1RDM of (fermionic or bosonic) ensembles, the conditions
are simple and they are most easily expressed in the eigenbases
$\gamma_{e/b}=\sum n_{i}^{e/b}\left(\phi^{i}_{e/b}\right)^{*}\phi^{i}_{e/b}$,
where the $\phi^{i}_{e/b}$ are called the natural orbitals and the
$n^{i}_{e/b}$ the natural occupation numbers. Then, the conditions are

$\displaystyle 0\leq\,$ $\displaystyle n_{i}^{e}\leq 1,$ (3.47a)
$\displaystyle 0\leq\,$ $\displaystyle n_{i}^{b},$ (3.47b)

for fermions and bosons, respectively. Let us recall our example of the
$n=n_{1}+...+n_{M}$-boson state that describes, e.g., photons occupying mode 1
with $n_{1}$, mode 2 with $n_{2}$ photons etc. From the definition of the
bosonic 1RDM it is obvious that $n_{i}\equiv n^{b}_{i}$ in the above
expression and thus, we see how the bosonic character of photons is directly
transferred to the N-representability condition. In practice, this means that
_any_ positive-semidefinite matrix can be a bosonic 1RDM, which makes the
conditions much less stringent than in the fermionic case, where the Pauli-
principle translates to the upper bound for the $n_{i}^{e}$. Additionally, if
the particle number $N_{e/b}$ of one species of the system is conserved, the
respective sum-rule

$\displaystyle\sum_{i=1}^{\infty}n_{i}^{e/b}=N_{e/b}$ (3.48)

becomes a second part of the $N$-representability conditions. Thus, in our
case, we have yet another condition for the electronic part of the system, but
no further bounds for the photons. We mention this here to stress that the
1RDM $N$-representability conditions for electrons reduce the configuration
space of valid 1RDMs quite strongly. If we want to build some kind of RDMFT,
this is on the one hand numerically challenging, because we have to test these
conditions, but on the other hand it is really helpful, because then any
approximation to the interaction functional needs to carry less “information.”
We discussed in Sec. 2.4 a similar example of the more stringent
$N$-representability conditions that refer to pure-states and not only
ensembles. Theophilou et al. [206] showed that enforcing the pure-state
conditions in RDMFT yield more accurate results than enforcing only the
ensemble conditions for the same functional. This indicates that the
functional construction in a theory that is based on the photonic 1RDM might
be more difficult.

###### The problem of the electron-photon coupling: we need a new type of RDM

However, the coupled electron-photon theory provides us even more intricate
problems in terms of RDMs and representability conditions. The bilinear
coupling term $\hat{H}_{int}$, which is the key quantity of the coupled theory
cannot be treated by the qRDMs, that we have defined in Def. (2.77) and Eq.
(3.43). Formally, we could write

$\displaystyle H_{I}[\Gamma^{(3/2)}_{e,b}]=$
$\displaystyle\braket{\Psi}{\left[\sum_{\alpha=1}^{M}-\omega_{\alpha}p_{\alpha}\bm{\lambda}_{\alpha}\cdot\hat{\bf
D}\right]\Psi}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{M}-\omega_{\alpha}\braket{\Psi}{\left[\sqrt{\frac{2}{\omega_{\alpha}}}\left(\hat{a}^{\dagger}_{\alpha}+\hat{a}_{\alpha}\right)\bm{\lambda}_{\alpha}\cdot\hat{\bf
D}\right]\Psi},$ (3.49)

with a new reduced quantity that mixes light and matter degrees of freedom and
can be interpreted as a “$3/2$-body” operator
$\Gamma^{(3/2)}_{e,b}(\alpha;\mathbf{r},\mathbf{r}^{\prime})$. This can be
best understood, if we also lift the continuous fermionic problem into its own
Fock space and introduce genuine field operators
$\hat{\psi}_{e}^{\dagger}(\mathbf{r}\sigma)$ and
$\hat{\psi}_{e}(\mathbf{r}\sigma)$ with the usual anti-commutation relations.
Similarly to the discussed bosonic case, the electronic RDMs can then be
written in terms of strings of creation and annihilation field operators
[254]. We re-express

$\displaystyle\braket{\Psi}{\left[\left(\hat{a}^{\dagger}_{\alpha}\dagger\hat{a}_{\alpha}\right)\bm{\lambda}_{\alpha}\cdot\hat{\bf
D}\right]\Psi}=\sum_{\sigma}\int\mathrm{d}^{3}r\braket{\Psi}{\left[\left(\hat{a}^{+}_{\alpha}+\hat{a}_{\alpha}\right)\hat{\psi}_{e}^{\dagger}(\mathbf{r}\sigma)\hat{\psi}_{e}(\mathbf{r}\sigma)\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}\right)\right]\Psi}$

and define

$\displaystyle\Gamma^{(3/2)}_{e,b}(\alpha;\mathbf{r},\mathbf{r}^{\prime})=\sum_{\sigma}\braket{\Psi}{\left[\left(\hat{a}^{\dagger}_{\alpha}+\hat{a}_{\alpha}\right)\hat{\psi}_{e}^{\dagger}(\mathbf{r}\sigma)\hat{\psi}_{e}(\mathbf{r}^{\prime}\sigma)\right]\Psi}.$
(3.50)

We can now re-write Eq. 3.3.1

$\displaystyle
H_{I}[\Gamma^{(3/2)}_{e,b}]=\sum_{\alpha=1}^{M}-\sqrt{2\omega_{\alpha}}\int\mathrm{d}^{3}r\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}\right)\Gamma^{(3/2)}_{e,b}(\alpha;\mathbf{r},\mathbf{r}).$
(3.51)

The “property” of the bilinear interaction term to create/annihilate bosons by
interacting with the electronic subsystem is thus directly connected to the
half-integer index of the corresponding RDM. This is clearly the straight-
forward generalization of the $q$RDMs. However, this definition is only useful
if we understand its properties. First of all, we note that the $3/2$-body RDM
has in general no simple connection to any $q$RDM, similarly to the connection
formula (2.84), even if we extend the definitions to include combined matter-
boson $q$RDMs. If we consider for example the integration
$\int\mathrm{d}\mathbf{r}\Gamma^{(3/2)}_{e,b}(\alpha;\mathbf{r},\mathbf{r})$,
we will have a function that depends only on $\alpha$ and thus will be related
to some kind of a bosonic $1/2$-body RDM. We refrain at this point from a
precise definition of such half-body objects since we will not further discuss
them. However, one obvious reason for this is that $q$RDMs conserve particle
numbers, while half-body RDMs do not.

We want to finish this discussion with a simple example, that shows that the
information of the photonic 1RDM
$\gamma_{b}(\alpha,\beta)=\braket{\Phi}{\hat{a}^{\dagger}_{\beta}\hat{a}_{\alpha}\Phi}$
is different from the photonic $1/2$RDM
$\Gamma^{(1/2)}_{b}(\alpha)=\braket{\Phi}{\hat{a}^{\dagger}_{\beta}\Phi}$. In
the special case that $\ket{\Phi}$ consists only of coherent states for each
mode (which essentially means that we have treated the photons in mean field)
and since the coherent states are eigenfunctions to the annihilation
operators, we find $\gamma_{b}(\alpha,\beta)=d_{\beta}^{*}d_{\alpha}$, where
$d_{\alpha}$ is the total displacement of the coherent state of mode $\alpha$.
In this case we also know
$\braket{\Phi}{\hat{a}^{\dagger}_{\beta}\Phi}=d^{*}_{\beta}$. If we now assume
all but one mode, say mode 1, having zero displacement, then we only know
$\gamma_{b}(1,1)=|d_{1}|^{2}$ from the bosonic 1RDM. We do, however, in
general not know what $d_{1}^{*}$ is. For other states, such a connection is
even less explicit.

###### Variational RDM theory for NR-QED: The representability problem

Putting the interrelations among the different RDMs aside for the moment, the
minimization for the coupled matter-boson problem can be reformulated by

$\displaystyle E_{0}=$ $\displaystyle\inf_{\Psi}\braket{\Psi}{\hat{H}\Psi}$
$\displaystyle=$
$\displaystyle\inf_{\begin{subarray}{c}\\{\gamma_{e},\Gamma_{e}^{(2)},\gamma_{b},\Gamma_{e,b}^{(3/2)}\end{subarray}\\}\rightarrow\Psi}\left\\{\left(T+V\right)[\gamma_{e}]+\left(W+H_{d}\right)[\Gamma_{e}^{(2)}]+H_{ph}[\gamma_{b}]+H_{I}[\Gamma_{e,b}^{(3/2)}]\right\\}.$
(3.52)

So in principle, we could replace the variation over all wave functions $\Psi$
by their respective set of RDMs needed to define the energy expectation
values. Instead of varying over the full configuration space
$(\mathbf{x}_{1},...,\mathbf{x}_{N},$ $p_{1},...,p_{M})$, the above
reformulation seems to indicate that we can replace this by varying over
$(\mathbf{r},\mathbf{r}^{\prime})$ for the diagonal of $\Gamma_{e}^{(2)}$ and
also for the 1RDM $\gamma_{e}$ together with a variation over $(\alpha,\beta)$
for $\gamma_{b}$ and over $(\alpha,\mathbf{r})$ for $\Gamma^{(3/2)}_{e,b}$.
Such a reformulation is the basis of any RDM theory, and for electronic
systems the properties of RDMs have been studied for more than 50 years [195].
As we have seen in Sec.2.4 for the electronic case, this seeming reduction of
complexity is deceptive. The many-body problem has merely been shifted to the
representability conditions of the RDMs. Although the conditions are simple
for 1RDMs of ensembles, for any higher order RDM, they are extremely
complicated.

Let us transfer the representability question to the present case. In order to
find physically sensible results, we cannot vary arbitrarily over the above
RDMs but need to ensure that they are consistent among each other and that
they are all connected to _the same_ physical wave function. This is indicated
in Eq. (3.3.1), where
$\\{\gamma_{e},\Gamma_{e}^{(2)},\gamma_{b},\Gamma_{e,b}^{(3/2)}\\}\rightarrow\Psi$
highlights that the RDMs are contractions of a common wave function. Thus, we
are confronted with two new difficulties in the coupled case: consistency
between the RDMs and representability with respect to one wave function. In
the purely electronic problem, all the occurring RDMs are connected and thus,
we can calculate the energy expectation value only by one single object, the
2RDM. Consequently, we just need to bother about the representability problem
$\Gamma_{e}^{(2)}\rightarrow\Psi$. For the coupled case, we still have the
connection between $\gamma_{e}$ and $\Gamma_{e}^{(2)}$, but we cannot connect
them to the other RDMs. There is also no kind of coupled higher order RDM that
by different contraction yields all the others as we discussed before. Thus,
we need to treat all the three objects
$\\{\Gamma_{e}^{(2)},\gamma_{b},\Gamma_{e,b}^{(3/2)}\\}$ together and find
representability conditions for this set. The insights from electronic theory
might be helpful here, but they are not generalizable in a straightforward
way. One of the crucial new problems is the underlying Fock-space for the
photon-degrees of freedom. As the name suggests, $N$-representability
conditions are always defined with respect to a constant particle number.
These considerations show yet from a different perspective how intricate a
many-body description of the coupled-electron photon problem is, and this even
in the simplest dipole-approximated case.

#####  3.3.2 QED-RDMFT

From the analysis of the last subsection, we gained yet another point of view
on the difficulties of the electron-photon interaction. Transferring concepts
from electronic-structure theory is also from an RDM perspective far from
trivial and provides us with many interesting research questions. In this last
part of the section, let us see how far we can get with a “conservative”
approach that is to generalize electronic RDMFT to its QED version in the same
way as Ruggenthaler [242] generalized DFT. To the best of our knowledge, the
proof of the according Hohenberg-Kohn like theorem has not been published in
the literature and thus, we present it in the following

We will show the proof for the full Pauli-Fierz theory with minimal coupling
and only afterwards show its form for the cavity setting of Hamiltonian (3.1).
We follow Gilbert’s proof [199] of a one-to-one mapping between the von-
Neumann density matrix of the system and the “internal” pair of
$(\gamma(\mathbf{r},\mathbf{r}^{\prime}),\mathbf{A}(\mathbf{r}))$, where
$\gamma$ is the electronic 1RDM and $\mathbf{A}$ the vector potential of the
coupled electron-photon ground state. The corresponding “external” pair is
$(v(\mathbf{r},\mathbf{r}^{\prime}),\mathbf{j}(\mathbf{r}))$, where v is a
possibly non-local external potential and $\mathbf{j}$ an external charge
current. As in electronic RDMFT, there is no unique $v$ corresponding to a
given $\gamma$. For the sake of generality, we will present the proof in SI
units.

###### Basic Definitions

We consider non-relativistic QED in the setting that we introduced in Sec.
1.3.1, i.e., we neglect any form of spin-dependence of the Hamiltonian (e.g.,
due to the Stern-Gerlach term) and quantize the theory in the Coulomb gauge
where for the vector potential holds $\nabla\cdot\mathbf{A}=0$. Thus, the
corresponding vector-potential operator is purely transversal.202020For
notational convenience, we use in this subsection the symbol $\mathbf{A}$ to
denote the transversal vector potential. We then expand

$\displaystyle\hat{\mathbf{A}}(\mathbf{r},t)=\sqrt{\tfrac{\hbar
c^{2}}{\epsilon_{0}L^{3}}}\sum_{\mathbf{k},\sigma}\tfrac{\bm{\epsilon}_{\mathbf{k},\sigma}}{\sqrt{2\omega_{q}}}\left(\hat{a}_{\mathbf{k},\sigma}e^{\mathbf{i}\mathbf{k}\cdot\mathbf{r}}+\hat{a}_{\mathbf{k},\sigma}^{+}e^{-\mathbf{i}\mathbf{k}\cdot\mathbf{r}}\right),$
(3.53)

where we assume a quantization box with side-length $L$, leading to allowed
wave vectors $\mathbf{k}=2\pi\mathbf{n}/L$ ($\mathbf{n}\in\mathbb{Z}^{3}$) and
corresponding frequencies $\omega_{\mathbf{k}}=c|\mathbf{k}|$. The index
$\sigma=\\{1,2\\}$ denotes the transversal polarization directions and the
corresponding polarization vectors obey
$\epsilon_{\mathbf{k},\sigma}\cdot\epsilon_{\mathbf{k},\sigma}=\delta_{\sigma,\sigma^{\prime}}$
and $\epsilon_{\mathbf{k},\sigma}\cdot\mathbf{k}=0$. The expansion
coefficients become the usual creation (annihilation) operators
$\hat{a}_{\mathbf{k},\sigma}^{(\dagger)}$ with (transversal) commutation
relations
$[\hat{a}_{\mathbf{k},\sigma},\hat{a}_{\mathbf{k},\sigma}^{\dagger}]=\delta_{\sigma,\sigma^{\prime}}\delta^{T}(\mathbf{k}-\mathbf{k}^{\prime})$,
where $\delta^{T}$ is the transversal delta-distribution. We denoted the
vacuum permittivity with $\epsilon_{0}$, the speed of light with $c$ and the
Planck constant with $\hbar$. We consider a classical external charge current
$\mathbf{j}$ that because of the transversality of $\mathbf{A}$ enters also
only by its transversal component. We expand $\mathbf{j}$ in the modes of the
quantization box, i.e.,

$\displaystyle\mathbf{j}(\mathbf{r})=\sqrt{\tfrac{\epsilon_{0}\hbar}{L^{3}}}\sum_{\mathbf{k},\sigma}\omega_{\mathbf{k}}\tfrac{\epsilon_{\mathbf{k},\sigma}}{\sqrt{2\omega_{q}}}\left(j_{\mathbf{k},\sigma}e^{\mathbf{i}\mathbf{k}\cdot\mathbf{r}}+j^{*}_{\mathbf{k},\sigma}e^{-\mathbf{i}\mathbf{k}\cdot\mathbf{r}}\right),$
(3.54)

with expansion coefficients
$\mathbf{j}_{\mathbf{k},\sigma}=\mathbf{j}_{-\mathbf{k},\sigma}^{*}=(2\omega_{\mathbf{k}}^{3}\epsilon_{0}\hbar
L^{3})^{-1/2}\int\mathrm{d}^{3}r\epsilon_{\mathbf{k},\sigma}\cdot\mathbf{j}(\mathbf{r})\exp(-\mathbf{i}\mathbf{k}\cdot\mathbf{r})$.
We couple this field to $N$ non-relativistic electrons. The Hamiltonian of the
full coupled system reads then (cf. Def. 2)

$\displaystyle\hat{H}=$
$\displaystyle\sum_{i=1}^{N}\left\\{\frac{1}{2m}\left(-\mathbf{i}\hbar\nabla_{i}+\frac{e}{c}\hat{\mathbf{A}}(\mathbf{r}_{i})\right)^{2}-ev(\mathbf{r}_{i};\mathbf{r}_{i}^{\prime})\right\\}+\frac{1}{2}\sum_{i,j=1}^{N}\frac{e^{2}}{|\mathbf{r}_{i}-\mathbf{r}_{k}|}$
$\displaystyle+\sum_{\mathbf{k},\sigma}\hbar\omega_{\mathbf{k}}(\hat{a}_{\mathbf{k},\sigma}^{\dagger}\hat{a}_{\mathbf{k},\sigma}+\tfrac{1}{2})-\sum_{\mathbf{k},\sigma}\hbar\omega_{\mathbf{k}}(\hat{a}_{\mathbf{k},\sigma}j^{*}_{\mathbf{k},\sigma}+\hat{a}_{\mathbf{k},\sigma}^{\dagger}j_{\mathbf{k},\sigma}),$
(3.55)

where $m$ is the electron mass, and $e$ the elementary charge. We introduced
now a general non-local external potential $v(\mathbf{r};\mathbf{r}^{\prime})$
that acts on a function $f(\mathbf{r})$ in an integral sense as
$\hat{v}f(\mathbf{r})=\int\mathrm{d}^{3}r^{\prime}v(\mathbf{r};\mathbf{r}^{\prime})f(\mathbf{r}^{\prime})$.
Note that for the time-independent case, introducing an external vector
potential is equivalent to an external current [242]. The physical charge
current that is preserved in a time-independent setting (because it obeys the
continuity equation) reads

$\displaystyle\hat{\mathbf{J}}(\mathbf{r})=\hat{\mathbf{j}}_{P}(\mathbf{r})-\frac{2}{mc}\hat{n}(\mathbf{r})\hat{\mathbf{A}}(\mathbf{r}),$
(3.56)

where

$\displaystyle\hat{\mathbf{j}}_{P}(\mathbf{r})=\frac{\mathbf{i}e\hbar}{2m}\sum_{i=1}^{N}(\delta(\mathbf{r}-\mathbf{r}_{k})\overset{\rightarrow}{\nabla}_{k}-\overset{\leftarrow}{\nabla}_{k}\delta(\mathbf{r}-\mathbf{r}_{k}))$
(3.57)

is the paramagnetic current and
$\hat{n}(\mathbf{r})=-e\sum_{i=1}^{N}\delta(\mathbf{r}-\mathbf{r}_{k})$ is the
charge density operator. Note that by introducing the term
$\frac{e^{2}}{|\mathbf{r}_{i}-\mathbf{r}_{k}|}$, we assumed already the
$L\rightarrow\infty$ limit.

The Hamiltonian (3.3.2) is a hermitian operator on the Hilbert space

$\displaystyle\mathfrak{H}=$
$\displaystyle\mathfrak{H}_{A}^{N}\otimes\mathfrak{F},$

where $\mathfrak{H}_{a}^{N}$ is the antisymmetric Hilbert space of $N$
electrons and $\mathfrak{F}$ is the photon Fock space (see Sec. 3.1.1). We
denote the eigenstates of $\hat{H}$ with $\ket{\Psi_{i}}$, i.e., they satisfy
the static Schrödinger equation $\hat{H}\ket{\Psi_{i}}=E_{i}\ket{\Psi_{i}}$.

The $\ket{\Psi_{i}}$ can be parametrized by $N$ 4-dimensional spin-spatial
coordinates $\mathbf{x}=(\mathbf{r},\sigma)$ and some parametrization of the
photon degrees of freedom that we denote by $\alpha=(\mathbf{k},\sigma)$. We
choose this representation in analogy to the real-space representation of the
electrons. Note however that this description is only well-defined in the
mode-space and the “back-transformation” to real space is
problematic.212121The concept of a real-space photon wave function is the
object of a long-standing debate and a good summary of this debate can be
found in [89, Chap. 1.5.4]. For the concrete example of a one-photon-wave
function in real-space and its connection to mode-space, the reader is
referred to, e.g., Ref. [255]. The author considers there the spontaneous
emission process of an atom and shows that the photon wave function diverges
at the position of the atom. However, for all other positions, it is well-
defined and makes physically sense. This is a good example of the problem:
there are several ways to define a photon wave function in real-space that are
reasonable, but they are all limited in a similar way as the example of Ref.
[255]. We then define

$\displaystyle\Psi(\mathbf{r}_{1},...,\mathbf{r}_{N},\alpha_{1},\alpha_{2},...)=\braket{\mathbf{r}_{1},...,\mathbf{r}_{N},\alpha_{1},\alpha_{2},...}{\Psi},$
(3.58)

where we leave the exact number of photons arbitrary, because it is not
determined by the Hamiltonian. More generally, we will not only consider pure
states, i.e., the eigenstates $\ket{\Psi_{i}}$ of $\hat{H}$, but ensembles of
such eigenstates that are described by the density matrix

$\displaystyle
D(\mathbf{r}_{1},...,\mathbf{r}_{N},\alpha_{1},\alpha_{2},...;\mathbf{r}_{1}^{\prime},...,\mathbf{r}_{N}^{\prime},\alpha_{1}^{\prime},\alpha_{2}^{\prime},...)=\sum_{i}w_{i}\Psi_{i}(\mathbf{r}_{1},...,\mathbf{r}_{N},\alpha_{1},\alpha_{2},...)\Psi^{*}_{i}(\mathbf{r}_{1}^{\prime},...,\mathbf{r}_{N}^{\prime},\alpha_{1}^{\prime},\alpha_{2}^{\prime},...),$
(3.59)

where $0\leq w_{i}\leq 1$, and $\sum_{i}w_{i}=1$, making the set
$P^{N}_{\infty}$ of all such $D$ convex (see Sec. 2.4.1). In the following, we
are interested in the electronic one-body reduced density matrix (1RDM)

$\displaystyle\gamma(\mathbf{r},\mathbf{r}^{\prime})=$ $\displaystyle
N\int\mathrm{d}\mathbf{x}_{2}\cdots\mathrm{d}\mathbf{x}_{N}\mathrm{d}\mathbf{q}_{1}\mathrm{d}\mathbf{q}_{2}\cdots\,D(\mathbf{r}_{1},...,\mathbf{r}_{N},\alpha_{1},\alpha_{2},...;\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}...,\mathbf{r}_{N},\alpha_{1},\alpha_{2},...)$
$\displaystyle=$ $\displaystyle\gamma[D],$ (3.60)

and the expectation value of the vector potential operator from Eq. (3.53),
i.e.,

$\displaystyle\mathbf{A}(\mathbf{r})=\operatorname{tr}[D\mathbf{A}].$ (3.61)

###### Gilbert’s theorem for NR-QED

With these definitions, we are now set to show the following correspondence

$\displaystyle(v,\mathbf{j})\longrightarrow
D_{0}\overset{1:1}{\longleftrightarrow}(\gamma_{0},\mathbf{A}_{0}).$ (3.62)

We will do this in two steps. First, we show that the pair $(v,\mathbf{j})$
determines the (non-degenerate) ground state density-matrix
$D_{0}=\Psi_{0}\Psi_{0}^{*}$ (though the opposite is not true in general) and
then, we will show the one-to-one correspondence between $D_{0}$ and the pair
$(\gamma_{0},\mathbf{A}_{0})$, where $\gamma_{0}=\gamma[D_{0}]$ and
$\mathbf{A}_{0}=\operatorname{Tr}[D_{0}\mathbf{A}]$. This proof is a direct
generalization of Gilbert’s proof [199] and establishes QED-RDMFT.

The first part is to show that we can associate to every pair $(v,j)$ a
corresponding ground state density matrix $D_{0}$. This is however simple,
because the pair $(v,\mathbf{j})$ consists of the only external parameters in
the Hamiltonian and thus we can label the Hamiltonian as

$\displaystyle H=H(v,\mathbf{j}).$

Obviously, we can also label the (unique) ground state of $H(v,\mathbf{j})$ in
the same way, i.e., $D_{0}=D_{0}(v,\mathbf{j})$. However, the opposite is not
true, i.e., we cannot label $v(D)$, because there are $v\neq v^{\prime}$ that
have the same ground state density matrix $D_{0}=D_{0}^{\prime}$. Take for
instance a non-interacting $N$-electron system that consists only of one-body
terms. The Hamiltonian
$\hat{H}=\sum_{i=1}^{N}=-1/2\nabla_{i}^{2}+v(\mathbf{r}_{i};\mathbf{r}_{i}^{\prime})$
of such a system per construction commutes with
$\gamma=\sum_{i}n_{i}\phi_{i}^{*}(\mathbf{r}^{\prime})\phi_{i}(\mathbf{r})$
and thus the ground state is simply the Slater determinant of the lowest $N$
natural orbitals $\phi_{1},...,\phi_{N}$. Consider now the nonlocal potential
$v+\delta v$, with $\delta
v(\mathbf{r};\mathbf{r}^{\prime})=\sum_{ij=N+1}^{\infty}v_{ij}\phi_{i}^{*}(\mathbf{r}^{\prime})\phi_{j}(\mathbf{r})$
that only depends on the unoccupied natural orbitals. Obviously, $v+\delta v$
is different from $v$, but both potentials lead to the same ground state,
since they differ only in the unoccupied subspace.222222Note that conversely,
such a construction is not possible, if _all_ natural orbitals have a non-zero
occupation number. In this case, the mapping is unique. It has been
conjectured that this is the case for Coulomb systems, but there is still no
general proof [256]. Thus, there are infinitely many non-interacting systems
with different nonlocal potential that have the same ground state and thus,
one cannot construct a KS-like auxiliary system in (zero-temperature)
RDMFT.232323 If we however generalize the description to grand-canonical
ensembles, the mapping between $\gamma$ and $v$ is indeed unique for all
systems [200]. In DFT, there is instead a unique KS system for every density,
since we consider only local potentials that cannot separate between the
occupied and unoccupied space [165]. Thus there is a one-to-one correspondence
between such potentials and the density. However, to define RDMFT (or QED-
RDMFT), we do not need this unique mapping as Gilbert has pointed out [199].

Let us now prove the second part of the theorem, i.e., the one-to-one
correspondence between ground state density matrix $D_{0}$ and the internal
variables $(\gamma_{0},\mathbf{A}_{0})$. We do so by reductio ad absurdum and
show that the opposite assumption leads to a contradiction. Thus, we
considering two Hamiltonians $H(v^{1},\mathbf{j}^{1}),H(v^{2},\mathbf{j}^{2})$
that have to two _different_ ground states $D_{0}^{1},D_{0}^{2}$ but the same
$(\gamma_{0}^{1},\mathbf{A}_{0}^{1})=(\gamma_{0}^{2},\mathbf{A}_{0}^{2})=(\gamma_{0},\mathbf{A}_{0})$.
We have then

$\displaystyle E_{0}^{1}=\text{Tr}[D_{0}^{1}H(v_{1},\mathbf{j}_{1})]<$
$\displaystyle\text{Tr}[D_{0}^{2}H(v^{1},\mathbf{j}^{1})]$ $\displaystyle=$
$\displaystyle
E_{0}^{2}+\int\mathrm{d}\mathbf{r}\mathrm{d}\mathbf{r}^{\prime}\gamma^{2}(\mathbf{r},\mathbf{r}^{\prime})\left(v^{1}(\mathbf{r},\mathbf{r}^{\prime})-v^{2}(\mathbf{r},\mathbf{r}^{\prime})\right)-\int\mathrm{d}\mathbf{r}\mathbf{A}^{2}(\mathbf{r})\left(\mathbf{j}^{1}(\mathbf{r})-\mathbf{j}^{2}(\mathbf{r})\right)$
$\displaystyle E_{0}^{2}=\text{Tr}[D_{0}^{2}H(v^{2},\mathbf{j}^{2})]<$
$\displaystyle E^{2}=\text{Tr}[D_{0}^{1}H(v^{2},\mathbf{j}^{2})]$
$\displaystyle=$ $\displaystyle
E_{0}^{1}+\int\mathrm{d}\mathbf{r}\mathrm{d}\mathbf{r}^{\prime}\gamma^{1}(\mathbf{r},\mathbf{r}^{\prime})\left(v^{2}(\mathbf{r},\mathbf{r}^{\prime})-v^{1}(\mathbf{r},\mathbf{r}^{\prime})\right)-\int\mathrm{d}\mathbf{r}\mathbf{A}^{1}(\mathbf{r})\left(\mathbf{j}^{2}(\mathbf{r})-\mathbf{j}^{1}(\mathbf{r})\right),$

Since $(\gamma^{1},\mathbf{j}^{1})=(\gamma^{2},\mathbf{j}^{2})$, the sum of
both inequalities leads to

$\displaystyle E_{1}+E_{2}<E_{1}+E_{2},$

which is obviously a contradiction. Thus, we have proven the assumption.

We can trivially transfer this proof to the dipole-Hamiltonian (3.1) including
an external current, cf. Eq. (3.30), that we have employed also for the QEDFT
mapping. Therefore, we have

$\displaystyle(v(\mathbf{r},\mathbf{r}^{\prime}),\\{p_{\alpha}\\})\longrightarrow
D_{0}\overset{1:1}{\longleftrightarrow}(\gamma(\mathbf{r},\mathbf{r}^{\prime}),\\{\dot{j}_{\alpha}\\})$
(3.63)

###### QED-RDMFT in practice

We want to conclude the subsection with a small discussion of the
applicability of QED-RDMFT. We employ our basic setting, i.e., Hamiltonian
(3.1) and consider again wave functions, since the ground state of a closed
system is pure and does not require the ensemble construction.

The proof of the last paragraph is the justification to employ QED-RDMFT in
practice. The first step for this is the definition of the energy functional
to identify the corresponding unknown parts, especially the photon exchange
correlation functional. We can directly conclude from our discussion in Sec.
3.3 that we will have to find approximations for

$\displaystyle E_{XC,ph}[\gamma,\\{\dot{j}_{\alpha}\\}]=$
$\displaystyle\braket{\Psi[\gamma,\\{\dot{j}_{\alpha}\\}]}{\left[\hat{H}_{int}+\hat{H}_{self,2}\right]\Psi[\gamma,\\{\dot{j}_{\alpha}\\}]},$
(3.64)

where

$\displaystyle\hat{H}_{int}=$
$\displaystyle-\sum_{i}^{N}\sum_{\alpha}^{M_{ph}}\omega_{\alpha}p_{\alpha}\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{i}$
(3.65) $\displaystyle\hat{H}_{self,2}=$ $\displaystyle\frac{1}{2}\sum_{i\neq
j=1}^{N_{e}}\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{i}\right)\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{j}\right)$
(3.66)

Since both these terms are local in space, the 1RDM perspective does not
directly provide an advantage in comparison to the alternative QEDFT
description. However, typical RDMFT functionals are more generic than DFT
functionals as we have discussed in Sec. 2.4.2. We can for example generalize
straightforwardly the Müller functional to describe the self-interaction term,
simply by exchanging
$w(\mathbf{r},\mathbf{r}^{\prime})\rightarrow\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}\right)\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}^{\prime}\right)$.
The reason is that the Müller functional (as many other RDMFT functionals) is
basically an approximation of the 2RDM in terms of the 1RDM.

Unfortunately, since we do not yet have a good understanding of the $3/2$-body
RDM (Eq. (3.50)), that is connected to the interaction term, we cannot use a
similar “trick” for the first term in Eq. (3.64). However, in contrast to
QEDFT, we have some concrete tools that we can try to apply to this problem.
For instance, this suggests to investigate the still unknown representability
conditions of coupled RDMs (see the Outlook in part IV).

## PART IIDressed Orbitals - Old Theory in a New Basis

> You can look at something with a microscope and see it a certain way, you
> can look at it with a naked eye and see it in a certain way, you look at it
> with a telescope and you see it in another way. Now, which level of
> magnification is the correct one? Well, obviously, they’re all correct, but
> they’re just different points of view.  (Alan Watts, 1960 [257])

We now have discussed at length the challenges of accurately describing
equilibrium many-body states of non-relativistic QED. Even if we consider only
one photon mode, a wave function description of the many-body space is already
infeasible for very small systems. Since the electronic many-body problem is a
very well-studied topic with a plethora of theories and optimized
implementations that are capable to treat many systems extremely accurately,
it would be very desirable to extend such methods to the coupled electron-
photon problem. However, we have shown in Sec. 3.1.2 with the example of
Hartree-Fock theory that such an extension is not as straightforward as one
might hope. A direct generalization of the Hartree-Fock (single-reference)
wave function ansatz is not capable to describe any quantum effects of the
photon field. When we try to go beyond that, even a seemingly very restrictive
ansatz leads to an (over-)exponentially growing problem.

In order to study to role of the quantum nature of the electron-photon
interaction, we cannot avoid dealing with the _multi-reference character_ of
the coupled states. This is an important difference to electronic problems
where a single-reference ansatz, i.e., one Slater determinant already included
quantum-mechanical exchange because of its anti-symmetry. Understanding the
role of such quantum effects of the photon field is especially important for
polaritonic ground states, because static classical electric and magnetic
fields are zero for many equilibrium settings. In the time-dependent case
instead, the classical description is very accurate, when it is combined with
a quantum mechanical description of the electrons. We made this observation,
when we were discussing KS-QEDFT (Sec. 3.2). QEDFT is in principle capable to
treat quantum fluctuations of the photon-field by means of corresponding
exchange-correlation potentials. However, it is very difficult to construct an
accurate functional and the only known functional that goes beyond the MF, the
“photon OEP,” can besides being numerically very challenging only account for
one-photon processes. For instance, we cannot apply the photon OEP to larger
systems and one-photon processes are not sufficient to accurately describe
polaritonic ground states in the strong-coupling regime [247]. Consequently,
we need new methods to appropriately study such settings.

Let us therefore take one step back and analyze again the origin of the multi-
reference character of polaritonic systems. We describe electrons and photons
as _separate species_ in their own antisymmetric and symmetric many-body
spaces. These are then coupled by the interaction operator, which introduces
most of the challenges of the coupled problem (see Ch. 3). There is a clear
analogy to correlated electronic systems, which we also describe starting from
non-interacting electrons that then are coupled by the interaction. However,
this coupling takes place _within_ one antisymmtric many-body space and this
is an important difference to the electron-photon coupling. The central idea
of the dressed(-orbital) construction is now to change the basic description
of coupled electron-photon systems such that there is only one (effective)
many-body space with one (effective) symmetry. This means specifically that we
define a single-particle space with already coupled electron-photon orbitals,
i.e., _polaritonic_ or _dressed_ orbitals and construct the many-body space as
a product of such orbital spaces. We can then describe non-interacting
polaritons that are coupled due to the corresponding dressed interaction. Such
a programme naturally defines polaritons as its own particle class with its
own statistics that needs to be respected to construct the many-polariton
space.

Clearly, restructuring the many-electron-photon space in such a way is not
straightforward. In fact, the dressed construction considers an auxiliary
system that is even higher-dimensional than the cavity-QED Hamiltonian.
However, this extension provides the necessary flexibility to introduce an
exact reformulation of the cavity-QED systems in terms of dressed orbitals.
Most importantly, this allows to approach the challenges of describing coupled
light-matter systems from a different perspective. Approximations to the wave
function, the density description and reduced density matrices can be defined
also in the dressed setting. This is yet another way to utilize these state
descriptors and they will have again different properties with respect to the
already discussed versions of the previous chapter. A practical consequence of
this different characterization of light-matter states is that it allows for
relatively simple approximation schemes. For instance, the HF approximation to
the many-polariton wave function explicitly accounts for correlation between
the electronic and photonic subsystems in the original description.

A key role is hereby played by the _new structure_ of the auxiliary
Hamiltonian with respect to the standard cavity-QED Hamiltonian. Expressed in
terms of polaritonic orbitals, the new Hamiltonian consists only of one-body
and two-body operators. It is thus structurally equivalent to the electronic-
structure Hamiltonian (Sec. 2.1). Consequently, we can generalize in principle
every electronic-structure method to a “polaritonic-structure method” to
describe cavity-QED problems. For instance, this allows to apply RDM methods
to coupled light-matter systems, without the need to deal with the little-
understood $3/2$-body RDM. Also functional construction for a polariton-based
QEDFT is simpler, since one can more directly apply the successful strategies
of electronic-structure theory.

Another important advantage of the structural resemblance between polaritonic
and electronic-structure methods is the numerical implementation. Since these
methods always require to solve nonlinear equations that do not allow for an
analytical treatment, we need numerical solvers to apply any method in
practice. By employing the polariton description, we can simply use the
existing electronic-structure codes as a basis and extend them accordingly. We
do not have to develop entirely new codes.

However, to be able to introduce polaritonic orbitals, we have to increase
their dimension by one with every photon mode that we take into account. For
instance, to describe the standard case of cavity QED that considers one
effective photon mode, the polaritonic orbitals depend on three electronic and
one photonic coordinate and are thus four-dimensional. Considering such
orbitals (and depending on the system size also with more modes) is still
numerically feasible for modern high-performance clusters. Thus, polaritonic-
structure methods will be especially useful in such cavity-QED settings.

A further important feature of the polariton description is the exchange-
symmetry and the according statistics of polaritonic orbitals, which have a
_fermion-boson hybrid_ character. Such a hybrid statistics is not usual in
quantum many-body theory, which leads to interesting new research questions at
a very basic level. In fact, we also had underestimated its role when we
started the investigation of dressed orbitals. Only after we had obtained
reproducible numerical results that violated the Pauli principle, we
understood the significance of the hybrid statistics (see part III). On the
one hand, this leads to new challenges especially with regard to enforcing the
statistics in practice. On the other hand, the hybrid statistics are a new
tool to describe many-particle spaces of coupled species in general, which
could be very valuable also for other coupled problems.

### Chapter 4 The dressed-orbital construction

In this chapter, we introduce the dressed-orbital construction that allows to
describe coupled electron-photon systems in terms of new orbitals. This is the
basis for all the following discussions. The construction requires several
steps, which are nontrivial and we try to explain each of them in a way, as
simple and slow as possible. Specifically, we start in Sec. 4.1 with a brief
motivation for a polariton-based description of cavity-QED problems. Then, we
introduce the construction with an example (Sec. 4.2), that we generalize in
Sec. 4.3. We note that the dressed construction has been discussed in three
publications, the first of which was written by my colleagues Nielsen et al.
[258] and in the other two I was the main author ([1, 2]). This chapter is
based on Ref. [2], which is the most recent publication and thus presents the
most complete picture of the construction.

####  4.1 Polaritons and correlation

We want to start the discussion with some general comments about polaritons
and the multi-reference character of coupled electron-photon systems. We have
discussed in the introduction (Sec. 1.1) that the hallmark of strong electron-
photon coupling is the emergence of light-matter hybrid states or polaritons.
The basic physics of such hybrid states is described by a minimal model (see
Sec. 1.2) that considers two relevant electronic states, labelled by
$\ket{m_{1}}$ (“ground”) and $\ket{m_{2}}$ (“excited” state), and two photonic
states, the vacuum $\ket{0}$ and one-photon state $\ket{1}$. The resulting
polaritons are then denoted as

$\displaystyle
P_{\pm}=\alpha\ket{m_{1}}\otimes\ket{1}\pm\beta\ket{m_{2}}\otimes\ket{0},$
(4.1)

with coefficients $\alpha,\beta$ that depend on details of the model (see Sec.
1.2 for further details about the family of cavity-QED models). From a first-
principles perspective this means that polaritons correspond to _multi-
reference_ or correlated wave functions, even if we described the matter
states within a single-reference approach. We have discussed this in detail in
Sec. 3.1.2. For strongly-coupled ground states which are the focus of this
thesis, the results of more general cavity-QED models [24] predict that an
accurate description requires even more terms, i.e., more references, than the
two of the “minimal” polariton (4.1). Systems that require multi-reference
states are well known in electronic-structure theory and their accurate
description is among the most difficult challenges in the field.111See Ch. 2,
especially the last paragraph of Sec. 2.3 on strong correlation. A classical
example is the dissociation of the Hydrogen molecule (Sec. 2.5), where the
wave function obtains a multi-reference character for larger bond distances.
This prototypical system is commonly used as a challenging test case for
first-principle methods [259, 221, 260, 261].

However, the term multi-reference depends crucially on the basic entities,
i.e., the “single references” that are used to build the “multi-reference”
state. Transferred to the problem of polaritonic physics, we face the
fundamental problem that separate matter and photon wave functions 222We want
to remind the reader on the subtleties regarding the definition of a photon
wave function that we have discussed in Sec. 3.1.1 (as in the simple example
above) can be very inefficient in describing polaritonic states of the coupled
system. We saw this for example in Sec. 3.1.2, when we were discussing a
possible generalized HF wave function. Even a very simplified ansatz brings
back the exponential wall. We also saw that the (single-reference) KS system
in equilibrium QEDFT reduces the description basically to an effective matter-
only system. All effects of the photon field (besides some trivial static
shift for matter systems with a permanent dipole) have to be carried by the
unknown exchange-correlation potential and there is no known approximation
that is accurate in the strong-coupling regime [247].

One important challenge for the equilibrium description is that the exchange-
correlation functional only depends on matter quantities, since the photon
displacement coordinate is trivial (cf. Eq. (3.37) of Ch. 3). From a general
perspective, this is rooted in the choice of a non-interacting auxiliary
system, where the electron and photon spaces decouple. Thus, we are confronted
with a dilemma: On the one hand, strongly-coupled light matter systems require
fundamentally a multi-reference description, because polaritons, i.e., hybrid
matter-photon particles emerge as the principal degree of freedom. On the
other hand, a multi-reference description is especially difficult for coupled
systems, because the product space of the separate electron and photon spaces
is considerably more difficult to describe than the single-species spaces.
Especially, typical (and straightforward) approximations such as an effective
single-reference description capture considerably less effects than the
equivalent approaches in a matter-only theory.

In this chapter, we want to propose a construction that tries to mediate
between the opposing sides of this dilemma. The key idea is to build the
theory not as usual from separate electron and photon states but somehow
introduce electron-photon hybrid states as basic entities. Before we lay out
our construction to introduce such a “many-polariton” theory, we want to
illustrate the challenges that arise, when we try to modify the standard
description. Let us therefore shortly recapitulate how many-body spaces are
constructed with the example of the electronic problem. We start with the
single-particle Hilbert-space $\tilde{\mathfrak{h}}^{1}$, that is spanned by a
corresponding single-particle basis $\\{\tilde{\psi}_{i}(\mathbf{r})\\}$
(where we always assume bases to be orthonormalized, i.e.,
$\int\mathrm{d}\mathbf{r}\tilde{\psi}_{i}^{*}(\mathbf{r})\tilde{\psi}_{j}(\mathbf{r})=\delta_{ij}$),
such that every single-particle state can be expressed by some superposition

$\displaystyle\Psi(\mathbf{r})=\sum_{i=1}^{\infty}c_{i}\tilde{\psi}_{i}(\mathbf{r}),$

where the $c_{i}$ are expansion coefficients that need to satisfy the sum rule
$\sum_{i}|c_{i}|^{2}=1$. In the case of electrons, we need to include also the
spin-degree of freedom, which we do by a tensor-product. The spin-space
$\mathfrak{S}$ for a one-electron problem is just two-dimensional and we
denote its basis-elements by $\\{\alpha(\sigma),\beta(\sigma)\\}$. The tensor-
product $\mathfrak{h}=\tilde{\mathfrak{h}}\otimes\mathcal{S}$ then combines
both sets to one larger basis, taking into account all possible combinations.
We subsume both together in a spin-spatial coordinate
$\mathbf{r}\rightarrow\mathbf{x}=(\mathbf{r},\sigma)$ where $\sigma$ is the
spin-coordinate. The expansion then becomes

$\displaystyle\Psi(\mathbf{x})=$
$\displaystyle\sum_{j}c_{j}\psi_{j}(\mathbf{x})$ $\displaystyle=$
$\displaystyle\sum_{s=\alpha,\beta}\sum_{i=1}^{\infty}c_{i,s}\psi_{i,s}(\mathbf{r},\sigma)$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{\infty}\left[c_{i,\alpha}\tilde{\psi}_{i}(\mathbf{r})\otimes\alpha(\sigma)+c_{i,\beta}\tilde{\psi}_{i}(\mathbf{r})\otimes\beta(\sigma)\right].$

When we now want to describe two particles, we actually go exactly the same
way as we have done for the spin and consider the tensor product of the
single-particle spaces

$\displaystyle\mathfrak{h}^{2}=\mathfrak{h}^{1}\otimes\mathfrak{h}^{1}.$

This two-particle Hilbert space is however still too general, because it
includes _distinguishable_ states (Sec. 2.1). Since electrons are fermions, we
have to constrain $\mathfrak{h}^{2}$ to only antisymmetric combinations, i.e.,

$\displaystyle\mathfrak{h}_{A}^{2}=\mathfrak{h}\wedge\mathfrak{h}.$

We can make this explicit in the state parametrization by employing Slater
determinants $\Psi_{i,j}(\mathbf{x}_{1},\mathbf{x}_{2})=$
$\frac{1}{\sqrt{2}}\left(\psi_{i}(\mathbf{x}_{1})\psi_{j}(\mathbf{x}_{2})-\psi_{j}(\mathbf{x}_{1})\psi_{i}(\mathbf{x}_{2})\right)=|\psi_{i}\,\psi_{j}|_{-}$
as basis of the two-particle space. If we label every combination of i and j
by a new index $I=(i,j)$, we can describe a general state
$\Psi\in\mathfrak{h}_{A}^{2}$ as

$\displaystyle\Psi(\mathbf{x}_{1},\mathbf{x}_{2})=\sum_{I}c_{I}\Psi_{I}(\mathbf{x}_{1},\mathbf{x}_{2}),$

where $\sum_{I}|c_{I}|^{2}=1$. The advantage of this construction is obvious:
the basis already takes care about the fundamental particle symmetry. We have
discussed at length that only this additional symmetry information includes
nontrivial quantum effects in the description.

When we want to describe coupled electron-photon problems, we lose this
advantage. Electrons and photons are two different species and thus
distinguishable. Accordingly, we have to consider the full tensor product
between the matter and photon many-body spaces without further symmetry
restrictions (Ch. 3.1.1). Let us contrast this to a hypothetical many-
polariton space. In analogy to before, we start with the single-polariton
Hilbert-space that is build by objects that have an electron coordinate
$\mathbf{x}$ and a photon coordinate $\mathbf{p}$. Similarly to the extension
of the spatial coordinates by a spin-component for electrons, we could define
$\mathfrak{h}_{p}=\mathfrak{h}_{e}\otimes\mathfrak{h}_{ph}$, where
$\mathfrak{h}_{e}$ is the electronic one-particle space from before and
$\mathfrak{h}_{ph}$ is some one-particle photon space. Introducing a photon
orbital basis $\\{\chi_{\alpha}(\mathbf{p})\\}$, the elements of
$\mathfrak{h}_{p}$ are thus given by

$\displaystyle\phi(\mathbf{x},\mathbf{p})=\sum_{i,\alpha}c_{i\alpha}\left(\psi_{i}(\mathbf{x})\otimes\chi_{\alpha}(\mathbf{p})\right),$
(4.2)

where $\sum_{i\alpha}|c_{i\alpha}|^{2}=1$. Since electrons and photons are
distinguishable, there is no further symmetrization required.

We have discussed in Sec. 2.1.2 that many electronic-structure problems can be
simplified, if we assume that the orbitals are the eigenfunctions of the one-
body part of the Hamiltonian,
($h(\mathbf{r})=-\frac{1}{2}\nabla_{\mathbf{r}}^{2}+v(\mathbf{r})$). To get a
feeling for the challenges of these new kind of orbitals (4.2), let us try to
do find a similar “one-body” Hamiltonian. The straightforward choice is to
consider the cavity-QED Hamiltonian (Eq. (1.37)) for one electron, i.e.,

$\displaystyle\hat{h}_{p}=\frac{1}{2}\nabla^{2}+v(\mathbf{r})+\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].$

The eigenfunctions $\phi(\mathbf{x},\mathbf{p})$ of $\hat{h}_{p}$ would define
our photon variable as M-dimensional, i.e., $\mathbf{p}=(p_{1},...,p_{M})$.
This dimension is determined by the system that we want to describe and as we
have discussed, we can often assume even $M=1$ for a cavity system.

The next step is to construct the many-body space from this. According to the
standard prescription, we consider

$\displaystyle\mathfrak{h}=\mathfrak{h}_{p}\otimes\mathfrak{h}_{p},$

with elements

$\displaystyle\Phi(\mathbf{x}_{1}\mathbf{p}_{1},\mathbf{x}_{2}\mathbf{p}_{2})=\sum_{I,J}c_{I,J}\phi_{I}(\mathbf{x}_{1}\mathbf{p}_{1})\phi_{J}(\mathbf{x}_{2}\mathbf{p}_{2}).$

Now we need to understand how we utilize $\Phi$ to describe the ground state
of the actual Hamiltonian that describes two electrons, coupled to $M$ modes,
i.e.,

$\displaystyle\hat{H}=\sum_{i=1}^{2}\frac{1}{2}\nabla_{\mathbf{r}_{i}}^{2}+v(\mathbf{r}_{i})+\frac{1}{2}\left[-\partial_{p}^{2}+\omega^{2}\left(p+\frac{\bm{\lambda}}{\omega}\cdot(\mathbf{r}_{1}+\mathbf{r}_{2})\right)^{2}\right].$
(cf. 3.14)

This confronts us with two fundamental questions:

  1. 1.

How can we enforce the correct symmetry on this level?

  2. 2.

How can we deal with the additional $\mathbf{p}_{2}$-coordinates that are not
present Eq. (3.14)?

In the following, we show how the dressed-orbital construction can in
principle resolve both of the above questions. This allows to construct a
many-polariton space (almost) of the form

$\displaystyle\mathfrak{h}_{p}^{N}=\mathcal{S}_{p}(\underbrace{\mathfrak{h}_{p}\otimes\cdots\otimes\mathfrak{h}_{p}}_{N\text{times}}),$

where $\mathcal{S}_{p}$ enforces the underlying exchange-symmetry that is of a
fermion-boson hybrid nature.

####  4.2 Polaritons in the dressed auxiliary system

Before we discuss the general case, we want to illustrate the dressed
construction with the example of the 2-electron-1-mode system that we have
considered already in Sec. 3.1.1.333Note that the contents of this section are
part of Ref. [2] Since the construction requires many (sometimes quite
technical) steps, we reserve this whole section only for the example system.
This allows us to go through every step in detail and motivate its purpose. In
the next section (Sec. 4.3), we then recapitulate all these steps for the
general case of $N$ electrons coupled to $M$ photon modes. The basic idea of
the construction is sketched in Fig. 4.1.

Figure 4.1: Sketch of the auxiliary construction for an example wave function
$\Psi(x,y,p)$ of two one-dimensional electrons $(x,y)$ coupled to one photon
mode with displacement coordinate p. The coupling is indicted by the double
arrows $x\leftrightarrow p,y\leftrightarrow p$. The electronic orbital wave
functions are symbolized by the ground and first excited state of a box with
zero-boundary conditions and the photon mode by a wiggly line. The
corresponding dressed wave function $\Psi^{\prime}(x,q_{1},y,q_{2})$ has
instead two photon-coordinates $q_{1},q_{2}$ that are related to the physical
coordinate p by $p=1/\sqrt{2}(q_{1}+q_{2})$. On the wave function level, this
connection can be utilized to introduce two polariton orbitals with
coordinates $(x,q_{1})$ and $(y,q_{2})$, respectively, that are interacting.
This new interaction is indicated by a double arrow between the two orbitals.

#####  4.2.1 The auxiliary Hamiltonian: Matching electronic and photonic
coordinates

The Hamiltonian of the 2-electron-1-mode system (introduced in Sec. 3.1.1)
reads

$\displaystyle\hat{H}=\sum_{k=1}^{2}\frac{1}{2}\Delta_{\mathbf{r}_{k}}+v(\mathbf{r}_{k})+\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}+\frac{1}{2}\left[-\partial_{p}^{2}+\omega^{2}\left(p+\frac{\bm{\lambda}}{\omega}\cdot(\mathbf{r}_{1}+\mathbf{r}_{2})\right)^{2}\right].$
(cf. 3.14)

We start by expanding $\hat{H}$ in its single constituents

$\displaystyle\begin{split}\hat{H}=&\overbrace{\sum_{k=1}^{2}\left[-\frac{1}{2}\Delta_{\mathbf{r}_{k}}+v(\mathbf{r}_{k})+\tfrac{1}{2}(\lambda\cdot\mathbf{r}_{k})^{2}\right]}^{\hat{H}^{e}=\sum_{k=1}^{2}h^{e}(\mathbf{r}_{k})}\overbrace{-\sum_{k=1}^{2}\omega
p\bm{\lambda}\\!\cdot\\!\mathbf{r}_{k}}^{\hat{H}^{ep}=\sum_{k=1}^{2}h^{ep}(\mathbf{r}_{k},p)}\\\
&+\underbrace{\tfrac{1}{2}\left(-\frac{\partial^{2}}{\partial
p^{2}}+\omega^{2}p^{2}\right)}_{\hat{H}^{ph}=h^{ph}(p)}+\underbrace{\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}+\tfrac{1}{2}\left(\bm{\lambda}\\!\cdot\\!\mathbf{r}_{1}\right)\left(\bm{\lambda}\\!\cdot\\!\mathbf{r}_{2}\right)}_{\hat{H}^{ee}=h^{ee}(\mathbf{r}_{1},\mathbf{r}_{2})},\end{split}$
(4.3)

where we grouped terms according to the coordinates that are appearing.444Note
that we employ in this section the symbol $\Delta\equiv\nabla^{2}$ for the
Laplacian. We observe that there is an asymmetry between electron- and photon-
coordinates: there are two electron coordinates that appear in a completely
symmetric way in all the terms. However there is only one mode-displacement
coordinate, independently on the number of electrons.

The basics idea of the dressed construction is now to simply introduce as many
_artificial_ photon coordinates as electron coordinates, which formally
“lifts” the asymmetry. In this case, this means to introduce one further
photon coordinate $p_{2}$.555For a system of $N$ electrons, we need to
introduce therefore $N-1$ extra photon coordinates (for every mode). This has
to be done such that the original system cannot be influenced. We can imagine
another photon cavity that is very far away from our system, but that we want
to describe at the same time. If we further assume that this artificial mode
has the same frequency as the physical one, the auxiliary system is described
by the Hamiltonian (we denote the quantities of the auxiliary system by a
prime-symbol)

$\displaystyle\hat{H}^{\prime}=$ $\displaystyle\hat{H}+h^{ph}(p_{2})$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{2}h^{e}(\mathbf{r}_{k})+h^{ee}(\mathbf{r}_{1},\mathbf{r}_{2})-\sum_{k=1}^{2}\omega
p\bm{\lambda}\\!\cdot\\!\mathbf{r}_{k}+h^{ph}(p)+h^{ph}(p_{2}).$ (4.4)

If we denote by $\chi(p_{2})$ an eigenstate of the additional photon
Hamiltonian $h^{ph}(p_{2})$ and $\Psi(\mathbf{x}_{1},\mathbf{x}_{2},p)$ is an
eigenstate of $\hat{H}$, then

$\displaystyle\Psi^{\prime}(\mathbf{x}_{1},\mathbf{x}_{2},p,p_{2})=\Psi(\mathbf{x}_{1},\mathbf{x}_{2},p)\chi(p_{2}),$
(4.5)

is an eigenstate of $\hat{H}^{\prime}$. The ground state of $\hat{H}^{\prime}$
is then uniquely defined and given by $\Psi^{\prime}=\Psi\chi_{0}$ with the
harmonic oscillator vacuum state $\chi_{0}$.

With $\hat{H}^{\prime}$, we almost achieved our goal to symmetrize electron
and photon coordinates, but there is still the electron-photon interaction
term $h^{ep}(\mathbf{r}_{k},p)$ that depends only on $p$, but not on $p_{2}$.
The basic idea to remedy this stems from the similarity of the Hamiltonian
(4.2.1) with a molecular center-of-mass Hamiltonian [262]. Transferred to this
case, we interpret $p$ as the “center-of-displacement” coordinate of two
independent displacements $q_{1},q_{2}$. The second key ingredient of the
dressed construction is therefore a coordinate transformation
$(p,p_{2})\rightarrow(q_{1},q_{2})$ that sets

$\displaystyle p\propto(q_{1}+q_{2}).$ (4.6)

Obviously, the transformation should be norm-conserving to not increase the
displacement artificially, but more importantly, it needs to keep the photon-
energy part form-invariant, i.e.,

$\displaystyle h^{ph}(p)+h^{ph}(p_{2})\rightarrow
h^{ph}(q_{1})+h^{ph}(q_{2}).$ (4.7)

Since
$h^{ph}(p)=\tfrac{1}{2}\left(-\frac{\mathrm{d}^{2}}{\mathrm{d}p}+\omega^{2}p^{2}\right)=\left[\tfrac{1}{\sqrt{2}}(\mathbf{i}\frac{\mathrm{d}}{\mathrm{d}p},\omega
p)\right]\cdot\left[\tfrac{1}{\sqrt{2}}(\mathbf{i}\frac{\mathrm{d}}{\mathrm{d}p},\omega
p)\right]$ has the form of an inner product, we can simply employ _an
orthogonal transformation_ for that. Specifically, a transformation for a set
of variables $(x_{1},x_{2},...)\rightarrow(x_{1}^{\prime},x_{2}^{\prime},...)$
is orthogonal, if it leaves any expression of the form
$\sum_{i}x_{i}^{2}\rightarrow\sum_{i}x_{i}^{\prime}{}^{2}$ invariant. In the
two-particle case, there is exactly one possible choice for that, which reads
(see also Fig. 4.1)

$\displaystyle\begin{split}p=&1/\sqrt{2}(q_{1}+q_{2})\\\
p_{2}=&1/\sqrt{2}(q_{1}-q_{2}).\end{split}$ (4.8)

As constructed, the transformation changes only
$\hat{H}^{ep}=\sum_{k=1}^{2}h^{ep}(\mathbf{r}_{k},p)$, which becomes

$\displaystyle\hat{H}^{ep}\rightarrow$
$\displaystyle-\sum_{k=1}^{2}\omega/\sqrt{2}(q_{1}+q_{2})\bm{\lambda}\\!\cdot\\!\mathbf{r}_{k}$
$\displaystyle=$
$\displaystyle-\sum_{k=1}^{2}\omega/\sqrt{2}q_{k}\bm{\lambda}\\!\cdot\\!\mathbf{r}_{k}-\omega/\sqrt{2}(q_{1}\bm{\lambda}\\!\cdot\\!\mathbf{r}_{2}+q_{2}\bm{\lambda}\\!\cdot\\!\mathbf{r}_{1}).$
(4.9)

Again, we grouped the terms according to the particle indices, which we can do
now also for the photon coordinates. Since the auxiliary system after the
transformation is symmetric under the exchange of photon indices, we can
arbitrarily exchange the indices of the $q$-coordinates. Putting all terms
together, we get the following new form of the auxiliary Hamiltonian

$\displaystyle\hat{H}^{\prime}=$
$\displaystyle\sum_{k=1}^{2}\left[-\tfrac{1}{2}\left(\Delta_{\mathbf{r}_{k}}+\tfrac{\mathrm{d}^{2}}{\mathrm{d}q_{k}}\right)+v(\mathbf{r}_{k})+\tfrac{\omega^{2}}{2}q_{k}^{2}-\tfrac{\omega}{\sqrt{2}}q_{k}\bm{\lambda}\\!\cdot\\!\mathbf{r}_{k}+\tfrac{1}{2}(\lambda\cdot\mathbf{r}_{k})^{2}\right]$
$\displaystyle+\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}-\tfrac{\omega}{\sqrt{2}}(q_{1}\bm{\lambda}\\!\cdot\\!\mathbf{r}_{2}+q_{2}\bm{\lambda}\\!\cdot\\!\mathbf{r}_{1})+\tfrac{1}{2}\left(\bm{\lambda}\\!\cdot\\!\mathbf{r}_{1}\right)\left(\bm{\lambda}\\!\cdot\\!\mathbf{r}_{2}\right).$
(4.10)

Interestingly, this Hamiltonian has the same _structure_ as the electronic
many-body Hamiltonian 2.1. We can see this more explicitly, when we introduce
the new _dressed_ coordinate

$\displaystyle\mathbf{z}=(\mathbf{r},q)$ (4.11)

and the according dressed Laplacian666We do not need further differential
operators here, but of course we can generalize every differential operator to
the new coodinates. The dressed nabla-operator would read for example
$\nabla_{\mathbf{z}}=(\nabla_{\mathbf{r}},\frac{\mathrm{d}}{\mathrm{d}q})$.

$\displaystyle\Delta_{\mathbf{z}}{}^{\prime}=\Delta_{\mathbf{r}}+\frac{\mathrm{d}^{2}}{\mathrm{d}q}.$
(4.12)

We can write then

$\displaystyle\hat{H}^{\prime}=$
$\displaystyle\sum_{k=1}^{2}\left[-\tfrac{1}{2}\Delta_{\mathbf{z}_{k}}^{\prime}+v^{\prime}(\mathbf{z}_{k})\right]+\tfrac{1}{2}\sum_{k\neq
l=1}^{2}w^{\prime}(\mathbf{z}_{k},\mathbf{z}_{l}),$ (4.13) where the dressed
local potential reads $\displaystyle v^{\prime}(\mathbf{z})=$ $\displaystyle
v(\mathbf{r})+\tfrac{\omega^{2}}{2}q_{k}^{2}-\tfrac{\omega}{\sqrt{2}}q_{k}\bm{\lambda}\\!\cdot\\!\mathbf{r}_{k}+\tfrac{1}{2}(\lambda\cdot\mathbf{r}_{k})^{2}$
(4.14) and the dressed interaction kernel is $\displaystyle
w^{\prime}(\mathbf{z}_{k},\mathbf{z}_{l})=$
$\displaystyle\tfrac{1}{|\mathbf{r}_{k}-\mathbf{r}_{l}|}-\tfrac{\omega}{\sqrt{2}}(q_{k}\bm{\lambda}\\!\cdot\\!\mathbf{r}_{l}+q_{l}\bm{\lambda}\\!\cdot\\!\mathbf{r}_{k})+\tfrac{1}{2}\left(\bm{\lambda}\\!\cdot\\!\mathbf{r}_{k}\right)\left(\bm{\lambda}\\!\cdot\\!\mathbf{r}_{l}\right).$
(4.15)

#####  4.2.2 The auxiliary wave function from a first glance

With this construction, we have found a way to introduce new coordinates and
rewrite our basic Hamiltonian in a form that is completely symmetric with
respect to the new coordinates. The next step is to define the corresponding
wave function. Let us therefore recall the Pauli citation of Sec. 2.1.3:

> _If we deal with many similar (indistinguishable, A/N) particles, special
> circumstances occur that arise from the Hamilton operator being invariant
> under any permutations of particles._ (Pauli, 1933 [140])

According to Pauli, exchange-symmetry is a consequence of the symmetry of the
Hamiltonian, which indicates that we are quite close to our goal to formulate
a many-polariton theory. We only need to derive the details of the “special
circumstances,” i.e., the _type_ of exchange symmetry that the auxiliary wave
function. In the standard description, there are only two choices: symmetric
and antisymmetric wave functions. However, we will see that the situation here
is a bit more involved.

###### Polaritonic symmetry

Let us briefly recapitulate the steps to arrive at the dressed Hamiltonian
(4.2.1). We start with the physical Hamiltonian (Eq. (4.3)) that depends on
the photon coordinate $p$. We then add another harmonic oscillator term,
depending on the auxiliary coordinate $p_{2}$ which is still distinguishable
from the original $p$ (because of $\hat{h}^{ep}$) and obtain Eq. (4.2.1).
Finally, the coordinate transformation (4.8) introduces new photon coordinates
$p_{1},p_{2}$, which enter the transformed Hamiltonian (4.2.1) in a completely
symmetric way. Since the Hamiltonian is also symmetric in the electronic
coordinates, we can deduce the polaritonic symmetry, i.e., the Hamiltonian
(4.2.1) is invariant under the exchange of the two polaritonic coordinates
$\mathbf{z}_{1}\leftrightarrow\mathbf{z}_{2}$.

With that we have defined all symmetries on the Hamiltonian level. Regarding
the wave function
$\Psi^{\prime}(\mathbf{z}_{1}\sigma_{1},\mathbf{z}_{2}\sigma_{2})$, we still
have to understand whether it is symmetric or antisymmetric under
$\mathbf{z}_{1}\sigma_{1}\leftrightarrow\mathbf{z}_{2}\sigma_{2}$. For that,
we first note that the photon ground state is inversion-symmetric, i.e.,

$\displaystyle\chi(p_{2})=\chi(-p_{2}).$ (4.16)

Further, we want to preserve the antisymmetry of the physical wave function in
the electronic coordinates, i.e.,
$\Psi(\mathbf{x}_{1},\mathbf{x}_{2},p)=-\Psi(\mathbf{x}_{2},\mathbf{x}_{1},p)$.
Taken both together, we deduce that

$\displaystyle\Psi^{\prime}(\mathbf{x}_{1},q_{1},\mathbf{x}_{2},p,q_{2})=$
$\displaystyle\hskip
9.0pt\Psi(\mathbf{x}_{1},\mathbf{x}_{2},1/\sqrt{2}(q_{1}+q_{2}))\chi(1/\sqrt{2}(q_{1}-q_{2}))$
$\displaystyle\overset{\mathbf{z}_{1}\sigma_{1}\leftrightarrow\mathbf{z}_{2}\sigma_{2}}{=}$
$\displaystyle-\Psi(\mathbf{x}_{2},\mathbf{x}_{1},1/\sqrt{2}(q_{2}+q_{1}))\chi(1/\sqrt{2}(q_{2}-q_{1}))$
$\displaystyle\overset{Eq.\leavevmode\nobreak\
\eqref{eq:dressed_construction_auxiliary_wf_requirement}}{=}$
$\displaystyle-\Psi(\mathbf{x}_{2},\mathbf{x}_{1},1/\sqrt{2}(q_{1}+q_{2}))\chi(1/\sqrt{2}(q_{1}-q_{2}))$
$\displaystyle=$
$\displaystyle-\Psi^{\prime}(\mathbf{x}_{2},q_{2}\mathbf{x}_{1},q_{1}).$
(4.17)

The dressed wave function $\Psi^{\prime}$ is _antisymmetric_ with respect to
the exchange of dressed coordinates
$(\mathbf{x},q)=(\mathbf{r},\sigma,q)=(\mathbf{z},\sigma)$.

Crucially, this allows us to describe the coupled electron-photon Hilbert
space in a different many-body basis, i.e., in terms of Slater determinants of
dressed orbitals (see also Fig. 4.1). Therefore, let us go back to our
polariton-orbital basis $\\{\phi(\mathbf{z},\sigma)\\}$ of Sec. 4.1. We can
now explicitly follow the steps of the electronic theory: we consider a basis
$\phi_{i}(\mathbf{z},\sigma)$ of $\mathfrak{h}_{p}$. For instance, this could
be the eigenstates of the single-particles Hamiltonian
$h(\mathbf{z})=-\tfrac{1}{2}\Delta_{\mathbf{z}}{}^{\prime}+v^{\prime}(\mathbf{z})$.
The two-particle space is then given as the antisymmetrized tensor-product

$\displaystyle\mathfrak{h}_{p}^{2}=\mathcal{A}(\mathfrak{h}_{p}\otimes\mathfrak{h}_{p})$
(4.18)

and the elements of $\mathfrak{h}_{p}^{2}$ can be parametrized as

$\displaystyle\begin{split}\Psi^{\prime}(\mathbf{z}_{1},\sigma_{1},\mathbf{z}_{2},\sigma_{2})=&\sum_{ij}c_{ij}\left[\vphantom{\sum}\phi_{i}(\mathbf{z}_{1},\sigma_{1})\phi_{j}(\mathbf{z}_{2},\sigma_{2})-\phi_{j}(\mathbf{z}_{1},\sigma_{1})\phi_{i}(\mathbf{z}_{2},\sigma_{2})\right]\\\
\overset{!}{=}&\sum_{i,j,\alpha,\beta}C^{\prime}{}_{ij}^{\alpha,\beta}\psi_{i}(\mathbf{x}_{1})\psi_{j}(\mathbf{x}_{2})\chi^{\alpha}(q_{1})\chi^{\beta}(q_{2}).\end{split}$
(4.19)

It is clear that since the auxiliary configuration space is much larger than
the original configuration space, we made an (beyond simple systems)
infeasible numerical problem even more infeasible. What the dressed
construction offers us in compensation, is a _new structure_ of the
configuration space. In terms of the above expansion, we hope that the new
coefficients $c_{ij}$ are easier to approximate than the original
$C_{ij}^{\alpha}$. The main advantage is that the auxiliary wave function can
be expanded in terms of Slater determinants and we have a lot of methods at
our disposal that are geared toward such a situation, i.e., first-principles
electronic structure theories. Yet the dressed formulation allows for
relatively simple approximation schemes, e.g., HF theory in terms of a single
polaritonic Slater determinant (see Sec. 2). Importantly, such simple wave
functions in terms of polaritonic orbitals correspond to correlated (multi-
determinant) wave functions in physical space. In Sec. 6, we present
calculations based on dressed orbitals which agree remarkably well with exact
references even for very strong coupling strengths. We also illustrate there
the explicitly correlated character of the electronic subsystem.

###### The explicit form of the auxiliary wave function

Let us now have a second closer look at the wave function. In fact, we can
explicitly construct the dressed wave function in the new coordinates

$\displaystyle\Psi^{\prime}(\mathbf{x}_{1},\mathbf{x}_{2},p,p_{2})=\Psi(\mathbf{x}_{1},\mathbf{x}_{2},p)\chi(p_{2})\rightarrow\Psi(\mathbf{x}_{1},\mathbf{x}_{2},1/\sqrt{2}(q_{1}+q_{2}))\chi(1/\sqrt{2}(q_{1}-q_{2}))$
(4.20)

if we know $\Psi$. To see this, we consider some electronic orbital-basis
$\\{\psi_{i}(\mathbf{x})\\}$ and the eigenfunctions $\\{\chi^{\alpha}(p)\\}$
of the photon Hamiltonian $\hat{h}^{ph}$ and expand

$\displaystyle\Psi(\mathbf{x}_{1},\mathbf{x}_{2},p)=\sum_{i,j,\alpha}C_{ij}^{\alpha}\psi_{i}(\mathbf{x}_{1})\psi_{j}(\mathbf{x}_{2})\chi^{\alpha}(p).$

The first step of the construction of the auxiliary ground state is simple. We
(tensor-)multiply $\Psi$ with the ground state $\chi^{0}(p_{2})$ of the extra
Harmonic oscillator that by construction is also an eigenstate of
$\hat{h}^{ph}$, i.e.,

$\displaystyle\Psi^{\prime}(\mathbf{x}_{1},\mathbf{x}_{2},p,p_{2})=\left[\sum_{i,j,\alpha}C_{ij}^{\alpha}\psi_{i}(\mathbf{x}_{1})\psi_{j}(\mathbf{x}_{2})\chi^{\alpha}(p)\right]\chi^{0}(p_{2}).$
(4.21)

The nontrivial step is the transformation (4.8) to the new coordinates.
However, since we know the analytical expression of the photon basis,

$\displaystyle\chi^{\alpha}(p)=\frac{1}{\sqrt{2^{\alpha}\alpha!}}\left(\frac{\omega}{\pi}\right)^{1/4}e^{-\omega
p^{2}/2}H_{\alpha}(\sqrt{\omega}p),$ (4.22)

where
$H_{\alpha}(z)=(-1)^{\alpha}e^{z^{2}}\frac{\mathrm{d}^{\alpha}}{\mathrm{d}z^{\alpha}}\left(e^{-z^{2}}\right)$
are Hermite-polynomials, we can perform the transformation (4.8) _explicitly_.
Specifically, we have to calculate terms of the form

$\displaystyle\chi^{\alpha}(1/\sqrt{2}(q_{1}+q_{2}))\chi^{0}(1/\sqrt{2}(q_{1}-q_{2}))\propto$
$\displaystyle
e^{-\omega(q_{1}+q_{2})^{2}/4}H_{\alpha}(\sqrt{\omega/2}(q_{1}+q_{2}))H_{0}e^{-\omega(q_{1}-q_{2})^{2}/4}$
$\displaystyle=$ $\displaystyle
e^{-\omega(q_{1}+q_{2})^{2}/4}H_{\alpha}(\sqrt{\omega/2}(q_{1}+q_{2}))e^{-\omega(q_{1}-q_{2})^{2}/4}$

for all $\alpha$, where we used that $H_{0}=1$. We start with the Gaussian
part of the oscillator states and calculate explicitly

$\displaystyle e^{-\omega p^{2}/2}e^{-\omega p_{2}^{2}/2}\rightarrow$
$\displaystyle e^{-\omega(q_{1}+q_{2})^{2}/4}e^{-\omega(q_{1}-q_{2})^{2}/4}$
$\displaystyle=$ $\displaystyle e^{-\omega q_{1}^{2}/2}e^{-\omega
q_{2}^{2}/2}.$

The transformation merely exchanges the coordinates, because it is orthogonal
(cf. Eq. (4.7)).777Thus, also for more than two particles, the Gaussian part
of the states remains form-invariant under the transformation.

For the remaining part involving the Hermite-polynomial, the transformation is
more involved, but we can calculate it still analytically. Making use of the
identity

$\displaystyle
H_{\alpha}(z_{1}+z_{2})=2^{-\alpha/2}\sum_{\beta=0}^{\alpha}\binom{\alpha}{\beta}H_{\alpha-\beta}(z_{1}\sqrt{2})H_{\beta}(z_{2}\sqrt{2}),$
(4.23)

with $z_{i}=\sqrt{\omega/2}q_{i}$, we arrive after some algebra at

$\displaystyle\Psi^{\prime}(\mathbf{x}_{1},\mathbf{x}_{2},q_{1},q_{2})=$
$\displaystyle\sum_{i,j,\alpha}C_{ij}^{\alpha}\psi_{i}(\mathbf{x}_{1})\psi_{j}(\mathbf{x}_{2})\sum_{\beta=0}^{\alpha}\sqrt{\frac{\alpha!}{\beta!(\alpha-\beta)!}}\frac{1}{\sqrt{2^{\alpha}}}\chi^{(\alpha-\beta)}(q_{1})\chi^{\beta}(q_{2})$
$\displaystyle=$
$\displaystyle\sum_{i,j,\alpha,\beta}C^{\prime}{}_{ij}^{\alpha,\beta}\psi_{i}(\mathbf{x}_{1})\psi_{j}(\mathbf{x}_{2})\chi^{\alpha}(q_{1})\chi^{\beta}(q_{2}).$
(4.24)

In the last line, we introduced the polaritonic expansion coefficient
$C^{\prime}{}_{ij}^{\alpha,\beta}$, which is uniquely determined by the above
equation.

First of all, Eq. (4.2.2) highlights how the value of $p$ is “distributed”
over the new coordinates $q_{1}$ and $q_{2}$. This shows how 2 (or in general
$N$) polaritonic orbitals can carry collectively the information of only one
mode ($M$ modes). Beyond that, this explicit construction illustrates how
intricate the coordinate transformation acts on the system. It is true that
Eq. (4.2.2) is antisymmetric under
$\mathbf{z}_{1}\sigma_{1}\leftrightarrow\mathbf{z}_{2}\sigma_{2}$ and thus it
is covered by our just derived ansatz
$\Psi^{\prime}(\mathbf{z}_{1}\sigma_{1},\mathbf{z}_{2}\sigma_{2})$. However,
the opposite is clearly not true, i.e., not all polaritonic wave functions
that are antisymmetric will also be of the form (4.2.2). This is a first
indication that the ansatz (4.19) might be limited. Nevertheless, there is one
generic feature in Eq. (4.2.2): it is symmetric with respect to the new
coordinates $q_{1},q_{2}$.888Note that we can in principle derive a similar
(but considerably more cumbersome) expression for $N$ electrons and $M$ modes.
To derive Eq. (4.2.2), we only have to make use of the orthogonality and the
explicit transformation of $p$. The auxiliary coordinates $p_{2},p_{3},...$ do
not play a role because only for $p$, we have to take Hermite-polynomials
$H_{\alpha}$ with $\alpha>0$ into account. We discuss this in the next
subsection.

#####  4.2.3 The auxiliary wave function on a second glance

In fact, the antisymmetric polariton space $\mathfrak{h}_{p}^{2}$ that we
constructed in the previous subsection, includes states that do not have a
correspondence in the physical system. Let us illustrate this with our two-
electron-one-mode example, neglecting no electron-electron and electron-photon
interactions. The corresponding _independent particles_ (IP) Hamiltonian reads

$\displaystyle\hat{H}_{IP}=\sum_{k=1}^{2}-\tfrac{1}{2}\Delta_{\mathbf{r}_{k}}+v(\mathbf{r}_{k})\tfrac{1}{2}\left(-\frac{\mathrm{d}^{2}}{\mathrm{d}p}+\omega^{2}p^{2}\right).$
(4.25)

The ground state of $\hat{H}_{IP}$ is simply

$\displaystyle\Psi(\mathbf{r}_{1}\sigma_{1},\mathbf{r}_{2}\sigma_{2},p)=\frac{1}{\sqrt{2}}\left[\vphantom{\sum}\psi_{1}(\mathbf{r}_{1}\sigma_{1})\psi_{2}(\mathbf{r}_{2}\sigma_{2})-\psi_{2}(\mathbf{r}_{1}\sigma_{1})\psi_{1}(\mathbf{r}_{2}\sigma_{2})\right]\chi_{0}(p),$

a product of a Slater determinant consisting of the two lowest eigenfunctions
$\psi_{1},\psi_{2}$ of $[-\tfrac{1}{2}\Delta_{\mathbf{r}}+v(\mathbf{r})]$ and
$\chi_{0}(p)=(\frac{\omega}{\pi})^{1/4}e^{-\omega p^{2}/4}$, which is the
ground state of a harmonic oscillator with frequency $\omega$. As shown in the
last subsection, we obtain the dressed version of $\Psi$ by multiplication
with another oscillator ground state $\chi_{0}(p_{2})$ with the same
frequency. Performing the coordinate transformation (4.8) in this case is
especially simple because also$\chi_{0}(p)$ is a ground state. Consequently,
the transformation does nothing else than replacing the coordinates
$(p,p_{2})$ with $(q_{1},q_{2})$. The complete auxiliary ground state reads
then

$\displaystyle\Psi^{\prime}(\mathbf{z}_{1}\sigma_{1},\mathbf{z}_{2}\sigma_{2})$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\vphantom{\sum}\psi_{1}(\mathbf{r}_{1}\sigma_{1})\psi_{2}(\mathbf{r}_{2}\sigma_{2})-\psi_{2}(\mathbf{r}_{1}\sigma_{1})\psi_{1}(\mathbf{r}_{2}\sigma_{2})\right]\chi_{0}(q_{1})\chi_{0}(q_{2})$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\vphantom{\sum}\psi_{1}(\mathbf{r}_{1}\sigma_{1})\chi_{0}(q_{1})\psi_{2}(\mathbf{r}_{2}\sigma_{2})\chi_{0}(q_{2})-\psi_{2}(\mathbf{r}_{1}\sigma_{1})\chi_{0}(q_{1})\psi_{1}(\mathbf{r}_{2}\sigma_{2})\chi_{0}(q_{2})\right]$
$\displaystyle=\frac{1}{\sqrt{2}}[\phi_{10}(\mathbf{z}_{1}\sigma_{1})\phi_{20}(\mathbf{z}_{2}\sigma_{2})-\phi_{20}(\mathbf{z}_{1}\sigma_{1})\phi_{10}(\mathbf{z}_{2}\sigma_{2})],$

where in the last line we subsumed the electronic and photonic orbitals with
the same coordinate index to a dressed orbital, i.e.,
$\psi_{i}(\mathbf{r}_{j}\sigma_{j})\chi_{n}(q_{j})\equiv\phi_{in}(\mathbf{z}_{j}\sigma_{j})$.
This wave function is obviously antisymmetric with respect to the exchange of
$\mathbf{z}_{1}\sigma_{1}$ and $\mathbf{z}_{2}\sigma_{2}$. Now let us consider
another wave function,

$\displaystyle\tilde{\Psi}^{\prime}(\mathbf{z}_{1}\sigma_{1},\mathbf{z}_{2}\sigma_{2})$
$\displaystyle=\phi_{10}(\mathbf{z}_{1},\sigma_{1})\phi_{11}(\mathbf{z}_{2},\sigma_{2})-\phi_{11}(\mathbf{z}_{1},\sigma_{1})\phi_{10}(\mathbf{z}_{2},\sigma_{2})$
$\displaystyle=\frac{1}{\sqrt{2}}\left[\vphantom{\sum}\psi_{1}(\mathbf{r}_{1}\sigma_{1})\chi_{0}(q_{1})\psi_{1}(\mathbf{r}_{2}\sigma_{2})\chi_{1}(q_{2})-\psi_{1}(\mathbf{r}_{1}\sigma_{1})\chi_{1}(q_{1})\psi_{1}(\mathbf{r}_{2}\sigma_{2})\chi_{0}(q_{2})\right]$
$\displaystyle=\psi_{1}(\mathbf{r}_{1}\sigma_{1})\psi_{1}(\mathbf{r}_{2}\sigma_{2})\left[\vphantom{\sum}\chi_{0}(q_{1})\chi_{1}(q_{2})-\chi_{1}(q_{1})\chi_{0}(q_{2})\right],$

that is a dressed Slater determinant and thus part of $\mathfrak{h}_{p}^{2}$.
This state is however unphysical, because is _violates the Pauli-principle_ :
both electrons occupy the same electronic orbital. Depending on their energy
eigenvalues
$E[\Psi^{\prime}]=\braket{\Psi^{\prime}}{\hat{H}_{IP}^{2}{}^{\prime}\Psi^{\prime}}$,
a minimization of the dressed Hamiltonian without ensuring the hybrid
statistics could determine either of the two wave functions as the ground
state. Only if
$E[\Psi^{\prime}]=\epsilon_{1}+\epsilon_{2}+\omega<2\epsilon_{1}+2\omega=E[\tilde{\Psi}^{\prime}]$,
where $\epsilon_{1}$ and $\epsilon_{2}$ are the eigenenergies corresponding to
$\psi_{1}$ and $\psi_{2}$, a simple minimization would yield the correct
solution. If, however, $\epsilon_{2}-\epsilon_{1}>\omega$, then a minimization
of the dressed problem would lead to the state $\tilde{\Psi}^{\prime}$ that
violates the Pauli principle. Within the dressed Slater determinant, the
antisymmetry can be “carried” by the electronic or the photonic part of the
orbital, which allows for states like $\tilde{\Psi}^{\prime}$. For an
interacting problem, both cases cannot be separated so easily, but the problem
remains in principle the same, as we show numerically in Sec. 6. Note that
this is very similar to the violations of the $N$-representability conditions
in variational 2RDM theory (see Sec. 2.4.1). This means that we either have to
make sure that $\omega$ is large compared to the electronic excitations, such
that the unrestricted minimization with _only_ fermionic symmetry in
$\mathbf{z}\sigma$ picks the right wave function [1] (fermion ansatz), or we
have to enforce the hybrid statistics.

To guarantee the Pauli principle, we actually have to build the dressed many-
body space by requiring antisymmetry only with respect to the electronic
coordinates together with symmetry in terms of the photonic coordinates. We
can summarize this as

$\displaystyle\mathbf{r}_{1}\sigma_{1}$
$\displaystyle\leftrightarrow\mathbf{r}_{2}\sigma_{2}\quad$
$\displaystyle\rightarrow\quad\Psi^{\prime}$
$\displaystyle\leftrightarrow-\Psi^{\prime}$ (4.26a) $\displaystyle\bm{q}_{1}$
$\displaystyle\leftrightarrow\bm{q}_{2}$
$\displaystyle\rightarrow\quad\Psi^{\prime}$
$\displaystyle\leftrightarrow\Psi^{\prime}.$ (4.26b)

This does not mean that we have to discard the wave function expansion in
terms of dressed Slater determinants (4.19), because from (4.26) follows also

$\displaystyle\mathbf{z}_{1}\sigma_{1}\leftrightarrow\mathbf{z}_{2}\sigma_{2}\quad$
$\displaystyle\rightarrow\quad\Psi^{\prime}\leftrightarrow-\Psi^{\prime}.$
(4.27)

If we require (4.27), which builds $\mathfrak{h}_{p}^{2}$ and then constrain
the space by _either_ (4.26a) or (4.26b), we fulfil (4.26) equivalently. In
Sec. 4.3, we show for the general case of $N$ electrons and $M$ modes that
(4.26) is indeed necessary and sufficient to guarantee a _one-to-one
correspondence_ between the physical and auxiliary ground state. This
establishes the polariton description of cavity QED.

Let us now illustrate how the extra condition indeed rules our the unphysical
state $\tilde{\Psi}^{\prime}$. We can enforce the constitutive relations on
$\tilde{\Psi}^{\prime}$ by adding two extra terms that exchange either
electronic or photonic coordinates, which leads to

$\displaystyle\tilde{\Psi}^{\prime}_{hybrid}(\mathbf{z}_{1}\sigma_{1},\mathbf{z}_{2}\sigma_{2}){}={}$
$\displaystyle\overbrace{\phi_{10}(\mathbf{r}_{1}q_{1}\sigma_{1})\phi_{11}(\mathbf{r}_{2}q_{2}\sigma_{2})-\phi_{11}(\mathbf{r}_{1}q_{1}\sigma_{1})\phi_{10}(\mathbf{r}_{2}q_{2}\sigma_{2})}^{\tilde{\Psi}^{\prime}}$
$\displaystyle+\phi_{10}(\mathbf{r}_{1}q_{2}\sigma_{1})\phi_{11}(\mathbf{r}_{2}q_{1}\sigma_{2})-\phi_{11}(\mathbf{r}_{1}q_{2}\sigma_{1})\phi_{10}(\mathbf{r}_{2}q_{1}\sigma_{2})$
$\displaystyle{}={}$
$\displaystyle\psi_{1}(\mathbf{r}_{1}\sigma_{1})\psi_{1}(\mathbf{r}_{2}\sigma_{2})\left[\vphantom{\sum}\chi_{0}(q_{1})\chi_{1}(q_{2})-\chi_{1}(q_{1})\chi_{0}(q_{2})\right]$
$\displaystyle+\psi_{1}(\mathbf{r}_{1}\sigma_{1})\psi_{1}(\mathbf{r}_{2}\sigma_{2})\left[\vphantom{\sum}\chi_{0}(q_{2})\chi_{1}(q_{1})-\chi_{1}(q_{2})\chi_{0}(q_{1})\right]$
$\displaystyle{}={}$ $\displaystyle 0.$

Thus, we obtain the desired result of excluding the solution that violates the
Pauli principle. We conclude that we have reached our goal and found an exact
representation of the coupled electron-photon ground state in terms of
polariton orbitals.

#####  4.2.4 Polariton symmetry and wave functions

The next step is to construct a theory based on these polaritons.
Unfortunately, the constitutive relation (4.26) cannot be enforced on the wave
function level in a practical way. For that, we would need to generalize the
concept of a Slater determinant to polaritonic _and_ electronic coordinates,
as we have done with $\tilde{\Psi}^{\prime}$ to construct
$\tilde{\Psi}^{\prime}_{hybrid}$. In the general case, such a “generalized
Slater determinant” for two particles is given by

$\displaystyle\Psi_{ab}^{\prime}(\mathbf{z}_{1},\sigma_{1},\mathbf{z}_{2},\sigma_{2})=$
$\displaystyle\phi_{a}(\mathbf{r}_{1},q_{1},\sigma_{1})\phi_{b}(\mathbf{r}_{2},q_{2},\sigma_{2})-\phi_{b}(\mathbf{r}_{1},q_{1},\sigma_{1})\phi_{a}(\mathbf{r}_{2},q_{2},\sigma_{2})$
$\displaystyle+$
$\displaystyle\phi_{a}(\mathbf{r}_{1},q_{2},\sigma_{1})\phi_{b}(\mathbf{r}_{2},q_{1},\sigma_{2})-\phi_{b}(\mathbf{r}_{1},q_{2},\sigma_{1})\phi_{a}(\mathbf{r}_{2},q_{1},\sigma_{2}),$
(4.28)

where $\phi_{a/b}$ are some (orthonormal) polariton orbitals. To see why such
an ansatz is problematic, we calculate the norm-square of
$\Psi_{ab}^{\prime}$, i.e.,

$\displaystyle||\Psi_{ab}^{\prime}||^{2}=$
$\displaystyle\sum_{\sigma_{1},\sigma_{2}}\int\mathrm{d}\mathbf{z}_{1}\mathbf{z}_{2}\Psi_{ab}^{\prime}{}^{*}(\mathbf{z}_{1},\sigma_{1},\mathbf{z}_{2},\sigma_{2})\Psi_{ab}^{\prime}(\mathbf{z}_{1},\sigma_{1},\mathbf{z}_{2},\sigma_{2})$
$\displaystyle=$ $\displaystyle
4\braket{\phi_{a}}{\phi_{a}}\braket{\phi_{b}}{\phi_{b}}-4\braket{\phi_{a}}{\phi_{b}}\braket{\phi_{b}}{\phi_{a}}$
$\displaystyle+\sum_{\sigma_{1},\sigma_{2}}\int\mathrm{d}\mathbf{z}_{1}\mathbf{z}_{2}\phi_{a}^{*}(\mathbf{r}_{1},q_{1},\sigma_{1})\phi_{b}^{*}(\mathbf{r}_{2},q_{2},\sigma_{2})\phi_{a}(\mathbf{r}_{1},q_{2},\sigma_{1})\phi_{b}(\mathbf{r}_{2},q_{1},\sigma_{2})+c.c.$
$\displaystyle-\sum_{\sigma_{1},\sigma_{2}}\int\mathrm{d}\mathbf{z}_{1}\mathbf{z}_{2}\phi_{a}^{*}(\mathbf{r}_{1},q_{1},\sigma_{1})\phi_{b}^{*}(\mathbf{r}_{2},q_{2},\sigma_{2})\phi_{b}(\mathbf{r}_{1},q_{2},\sigma_{1})\phi_{a}(\mathbf{r}_{2},q_{1},\sigma_{2})-c.c.$
$\displaystyle+\sum_{\sigma_{1},\sigma_{2}}\int\mathrm{d}\mathbf{z}_{1}\mathbf{z}_{2}\phi_{b}^{*}(\mathbf{r}_{1},q_{1},\sigma_{1})\phi_{a}^{*}(\mathbf{r}_{2},q_{2},\sigma_{2})\phi_{b}(\mathbf{r}_{1},q_{2},\sigma_{1})\phi_{a}(\mathbf{r}_{2},q_{1},\sigma_{2})+c.c.$
$\displaystyle-\sum_{\sigma_{1},\sigma_{2}}\int\mathrm{d}\mathbf{z}_{1}\mathbf{z}_{2}\phi_{b}^{*}(\mathbf{r}_{1},q_{1},\sigma_{1})\phi_{a}^{*}(\mathbf{r}_{2},q_{2},\sigma_{2})\phi_{a}(\mathbf{r}_{1},q_{2},\sigma_{1})\phi_{b}(\mathbf{r}_{2},q_{1},\sigma_{2})-c.c.$

In this expression, only the first line is what would appear in a standard
Slater determinant. Since we consider orthonormal orbitals, i.e.,
$\braket{\phi_{a}}{\phi_{b}}=\delta_{ab}$, the integrals of the first term are
all one and the integrals of the second term are all zero. Thus for an
ordinary Slater determinant, the orthonormality condition is sufficient to fix
its norm, independently of the orbitals. However, with the terms that stem
from the additional symmetry requirements, new “mixed-index” terms arise,
whenever one or more of the according orbitals have coordinates with different
indices. All these (in this case) 8 terms are actually _two-body_ like
integrals, because we cannot separate the integrations of one set of
coordinates as we have done in the first line. The computation of these
integrals is nontrivial and cannot be defined by an orthonormality condition.
This means that the norm of $\Psi_{ab}^{\prime}$ depends on the specific form
of $\phi_{a/b}$ and thus has to be calculated explicitly to normalize
$\Psi_{ab}^{\prime}$.

The number of such terms for an $N$-body generalized Slater determinant is
given by all permutations of polaritonic and electronic coordinates, which is
$N!^{2}$, minus the $N!$ “ordinary” terms. The total number of “mixed-index”
terms, $N!^{2}-N!$, grows therefore factorial with the number of particles.
The normalization (and in the same way also the calculation of expectation
values) of a wave function that explicitly exhibits the symmetry (4.26)
requires the numerical calculation of (over)exponentially many nontrivial
terms. Such an explicit Slater determinant ansatz is thus infeasible in
practice. We have found yet another exponential wall.

#####  4.2.5 Polariton ensembles

Comparing the challenges of the dressed wave function with the difficulties to
construct a “polariton-like” wave function within generalized MF theory (Sec.
3.1.2), it seems that any many-body description based on polaritons leads not
only to an exponential but even a “factorial” wall. However, the dressed
construction is more flexible than the standard description and we show in the
following how we can exploit this flexibility to construct simple but accurate
approximation strategies based on polaritonic orbitals. The key-role hereby
will be played by the (ensemble) 1RDM that allows the enforce exchange-
symmetry without the explicit use of, e.g., Slater determinants.

For that, we generalize the description from pure states (described by single
wave functions) to ensembles, characterized by the $N$-body density matrix
$\Gamma_{E}^{N}=\sum_{i}w_{i}\Psi_{i}^{*}\Psi_{i}$, cf. Eq. (2.89), where the
$\Psi_{i}$ are eigenstates of the Hamiltonian and the $w_{i}$ are weight
coefficients with $0\leq w_{i}\leq 1$ and $\sum_{i}w_{i}=1$ (see Sec. 2.4.1).
One important feature of the ensemble description is that it is easier to
connect RDMs (Eq. (2.77)) to $\Gamma_{E}^{N}$ than to a pure-state density
matrix $\Gamma^{N}=\Psi^{*}\Psi$. This connection is given by the (ensemble)
$N$-representability conditions. In other words, if a given matrix is bosonic
(fermionic) $N$-representable than a corresponding $N$-body density matrix
exists that is composed out of only (anti)symmetric pure state wave functions.
The idea is now to rule out unphysical states such as $\tilde{\Psi}^{\prime}$
by utilizing the $N$-representability conditions of the 1RDM instead of an
explicit ansatz such as the generalized Slater determinant (Eq. (4.2.4)).

We start by noting that we can straightforwardly generalize the ensemble
description to the dressed setting, if we utilize the eigenstates
$\Psi^{\prime}_{i}$ of $\hat{H}^{\prime}$. However, as we are only interested
in the pure ground state $\Psi^{\prime}$, we will only implicitly make use of
the ensemble description. This is comparable to the use of the ensemble
$N$-representability conditions in RDMFT (see Sec. 2.4.2). We therefore define
the _dressed 1RDM_

$\displaystyle\gamma[\Psi_{i}^{\prime}](\mathbf{z}\sigma,\mathbf{z}^{\prime}\sigma^{\prime})=\sum_{\sigma_{2}}\int\mathrm{d}\mathbf{z}_{2}\Psi^{\prime}{}^{*}(\mathbf{z}_{1}^{\prime}\sigma_{1}^{\prime},\mathbf{z}_{2}\sigma_{2})\Psi^{\prime}(\mathbf{z}_{1}\sigma_{1},\mathbf{z}_{2}\sigma_{2})$
(4.29)

in terms of the dressed wave function $\Psi^{\prime}$.999The generalization to
ensembles would simply consist in exchanging
$\Gamma^{N}{}^{\prime}=\Psi^{\prime}{}^{*}\Psi^{\prime}$ by
$\Gamma_{E}^{N}{}^{\prime}=\sum_{i}w_{i}\Gamma_{i}^{N}{}^{\prime}$ for a set
of pure state density matrices $\Gamma_{i}^{N}{}^{\prime}$. Since the pure
state $\Psi^{\prime}$ is an extreme case of an ensemble state, its Fermi
statistics _with respect to the polaritonic coordinates_ $\mathbf{z}\sigma$
(cf. Eq. (4.44)) is also apparent in $\gamma[\Psi^{\prime}]$ in the form of
fermionic $N$-representability conditions. By using the natural orbitals
$\phi_{i}$ and the natural occupation numbers $n_{i}$, which are defined by
the eigenvalue equation $n_{i}\phi_{i}=\hat{\gamma}\psi_{i}$, we represent
$\gamma$ in its diagonal form

$\displaystyle\gamma[{\Psi^{\prime}}](\mathbf{z}\sigma,\mathbf{z}^{\prime}\sigma^{\prime})=\sum_{i=1}^{\infty}n_{i}\,\phi_{i}^{*}(\mathbf{z}^{\prime}\sigma^{\prime})\phi_{i}(\mathbf{z}\sigma).$
(4.30)

The fermionic $N$-representable conditions become especially simple in this
representation and are given by (cf. Eq. (2.91))

$\displaystyle\begin{split}&0\leq\,n_{i}\leq 1,\quad\forall i\\\
&\textstyle\sum_{i}n_{i}=N.\end{split}$ (4.31)

However, we still need to take care that the antisymmetry remains in the
electronic subsystem. For that, we now go one step further and define from the
dressed 1RDM the _electronic 1RDM_

$\displaystyle\begin{split}\gamma_{e}[{\Psi^{\prime}}](\mathbf{r}\sigma,\mathbf{r}^{\prime}\sigma^{\prime})=&\int\mathrm{d}^{M}\mathbf{q}\;\gamma[{\Psi^{\prime}}](\mathbf{r}\mathbf{q}\sigma,\mathbf{r}^{\prime}\mathbf{q}\sigma^{\prime})\\\
=&\sum_{i=1}^{\infty}n_{i}^{e}\,\psi_{i}^{e}{}^{*}(\mathbf{r}^{\prime}\sigma^{\prime})\psi_{i}^{e}(\mathbf{r}\sigma)\end{split}$
(4.32)

by an integration over the photonic coordinates. In the second line, we
introduced the diagonal representation of $\gamma_{e}$ with the according
natural orbitals $\psi_{i}^{e}$ and the natural occupation numbers
$n_{i}^{e}$. The Fermi statistics with regard to _only the electronic
coordinates_ $\mathbf{r}\sigma$ thus becomes apparent by considering the
electronic natural occupation numbers $n^{e}_{i}$101010Note that the
normalization of $\gamma_{e}$ to the electron number $N$ is a direct
consequence of the dressed construction that considers exactly $N$ polaritons
for system with $N$ electrons.

$\displaystyle 0\leq n_{i}^{e}\leq\,1,\quad\forall i$ (4.33a)
$\displaystyle\textstyle\sum_{i}n_{i}^{e}=N.$ (4.33b)

With Eqs. (4.31) and (4.33) together with the definition of $\gamma_{e}$ in
terms of $\gamma$ (Eq. (4.32)), we have derived a set of _necessary
conditions_ to guarantee the hybrid statistics. Further conditions would be
necessary for a sufficient characterization of the set of polaritonic wave
functions that exhibit the exchange symmetry (4.26). This is comparable to the
$N$-representable conditions of the 2RDM [214], since the dressed 1RDM
$\gamma(\mathbf{z};\mathbf{z}^{\prime})=\gamma(\mathbf{r},q,\mathbf{r}^{\prime}q^{\prime})$
written in electronic and photonic variables depends on _four_ particle
coordinates exactly as the 2RDM. Importantly, the conditions (4.33) are
sufficient to guarantee the antisymmetry in the fermionic sector, which is the
only symmetry of the physical systems and thus the most important one. We
therefore focus in the following on the conditions (4.33).

Let us illustrate this new point of view with the example from before. The
correct auxiliary ground state $\Psi^{\prime}$ satisfies the conditions of Eq.
(4.33), since

$\displaystyle\gamma_{e}[{\Psi^{\prime}}]=$
$\displaystyle\sum_{\sigma_{2}}\int\mathrm{d}\mathbf{z}_{2}\mathrm{d}q_{1}{\Psi^{\prime}}^{*}(\mathbf{r}^{\prime}q_{1}\sigma^{\prime},\mathbf{z}_{2}\sigma_{2})\Psi^{\prime}(\mathbf{r}q_{1}\sigma,\mathbf{z}_{2}\sigma_{2})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{2}\psi_{i}^{*}(\mathbf{r}^{\prime}\sigma^{\prime})\psi_{i}(\mathbf{r}\sigma).$

The two electronic orbitals are the natural orbitals of
$\gamma_{e}[\Psi^{\prime}]$ with natural occupation numbers $n_{1}=n_{2}=1$.
If we do the same calculation with $\tilde{\Psi}^{\prime}$, we get instead

$\displaystyle\gamma_{e}[{\tilde{\Psi}^{\prime}}]=2\psi_{1}^{*}(\mathbf{r}^{\prime}\sigma^{\prime})\psi_{1}(\mathbf{r}\sigma),$

which violates the $N$-representability conditions of Eq. (4.33) and thus the
Pauli principle. For the wave functions of an interacting system, the
diagonalization of $\gamma_{e}[\Psi^{\prime}]$ will not be as trivial as for
this simple example, but nevertheless the conditions (4.33) are _sufficient_
to ensure the Pauli exclusion principle in the sense that maximally one
fermion can occupy a single quantum state (upper bound in Eq. (4.33a)).

Most importantly, we can employ these conditions to obtain a computationally
tractable procedure (presented in Sec. 5.1) to approximately ensure the
polariton statistics implied by Eq. (4.43), instead of the _factorially_
growing number of “mixed-index” orbitals.

####  4.3 Dressed construction: the general case

We recapitulate now the dressed construction for the general case of $N$
electrons that are coupled to $M$ photon modes. We follow hereby Ref. [2, Sec.
2]. We consider the setting of cavity QED with according Hamiltonian
(reordering its terms for the purposes of this section)

$\displaystyle\begin{split}\hat{H}=&\underbrace{\sum_{k=1}^{N}\left[-\tfrac{1}{2}\Delta_{\mathbf{r}_{k}}+v(\mathbf{r}_{k})\right]+\tfrac{1}{2}\sum_{k\neq
l}^{N}w(\mathbf{r}_{k},\mathbf{r}_{l})}_{\hat{H}_{m}=\hat{T}[t]+\hat{V}[v]+\hat{W}[w]}+\underbrace{\sum_{\alpha=1}^{M}\left(-\tfrac{1}{2}\tfrac{\partial^{2}}{\partial
p_{\alpha}^{2}}+\tfrac{\omega_{\alpha}^{2}}{2}p_{\alpha}^{2}\right)}_{\hat{H}_{ph}}\\\
&{}+\underbrace{\sum_{\alpha=1}^{M}-\omega_{\alpha}p_{\alpha}\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}}}_{\hat{H}_{I}}+\underbrace{\sum_{\alpha=1}^{M}\tfrac{1}{2}\left(\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}}\right)^{2}}_{\hat{H}_{d}}.\end{split}$
(cf. 1.37)

The first three terms constitute the usual matter Hamiltonian of quantum
mechanics, with the kinetic and external one-body parts, $\hat{T}[t]$ and
$\hat{V}[v]$, respectively, and the two-body interaction term $\hat{W}[w]$.
Here the kinetic term is the usual Laplacian
$t(\mathbf{r})=-\tfrac{1}{2}\Delta_{\mathbf{r}}$, the external potential
$v(\mathbf{r})$ is due to the attractive nuclei/ions and
$w(\mathbf{r},\mathbf{r}^{\prime})$ is the electron-electron repulsion.
Usually this is just taken as the free-space Coulomb interaction
$w(\mathbf{r},\mathbf{r}^{\prime})=1/|\mathbf{r}-\mathbf{r}^{\prime}|$, but
inside a cavity the interaction can be modified [127]. The fourth term
$\hat{H}_{ph}$ is the free field-energy of $M$ effective modes of the cavity.
The effective modes are characterized by their displacement coordinate
$p_{\alpha}$, frequency $\omega_{\alpha}$ and polarization vectors
$\bm{\lambda}_{\alpha}$. The latter include already the effective coupling
strength
$g_{\alpha}=|\bm{\lambda}_{\alpha}|\sqrt{\frac{\omega_{\alpha}}{2}}\propto
1/\sqrt{V}$ [48, 14] that is proportional to the inverse square-root of the
cavity mode volume V. In the dipole approximation, the coupling between light
and matter is described by the bilinear term $\hat{H}_{int}$ together with the
dipole self-energy term $\hat{H}_{self}$. Here the dipole operator is defined
by $\hat{\mathbf{D}}=\sum_{k=1}^{N}\mathbf{r}_{k}$.

The ground state of $\hat{H}$ is given by a wave function

$\displaystyle\Psi(\mathbf{x}_{1},...,\mathbf{x}_{N},p_{1},...,p_{M})$ (4.34)

that depends on N spin-spatial electron coordinates and M photon mode-
displacement coordinates and that is antisymmetric with respect to the
exchange of any two electron coordinates.

#####  4.3.1 The dressed auxiliary system

To turn this coupled electron-photon problem into an equivalent and exact
dressed problem we will follow three steps:

  1. 1.

For each mode $\alpha=1,\dotsc,M$ and all but the first electron
$i=2,\dotsc,N$, we introduce extra auxiliary coordinates $p_{\alpha,i}$. This
adds $(N-1)M$ extra degrees of freedom to the problem. In this higher-
dimensional auxiliary configuration space, we now consider wave functions
depending on $4N+NM$ coordinates, i.e.,

$\displaystyle\Psi^{\prime}(\mathbf{r}_{1}\sigma_{1},\dotsc,\mathbf{r}_{N}\sigma_{N},p_{1},\dotsc,p_{M},p_{1,2},\dotsc,p_{1,N},\dotsc,p_{M,2},\dotsc,p_{M,N}).$

Here and in the following, we will denote all quantities in the auxiliary
configuration space with a prime.

  2. 2.

We next construct an auxiliary Hamiltonian in the extended configuration space
of the form

$\displaystyle\hat{H}^{\prime}=\hat{H}+\sum_{\alpha=1}^{M}\hat{\Pi}_{\alpha}$

where

$\displaystyle\hat{\Pi}_{\alpha}=\sum_{i=2}^{N}\biggl{(}-\frac{1}{2}\frac{\partial^{2}}{\partial
p_{\alpha,i}^{2}}+\frac{\omega_{\alpha}^{2}}{2}p_{\alpha,i}^{2}\biggr{)}$
(4.35)

depends only on these new auxiliary coordinates. This construction guarantees
that the auxiliary degrees of freedom do not mix with the physical ones, which
will ensure a simple and explicit connection between the physical and
auxiliary system.

  3. 3.

Finally we perform an orthogonal coordinate transformation of the physical and
auxiliary photon coordinates
$(p_{1},...,p_{M,N})\rightarrow(q_{1,1},...,q_{M,N})$ such that

$\displaystyle p_{\alpha}=\tfrac{1}{\sqrt{N}}$
$\displaystyle\left(q_{\alpha,1}+\dotsb+q_{\alpha,N}\right),$
$\displaystyle-\frac{1}{2}\frac{\partial^{2}}{\partial
p_{\alpha}^{2}}+\frac{\omega_{\alpha}^{2}}{2}p_{\alpha}^{2}+\hat{\Pi}_{\alpha}$
$\displaystyle=\sum_{i=1}^{N}\biggl{(}-\frac{1}{2}\frac{\partial^{2}}{\partial
q_{\alpha,i}^{2}}+\frac{\omega_{\alpha}^{2}}{2}q_{\alpha,i}^{2}\biggr{)}.$
(4.36)

Note that the second line is automatically satisfied for any orthogonal
transformation and the first line defines $p_{\alpha}$ as the “center-of-mass”
of all the $q_{\alpha,i}$ with uniform relative masses $1/\sqrt{N}$.

In total, we then find the auxiliary Hamiltonian in the higher-dimensional
configuration space given as

$\displaystyle\hat{H}^{\prime}{}={}$
$\displaystyle\sum_{k=1}^{N}\left[-\tfrac{1}{2}\Delta_{\mathbf{r}_{k}}^{2}+v(\mathbf{r}_{k})\right]+\tfrac{1}{2}\sum_{k\neq
l}w(\mathbf{r}_{k},\mathbf{r}_{l})-\sum_{\alpha=1}^{M}\omega_{\alpha}p_{\alpha}\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}}+\sum_{\alpha=1}^{M}\tfrac{1}{2}\left(\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}}\right)^{2}$
$\displaystyle+\sum_{\alpha=1}^{M}\left(-\tfrac{1}{2}\tfrac{\partial^{2}}{\partial
p_{\alpha}^{2}}+\tfrac{\omega_{\alpha}^{2}}{2}p_{\alpha}^{2}\right)+\sum_{\alpha=1}^{M}\sum_{i=2}^{N}\left(-\tfrac{1}{2}\tfrac{\partial^{2}}{\partial
p_{\alpha,i}^{2}}+\tfrac{\omega_{\alpha}^{2}}{2}p_{\alpha,i}^{2}\right)$
$\displaystyle{}\overset{\eqref{eq:CenterOfMass}}{=}{}$
$\displaystyle\sum_{k=1}^{N}\left\\{-\tfrac{1}{2}\Delta_{\mathbf{r}_{k}}+v(\mathbf{r}_{k})+\sum_{\alpha=1}^{M}\\!\left[-\tfrac{1}{2}\tfrac{\partial^{2}}{\partial
q_{\alpha,k}^{2}}+\tfrac{1}{2}\omega_{\alpha}^{2}q_{\alpha,k}^{2}-\tfrac{\omega_{\alpha}}{\sqrt{N}}q_{\alpha,k}(\bm{\lambda}_{\alpha}\\!\cdot\mathbf{r}_{k})+\tfrac{1}{2}(\bm{\lambda}_{\alpha}\\!\cdot\mathbf{r}_{k})^{2}\right]\right\\}$
$\displaystyle+\tfrac{1}{2}\sum_{k\neq
l}\left[w(\mathbf{r}_{k},\mathbf{r}_{l})+\sum_{\alpha=1}^{M}\left(-\tfrac{\omega_{\alpha}}{\sqrt{N}}q_{\alpha,k}\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{l}-\tfrac{\omega_{\alpha}}{\sqrt{N}}q_{\alpha,l}\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{k}+\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{k}\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{l}\right)\right],$

where we inserted the definition of the total dipole operator and reordered
the expressions, such that the terms with only one index and the terms with
two different indices are grouped together. Introducing then a
$(3+M)$-dimensional polaritonic vector of space and transformed photon
coordinates $\mathbf{z}=\mathbf{r}\mathbf{q}$ with
$\mathbf{q}\equiv(q_{1},\dotsc,q_{M})$, we can rewrite the above Hamiltonian
as

$\displaystyle\hat{H}^{\prime}$
$\displaystyle=\sum_{k=1}^{N}\left[-\tfrac{1}{2}\Delta_{k}^{\prime}+v^{\prime}(\mathbf{z}_{k})\right]+\tfrac{1}{2}\sum_{k\neq
l}w^{\prime}(\mathbf{z}_{k},\mathbf{z}_{l})$ (4.37)
$\displaystyle=\hat{T}[t^{\prime}]+\hat{V}[v^{\prime}]+\hat{W}[w^{\prime}],$

where we introduced the dressed one-body terms

$\displaystyle t^{\prime}(\mathbf{z})$
$\displaystyle=-\tfrac{1}{2}\Delta_{k}^{\prime}\equiv-\tfrac{1}{2}\sum_{i=1}^{3}\tfrac{\partial^{2}}{\partial
r_{i}^{2}}-\tfrac{1}{2}\sum_{\alpha=1}^{M}\tfrac{\partial^{2}}{\partial
q_{\alpha}^{2}},$ (4.38) $\displaystyle v^{\prime}(\mathbf{z})$
$\displaystyle=v(\mathbf{r})+\sum_{\alpha=1}^{M}\left[\tfrac{1}{2}\omega_{\alpha}^{2}q_{\alpha}^{2}-\tfrac{\omega_{\alpha}}{\sqrt{N}}q_{\alpha}\bm{\lambda}_{\alpha}\\!\\!\cdot\\!\mathbf{r}+\tfrac{1}{2}(\bm{\lambda}_{\alpha}\\!\\!\cdot\\!\mathbf{r})^{2}\right],$
(4.39) and the dressed two-body interaction term $\displaystyle
w^{\prime}(\mathbf{z},\mathbf{z}^{\prime})=$ $\displaystyle
w(\mathbf{r},\mathbf{r}^{\prime})+\sum_{\alpha=1}^{M}\left[-\tfrac{\omega_{\alpha}}{\sqrt{N}}q_{\alpha}\bm{\lambda}_{\alpha}\\!\\!\cdot\\!\mathbf{r}^{\prime}-\tfrac{\omega_{\alpha}}{\sqrt{N}}q_{\alpha}^{\prime}\bm{\lambda}_{\alpha}\\!\\!\cdot\\!\mathbf{r}+\bm{\lambda}_{\alpha}\\!\\!\cdot\\!\mathbf{r}\bm{\lambda}_{\alpha}\\!\\!\cdot\\!\mathbf{r}^{\prime}\right].$
(4.40)

We see here that only the conditions (3), but not the details of the
coordinate transformation of step 3 are important for our construction [1].
The crucial part of this coordinate transformation is the replacement of
$p_{\alpha}$ in the interaction terms
$p_{\alpha}\bm{\lambda}_{\alpha}\\!\cdot\\!\hat{\mathbf{D}}$. Instead of
$p_{\alpha}$ only, now all $q_{\alpha,i}$ couple to the dipole of the matter
system just with a rescaled coupling-strength by the factor $1/\sqrt{N}$.

###### The hybrid statistics of the dressed wave function

Let us next discuss the wave function $\Psi^{\prime}$ in the auxiliary
configuration space. The wave function $\Psi$ in the usual configuration space
is a (normalized) solution of the (time-independent) Schrödinger equation
$E_{0}\Psi=\hat{H}\Psi$. Since
$\hat{H}^{\prime}=\hat{H}+\sum_{\alpha=1}^{M}\hat{\Pi}_{\alpha}$ and
$\hat{\Pi}_{\alpha}$ acts only on the auxiliary coordinates, we can simply
construct

$\displaystyle\Psi^{\prime}(\mathbf{r}_{1}\sigma_{1},\dotsc,p_{M,N})=\Psi(\mathbf{r}_{1}\sigma_{1},\dotsc,\mathbf{r}_{N}\sigma_{N};p_{1},\dotsc,p_{M})\chi(p_{1,2},\dotsc,p_{M,N}),$
(4.41)

with $\chi$ being the (normalized) ground state of
$\sum_{\alpha=1}^{M}\hat{\Pi}_{\alpha}$, which is a product of individual
harmonic-oscillator ground states. Clearly, $\Psi^{\prime}$ is a normalized
solution of the auxiliary Schrödinger equation
$E_{0}^{\prime}\Psi^{\prime}=\hat{H}^{\prime}\Psi^{\prime}$. In principle any
combination of eigenstates of the auxiliary harmonic oscillators would lead to
a new eigenfunction for $\hat{H}^{\prime}$ but since we here focus on the
ground state the natural choice is the lowest-energy solution. Rewriting this
wave function in the new coordinates and employing the polaritonic coordinates
$\mathbf{z}\sigma\equiv\mathbf{r}\bm{q}\sigma$, we arrive at

$\displaystyle\Psi^{\prime}(\mathbf{r}_{1}\sigma_{1},\dotsc,\mathbf{r}_{N}\sigma_{N},q_{1,1},\dotsc,q_{M,N})=\Psi^{\prime}(\mathbf{z}_{1}\sigma_{1},\dotsc,\mathbf{z}_{N}\sigma_{N}).$
(4.42)

This polaritonic wave function as the ground state of (4.37) is the
reformulation of the original electron-photon problem of (1.37) we were
looking for. Since all the new photonic coordinates belong to harmonic
oscillator ground states, exchanging $p_{\alpha,i}$ with $p_{\alpha,j}$ does
not change the total wave function $\Psi^{\prime}$ and this property transfers
to the exchange of any coordinate $q_{\alpha,i}$ and $q_{\alpha,j}$. Hence we
have now a _bosonic_ symmetry with respect to the $\bm{q}$ coordinates. Since
the electronic part of the auxiliary system is not affected by the coordinate
transformation, the electronic symmetries are the same in the physical and
auxiliary system, i.e., we have a _fermionic_ symmetry with respect to
$\mathbf{r}\sigma$. Together these two fundamental symmetries imply that the
polaritonic coordinates $\mathbf{z}\sigma$ have fermionic character. The
symmetries of the polaritonic wave function $\Psi^{\prime}$ can be summarized
as

$\displaystyle\mathbf{r}_{k}\sigma_{k}$
$\displaystyle\leftrightarrow\mathbf{r}_{l}\sigma_{l}\quad$
$\displaystyle\rightarrow\quad\Psi^{\prime}$
$\displaystyle\leftrightarrow-\Psi^{\prime}$ (4.43a) $\displaystyle\bm{q}_{k}$
$\displaystyle\leftrightarrow\bm{q}_{l}$
$\displaystyle\rightarrow\quad\Psi^{\prime}$
$\displaystyle\leftrightarrow\Psi^{\prime}$ (4.43b)

from which follows

$\displaystyle\mathbf{z}_{k}\sigma_{k}\leftrightarrow\mathbf{z}_{l}\sigma_{l}\quad$
$\displaystyle\rightarrow\quad\Psi^{\prime}\leftrightarrow-\Psi^{\prime}$
(4.44)

This means that though the dressed wave function has fermionic statistics
(4.44) in terms of the polaritonic coordinates $\mathbf{z}\sigma$, due to the
_constitutive relations_ (4.43) it actually consists of two types of
particles: one with fermionic character and another with bosonic character.
Consequently, the polariton wave function $\Psi^{\prime}$ has a hybrid Fermi-
Bose statistics. As consequences of these symmetries we find the Pauli
exclusion principle for the electrons, yet for the auxiliary photon
coordinates we find that many photonic auxiliary entities can occupy the same
quantum state.

Indeed, we can prove that the conditions of Eq. (4.43) are necessary and
sufficient to establish a one-to-one mapping between the dressed
$\Psi^{\prime}$ and the physical ground states $\Psi$. The first part
“$\Psi\rightarrow\Psi^{\prime}$” is given by the dressed construction. For the
second part “$\Psi^{\prime}\rightarrow\Psi$” we show that the minimal energy
state $\Psi^{\prime}$ has the form $\Psi^{\prime}=\Psi\chi$, cf. Eq. (4.41).
For that, we consider a trial wave function in the dressed space

$\displaystyle\Upsilon^{\prime}(\mathbf{z}_{1}\sigma_{1},...,\mathbf{z}_{N}\sigma_{N})\equiv\Upsilon^{\prime}(\mathbf{r}_{1}\sigma_{1},...,p_{M};p_{1,2},...,p_{M,N})$

that fulfils (4.43). Then it holds since $\hat{H}$ only acts on
$(\mathbf{r}_{1},\dotsc,\mathbf{r}_{N},p_{1},\dotsc,p_{M})$ and
$\sum_{\alpha=1}^{M}\hat{\Pi}_{\alpha}$ only acts on $p_{1,2},...,p_{M,N}$
that

$\displaystyle\inf_{\Upsilon^{\prime}}\braket{\Upsilon^{\prime}}{\hat{H}^{\prime}\Upsilon^{\prime}}$
$\displaystyle\,\,\,\,\,\geq\inf_{\Upsilon^{\prime}}\braket{\Upsilon^{\prime}}{\hat{H}\Upsilon^{\prime}}+\inf_{\Upsilon^{\prime}}\braket{\Upsilon^{\prime}}{\sum_{\alpha=1}^{M}\hat{\Pi}_{\alpha}\Upsilon^{\prime}}$
$\displaystyle\overset{\eqref{eq:SymmetryPauliDressedFermion}}{=}\braket{\Psi}{\hat{H}\Psi}+\braket{\chi}{\sum_{\alpha=1}^{M}\hat{\Pi}_{\alpha}\chi}$
$\displaystyle\,\,\,\,\,=\braket{\Psi^{\prime}}{\hat{H}^{\prime}\Psi^{\prime}}.$

We have thus proven the assumption. For excited states the constitutive
relations are necessary but not sufficient to single out the eigenfunctions of
$\hat{H}^{\prime}$ that correspond to the original Hamiltonian in terms of
simple products. The reason is that we cannot make use of the variational
principle as above. Let us illustrate this for the first excited state
$\Psi_{1}$ with corresponding dressed version
$\Psi_{1}^{\prime}=\Psi_{1}\chi$. We could try to follow a similar line of
proof as above by considering only trail wave functions
$\tilde{\Upsilon}^{\prime}\perp\Psi^{\prime}$ that are orthogonal to the
ground state. If we now search for the lowest energy solution among the
$\tilde{\Upsilon}^{\prime}$, the constitutive relations cannot differentiate
the correct state $\Psi_{1}^{\prime}=\Psi_{1}\chi$ from, e.g., the state
$\tilde{\Psi}_{1}^{\prime}=\Psi\chi_{1}$, where $\chi_{1}$ is the first
excited state of the auxiliary modes.111111The reader is referred to the work
of Nielsen et al. [258] for further details on excited states and the time-
dependent case in general. However, for ground states the above proof and thus
the one-to-one correspondence between physical and dressed system holds.

#####  4.3.2 Enforcing the hybrid statistics in practice: the Pauli principle
and $N$-representability

We have discussed in the previous section that enforcing the physical
conditions (4.43) on the polaritonic wave function directly, is not practical.
As an alternative, the physical conditions (4.43) are visible in the dressed
1RDM, which is given explicitly by

$\displaystyle\gamma[{\Psi^{\prime}}](\mathbf{z}\sigma,\mathbf{z}^{\prime}\sigma^{\prime})=\sum_{\sigma_{2},\dotsc,\sigma_{N}}\int\mathrm{d}^{3(N-1)}\mathbf{z}\Psi^{\prime}{}^{*}(\mathbf{z}^{\prime}\sigma^{\prime},\dotsc,\mathbf{z}_{N}\sigma_{N})\Psi^{\prime}(\mathbf{z}\sigma,\dotsc,\mathbf{z}_{N}\sigma_{N}).$
(4.45)

The Fermi statistics of the wave function $\Psi^{\prime}$ with respect to the
polaritonic coordinates $\mathbf{z}\sigma$ (4.44) is also apparent in
$\gamma[\Psi^{\prime}]$ in the form of the $N$-representability conditions. By
using the _natural orbitals_ $\phi_{i}$ and the _natural occupation numbers_
$n_{i}$, which are defined by the eigenvalue equation
$n_{i}\,\phi_{i}=\hat{\gamma}\,\psi_{i}$, we represent $\gamma$ in its
diagonal form
$\gamma[{\Psi^{\prime}}](\mathbf{z}\sigma,\mathbf{z}^{\prime}\sigma^{\prime})=\sum_{i=1}^{\infty}n_{i}\,\phi_{i}^{*}(\mathbf{z}^{\prime}\sigma^{\prime})\phi_{i}(\mathbf{z}\sigma)$.
The fermionic $N$-representability conditions become especially simple in this
representation and are given by

$\displaystyle\begin{split}&0\leq\,n_{i}\leq 1,\quad\forall i\\\
&\textstyle\sum_{i}n_{i}=N.\end{split}$ (4.46)

From the dressed 1RDM, we can define the electronic 1RDM

$\displaystyle\begin{split}\gamma_{e}[{\Psi^{\prime}}](\mathbf{r}\sigma,\mathbf{r}^{\prime}\sigma^{\prime})=\int\mathrm{d}^{M}\mathbf{q}\;\gamma[{\Psi^{\prime}}](\mathbf{r}\mathbf{q}\sigma,\mathbf{r}^{\prime}\mathbf{q}\sigma^{\prime})\end{split}$
(4.47)

and the auxiliary photonic 1RDM

$\displaystyle\gamma_{p}[{\Psi^{\prime}}](\mathbf{q},\mathbf{q}^{\prime})=\sum_{\sigma}\int\mathrm{d}^{3}\mathbf{r}\,\gamma[{\Psi^{\prime}}](\mathbf{r}\mathbf{q}\sigma,\mathbf{r}\mathbf{q}^{\prime}\sigma).$
(4.48)

Again, we can define the according natural orbitals $\psi_{i}^{e/p}$ and the
natural occupation numbers $n_{i}^{e/p}$ by the eigenvalue equations
$n_{i}\psi_{i}^{e/p}=\hat{\gamma}_{e/p}\psi_{i}^{e/p}$ and go into their
diagonal representations

$\displaystyle\gamma_{e}[{\Psi^{\prime}}](\mathbf{r}\sigma,\mathbf{r}^{\prime}\sigma^{\prime})=\sum_{i=1}^{\infty}n_{i}^{e}\,\psi_{i}^{e}{}^{*}(\mathbf{r}^{\prime}\sigma^{\prime})\psi_{i}^{e}(\mathbf{r}\sigma)$
(4.49)

and

$\displaystyle\gamma_{p}[{\Psi^{\prime}}](\mathbf{q},\mathbf{q}^{\prime})=\sum_{i=1}^{\infty}n_{i}^{p}\,\psi_{i}^{p}{}^{*}(\mathbf{q}^{\prime})\psi_{i}^{p}(\mathbf{q}).$
(4.50)

The Fermi statistics with regard to only the electronic coordinates
$\mathbf{r}\sigma$ thus becomes apparent by considering the electronic natural
occupation numbers $n^{e}_{i}$

$\displaystyle n_{i}^{e}\geq\,0,\quad\forall i$ (4.51a) $\displaystyle
n_{i}^{e}\leq 1,\quad\forall i$ (4.51b)
$\displaystyle\textstyle\sum_{i}n_{i}^{e}=N,$ (4.51c)

where we split the conditions in three parts for later convenience. The
equivalent bosonic symmetry of the auxiliary photonic coordinates leads
instead to the conditions

$\displaystyle 0\leq\,n_{i}^{p},\quad\forall i$ (4.52a)
$\displaystyle\textstyle\sum_{i}n_{i}^{p}=N.$ (4.52b)

Note that the normalization of $\gamma_{e/b}$ to the electron number $N$ is a
direct consequence of the auxiliary construction that considers exactly $N$
polaritons for a system with $N$ electrons. This becomes explicitly visible in
the fact that the normalization of $\gamma$ by definition transfers to
$\gamma_{e/b}$, since
$N=\sum_{\sigma}\int\mathrm{d}\mathbf{z}\gamma(\mathbf{z}\sigma,\mathbf{z}\sigma)=\sum_{\sigma}\int\mathrm{d}\mathbf{r}\gamma_{e}(\mathbf{r}\sigma,\mathbf{r}\sigma)=\int\mathrm{d}\mathbf{q}\gamma_{p}(\mathbf{q},\mathbf{q})$.
Additionally, the lower bounds of $\gamma_{e/b}$, cf. Eqs. (4.51a) and
(4.52a), transfer from $\gamma$, because the partial trace operation is a
completely positive map [263]. We can conclude that if (4.46) is enforced,
only the upper bound of the electronic 1RDM, cf. Eq. (4.51b) provides a
nontrivial additional constrained.

This now shows explicitly also for an interacting wave function that at most
one electron can occupy a specific quantum state, while many auxiliary photon
quantities can occupy a single quantum state. Further, the dressed 1RDM
$\gamma[\Psi^{\prime}]$ itself has only natural occupation numbers between
zero and one and is therefore fermionic, yet it contains a fermionic and a
bosonic subsystem. It is important to note that the original wave function
$\Psi$ did not have this simple hybrid statistics but only fermionic symmetry,
since the physical $p_{\alpha}$ did not follow any specific statistics.
Further, that we genuinely have formulated the coupled electron-photon problem
in terms of hybrid quasi-particles becomes most evident by actually using
single-particle (polariton) orbitals $\phi_{i}(\mathbf{z}\sigma)$ to expand
the dressed 1RDM of $\Psi^{\prime}$.

To obtain a computationally tractable procedure, we therefore use the
construction presented in Sec. 5.1 to (approximately) ensure the polariton
statistics implied by Eq. (4.43), instead of the _factorially_ growing number
of “mixed-index” orbitals. We will consider all fermionic density matrices in
the auxiliary configuration space, which we characterize by the conditions of
Eq. (4.46) in terms of polaritonic orbitals $\phi_{i}(\mathbf{z}\sigma)$. We
then constrain this space by enforcing the $N$-representability conditions of
Eq. (4.51) for the 1RDM of the electronic subsystem. Since this guarantees
that only (ensembles of) fermionic wave functions are allowed, also the
minimal energy solution (corresponds to an ensemble that) has fermionic
symmetry with respect to $\mathbf{r}\sigma$. We remind the reader that this is
the only symmetry of the physical system. Thus, it is especially important to
enforce this symmetry also within an accurate approximation scheme in the
dressed system. This together with the $\mathbf{z}\sigma$ antisymmetry
guarantees additionally the correct zero-coupling limit. We call this
construction the polariton ansatz for strong light-matter interaction. In the
next section, we will based on the polariton ansatz provide a detailed
prescription to generalize a given electronic-structure theory to treat ground
states of coupled electron-photon systems from first principles.

#####  4.3.3 Observables in the dressed system

To conclude we want to briefly discuss the role of observables in the
auxiliary system. Although we constructed the auxiliary space explicitly in a
way that the physical wave function $\Psi$ can be reconstructed exactly from
its dressed counterpart $\Psi^{\prime}$ by integration of all auxiliary
coordinates, this does not hold for all types of operators. For operators that
depend only on electronic coordinates, there is no difference and we have
$\braket{\Psi}{\hat{O}\Psi}=\braket{\Psi^{\prime}}{\hat{O}\Psi^{\prime}}$.
This is not surprising because the coordinate transformation (3) acts only on
the photonic part of the system. For photonic observables instead, the
transformation changes the respective operators and thus, the connection
between physical and auxiliary space becomes nontrivial in general. However,
at least for all observables that depend on photonic $1/2$\- or 1-body
expressions, there is an analytical connection. For half-body operators, i.e.,
any operator that depends only on the displacement of

$\displaystyle p_{\alpha}=\frac{1}{\sqrt{N}}(q_{\alpha,1}+...+q_{\alpha,N})$
(4.53)

and its conjugate

$\displaystyle\frac{\partial}{\partial
p_{\alpha}}=\frac{1}{\sqrt{N}}\left(\frac{\partial}{\partial
q_{\alpha,1}}+...+\frac{\partial}{\partial q_{\alpha,N}}\right),$ (4.54)

the coordinate transformation itself provides us with the connection. For
1-body operators, this becomes slightly more involved. For example, consider
the mode energy operator
$\hat{H}_{ph}=\sum_{\alpha=1}^{M}\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
p_{\alpha}^{2}}+\frac{\omega_{\alpha}^{2}}{2}p_{\alpha}^{2}\right]\equiv\sum_{\alpha=1}^{M}\hat{h}_{\alpha}$,
that we can straightforwardly generalize in the auxiliary space to
$\hat{H}_{ph}^{\prime}=\sum_{\alpha=1}^{M}\sum_{i=1}^{N}\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
q_{\alpha,i}^{2}}+\frac{\omega_{\alpha}^{2}}{2}q_{\alpha,i}^{2}\right]$
$\equiv\sum_{\alpha=1}^{M}\hat{h}_{\alpha}^{\prime}$. The connection between
$\hat{H}_{ph}$ and $\hat{H}_{ph}^{\prime}$ is given by the definition of the
coordinate transformation (3),

$\displaystyle\hat{H}_{ph}=\hat{H}_{ph}^{\prime}-\sum_{\alpha=1}^{M}\hat{\Pi}_{\alpha}.$
(4.55)

Since the expectation value of $\hat{\Pi}_{\alpha}$ is known analytically,
$\braket{\Psi^{\prime}}{\hat{\Pi}_{\alpha}(p_{\alpha,2},...,p_{\alpha,N})\Psi^{\prime}}=\braket{\chi}{\hat{\Pi}_{\alpha}\chi}=(N-1)\sum_{\alpha=1}^{M}\frac{\omega_{\alpha}}{2}$,
we have

$\displaystyle E_{ph}=$
$\displaystyle\braket{\Psi}{\hat{H}_{ph}\Psi}=\braket{\Psi^{\prime}}{\hat{H}_{ph}^{\prime}\Psi^{\prime}}-(N-1)\sum_{\alpha=1}^{M}\frac{\omega_{\alpha}}{2},$
(4.56) that enters the total energy $\displaystyle E=$
$\displaystyle\braket{\Psi}{\hat{H}\Psi}=\braket{\Psi^{\prime}}{\hat{H}^{\prime}\Psi^{\prime}}-(N-1)\sum_{\alpha=1}^{M}\frac{\omega_{\alpha}}{2}.$
(4.57)

This can be generalized to any operator that contains terms of the form
$p_{\alpha}p_{\beta},p_{\alpha}\frac{\partial}{\partial p_{\beta}},$ and
$\frac{\partial}{\partial p_{\alpha}}\frac{\partial}{\partial p_{\beta}}$,
where $\alpha$ and $\beta$ denote any two modes. The reason is that the
transformation $\eqref{eq:CenterOfMass}$ preserves the standard inner product
of the Euclidean space of the mode plus extra coordinates. This transfers also
to their conjugates and combinations of both. From the above, it is
straightforward to derive also the expression for the occupation of mode
$\alpha$,
$\hat{N}^{\alpha}_{ph}=\frac{1}{\omega_{\alpha}}\hat{h}_{\alpha}-\frac{1}{2}$.
In the auxiliary system, we have

$\displaystyle\hat{N}^{\alpha}_{ph}=\frac{1}{\omega_{\alpha}}\hat{h^{\prime}}_{\alpha}-\frac{N}{2}.$
(4.58)

### Chapter 5 Polaritons from first principles

In the last chapter, we have defined the theoretical framework to describe
coupled matter-photon systems with polaritons as the fundamental entity. Now
we need to discuss how to use the framework in practice. The key behind the
application of the dressed construction lies in the fact that the polaritons
are “almost electrons” - with an additional symmetry - and the Hamiltonian in
terms of these polaritons is also almost the Hamiltonian of standard
electronic-structure theory. In this chapter, we lay out how we can exploit
this by providing a general prescription of how to turn _any_ given
electronic-structure theory into a “polaritonic-structure theory” (Sec. 5.1).
This means that we keep the approximations to the Coulomb-interaction within
the electronic theory exactly as they are but use them to approximate the in
the previous section derived interaction between the polaritons. Sec. 5.1 is
based on chapter 3 of [2]. Then, we exemplify this prescription for two
examples. In Sec. 5.2, we introduce polaritonic HF, followed by polaritonic
RDMFT in Sec. 5.3.

####  5.1 Polaritonic structure theory

In this section, we lay out in detail how one can transform a given
electronic-structure theory that meets some minimal requirements into its
polaritonic version. The goal of such a “polaritonic-structure theory” is to
find the ground state of the cavity-QED Hamiltonian of Eq. (1.37) by
considering the ground state of the auxiliary Hamiltonian of Eq. (4.37). Let
us start by defining the according variational principle for the ground-state
energy $E_{0}^{\prime}$ as

$\displaystyle
E_{0}^{\prime}=\inf_{\Psi^{\prime}\,\in\,\mathfrak{P}}\braket{\Psi^{\prime}}{\hat{H}^{\prime}\Psi^{\prime}},$
(5.1)

where
$\mathfrak{P}=\\{\Psi^{\prime}:\Psi^{\prime}\leftrightarrow\text{\eqref{eq:SymmetryPauliDressedFermion}}\\}$
is the set of all normalized many-polariton wave functions that obey the
constitutive relations of Eq. (4.43). For our purposes, as explained in Sec.
4.3, we will instead consider the larger set $\mathfrak{M}$ of all (mixed-
state) density matrices
$\Gamma=\sum_{j}w_{j}\ket{\Psi^{\prime}_{j}}\bra{\Psi_{j}^{\prime}}$ with
$\sum_{j}w_{j}=1$, that obey the hybrid Fermi-Bose statistics. The minimal
energy also in this more general set corresponds to the pure state of Eq.
(5.1), i.e.,

$\displaystyle
E_{0}^{\prime}=\inf_{\Gamma\,\in\,\mathfrak{M}}\operatorname{Tr}\\{\Gamma\hat{H}\\}.$
(5.2)

The main trick now is in how we approximate this set. We do so by first
considering the yet larger set
$\tilde{\mathfrak{M}}=\\{\tilde{\Gamma}:\tilde{\Psi^{\prime}_{j}}=\sum
C_{j,K}\Phi_{K}\\}$, i.e., density matrices made of superpositions of Slater
determinants $\Phi_{K}=\text{det}(\phi_{K,1}\cdots\phi_{K,N})/\sqrt{N!}$ of
polariton orbitals $\phi_{K,i}$. This guarantees the overall Fermi statistics
in terms of the polaritonic coordinates $\mathbf{z}\sigma$. We then constrain
this larger set to

$\displaystyle\mathcal{M}^{\prime}=\\{\tilde{\Gamma}\in\tilde{\mathfrak{M}}:n^{e}_{i}[{\tilde{\Gamma}}]\leq
1\\},$ (5.3)

where $n^{e}_{i}[{\tilde{\Gamma}}]$ are the natural occupation numbers, cf.
Eq. (4.49), of the electronic 1RDM
$\gamma_{e}[\tilde{\Gamma}]=\sum_{j}w_{j}\gamma_{e}[\tilde{\Psi^{\prime}_{j}}]$
that depend on $\tilde{\Gamma}$. This enforces the fermionic statistics with
respect to the electronic coordinates $\mathbf{r}\sigma$. The rest of the
$N$-representability conditions (Eqs. (4.51a) and (4.51c)) are satisfied
automatically by choosing $\Psi^{\prime}\in\tilde{\mathcal{P}}$ as fermionic
with respect to polariton coordinates and thus the corresponding dressed 1RDM
satisfies the $N$-representability conditions (Eq. (4.46)). This guarantees
that the electronic and photonic 1RDMs of the system are $N$-representable,
and thus, e.g., the electronic Pauli-principle is enforced. However, higher
order RDMs are not treated exactly, which is an interesting topic for future
research. We thus avoid the direct construction of the exponentially growing
correlated electron-photon states.

As we have pointed out already, the polariton picture gives any coupled
problem that is described by Hamiltonian (1.37) the _same structure_ as a
purely electronic problem with two-body interactions, i.e., the electronic-
structure Hamiltonian (Eq. 2.1). Consequently, we can transfer every type of
electronic-structure theory to the coupled electron-photon problem, if the
theory provides an expression for the 1RDM (since we need the 1RDM to test the
$N$-representability constraints).

Hamiltonian Basis States Lagrangian of the EST
$\hat{H}_{m}=\hat{T}[t]+\hat{V}[v]+\hat{W}[w]$
$\\{\phi_{k}(\mathbf{r}\sigma)\\}$ $\rightarrow$
$L_{m}[\Psi]=E_{m}^{t,v,w}[\Psi]+\mathcal{C}[\Psi]$ $\downarrow$
$\hat{H}=\hat{H}_{m}+\hat{H}_{ph}+\hat{H}_{I}+\hat{H}_{d}$
$\\{\phi_{k}(\mathbf{r}\sigma),\chi_{\alpha}(p)\\}$ $\rightarrow$ “new theory
necessary” $\downarrow$
$\hat{H}^{\prime}=\hat{T}[t^{\prime}]+\hat{V}[v^{\prime}]+\hat{W}[w^{\prime}]$
$\\{\phi_{k}^{\prime}(\mathbf{r}\sigma,q)\\}$ $\rightarrow$
$L^{\prime}[\Psi^{\prime}]=E_{m}^{t^{\prime},v^{\prime},w^{\prime}}[\Psi^{\prime}]+\mathcal{C}[\Psi^{\prime}]+\mathcal{G}[\Psi^{\prime}]$

Figure 5.1: Graphical illustration of the polariton construction and its
connection to an electronic-structure theory (EST). Here $E_{m}$ indicates the
energy expression of the EST, such as the HF, configuration interaction, or
coupled cluster energy functional, and $\mathcal{C}$ indicates the constraints
of the EST, such as orthonormality of the orbitals. They are enforced on a
(possibly multideterminantal) wave function $\Psi$ constructed from an
electronic single-particle basis $\phi_{k}$. Further, $\mathcal{G}$ indicates
the new constraints that arise due to the hybrid statistics of the polaritons,
which are now enforced on a (possibly multideterminantal) wave function
$\Psi^{\prime}$ of polaritonic single-particle orbitals $\phi^{\prime}_{k}$.
The usual coupled electron-photon problem (second line) has a Hamiltonian with
a different structure and is built on separate orbitals $\phi_{k}$ and
$\chi_{\alpha}$. Thus, a new (efficient and accurate) approximate energy
expression would be needed.

The main steps how to do so are depicted in Fig. 5.1. We assume that the
theory provides us an energy expression $E_{m}$ with respect to a set of
electronic basis states $\\{\phi_{k}\\}$. This requirement is met by basically
every electronic-structure theory, as for instance the single-reference
methods KS-DFT or HF but also more involved (and numerically expensive)
approaches like coupled cluster, valence bond theory or configuration
interaction.111A standard reference for these and other quantum chemical
methods is Helgaker et al. [132]. Depending on the specific theory, $E_{m}$
might have quite different forms, but it is always derived from some many-body
Hamiltonian $\hat{H}_{m}=\hat{T}[t]+\hat{V}[v]+\hat{W}[w]$. More specifically,
the connection between $\hat{H}_{m}$ and $E_{m}$ is given by the particle
number $N$ and the integral kernels $(t,v,w)$ of the three energy operators.
For the matter Hamiltonian $\hat{H}_{m}$ of Eq. (1.37) for example, these
kernels are given by $t=-\frac{1}{2}\nabla^{2}_{\mathbf{r}}$,
$v=v(\mathbf{r})$ and $w=w(\mathbf{r},\mathbf{r}^{\prime})$. The goal of any
electronic-structure theory is then to find the minimum of
$E_{m}^{t,v,w}[\Psi]$, where $\Psi$ is a (possibly multi-determinantal) wave
function constructed from the orbital set $\\{\phi_{k}\\}$. Typically, one
needs to impose some constraints on the parametrization of the wave function
$\Psi$ to make it physical

$\displaystyle c_{k}[\Psi]=0,$ (5.4)

e.g., orthonormality of orbitals or the norm of the CI coefficients. The
_generic electronic-structure minimization problem_ is formulated as

minimize $\displaystyle E_{m}^{t,v,w}[\Psi]$ subject to $\displaystyle