Reduced Density-Matrix Approach to Strong Matter-Photon Interaction

TL;DR

This paper introduces a non-perturbative, first-principles reduced density-matrix approach to model strongly coupled matter-photon systems, enabling accurate analysis of complex electronic-photon interactions.

Contribution

It develops a novel higher-dimensional auxiliary system to map coupled fermion-boson systems into a dressed fermionic problem, facilitating the use of conventional density-matrix functional theory.

Findings

Simple density-matrix approximations are accurate across coupling regimes.

The influence of the photon field depends on electronic structure details.

Method is applicable to complex systems beyond exact calculation capabilities.

Abstract

We present a first-principles approach to electronic many-body systems strongly coupled to cavity modes in terms of matter-photon one-body reduced density matrices. The theory is fundamentally non-perturbative and thus captures not only the effects of correlated electronic systems but accounts also for strong interactions between matter and photon degrees of freedom. We do so by introducing a higher-dimensional auxiliary system that maps the coupled fermion-boson system to a dressed fermionic problem. This reformulation allows us to overcome many fundamental challenges of density-matrix theory in the context of coupled fermion-boson systems and we can employ conventional reduced density-matrix functional theory developed for purely fermionic systems. We provide results for one-dimensional model systems in real space and show that simple density-matrix approximations are accurate from the weak to the deep-strong coupling regime. This justifies the application of our method to systems that are too complex for exact calculations and we present first results, which show that the influence of the photon field depends sensitively on the details of the electronic structure.

Figures

Equations

$$\hat{H}$$

$$\Psi ({\mathbf{r}}_1{\sigma }_1,...,{\mathbf{r}}_N{\sigma }_N;p_1,...,p_M),$$

$${\hat{H}}_{ph}=\sum _{\alpha =1}^M{\omega }_{\alpha }\left({\hat{a}}_{\alpha }^{+}{\hat{a}}_{\alpha }^{}+\frac{1}{2}\right).$$

$$E_0[\Psi ]=\underset{\Psi }{\inf }⟨\Psi ⟩\hat{H}\Psi .$$

$${\Gamma }_e^{(q)}$$

$${\psi }_e^{\ast }({\mathbf{r}}_1^{\prime }{\sigma }_1,...,{\mathbf{r}}_q^{\prime }{\sigma }_q,{\mathbf{r}}_{q+1}{\sigma }_{q+1},...,{\mathbf{r}}_N{\sigma }_N){\psi }_e({\mathbf{r}}_1{\sigma }_1,...,{\mathbf{r}}_q{\sigma }_q,{\mathbf{r}}_{q+1}{\sigma }_{q+1},...,{\mathbf{r}}_N{\sigma }_N).$$

$${\Gamma }_b^{(q)}({\alpha }_1,...,{\alpha }_q;{\alpha }_1^{\prime },...{\alpha }_q^{\prime })=\frac{N_b!}{(N_b-q)!}\sum _{{\alpha }_{q+1},...,{\alpha }_{N_b}=1}^M$$

$${\psi }_b^{\ast }({\alpha }_1^{\prime },...,{\alpha }_q^{\prime },{\alpha }_{q+1},...,{\alpha }_{N_b}){\psi }_b({\alpha }_1,...,{\alpha }_q,{\alpha }_{q+1},...,{\alpha }_{N_b}).$$

$${\gamma }_b(\alpha ,\beta )=⟨\Phi ⟩{\hat{a}}_{\beta }^{+}{\hat{a}}_{\alpha }\Phi =\sum _{N_b=0}^{\infty }{N}_b\left(\sum _{{\alpha }_2,...,{\alpha }_{N_b}=1}^M{\psi }_b^{\ast }(\beta ,{\alpha }_2,...,{\alpha }_{N_b}){\psi }_b(\alpha ,{\alpha }_2,...,{\alpha }_{N_b})\right).$$

$$\left(T+V\right)[{\gamma }_e]$$

$$H_d^{(1)}[{\gamma }_e]$$

$$W[{\Gamma }_e^{(2)}]=$$

$$H_d^{(2)}[{\Gamma }_e^{(2)}]=$$

$$H_{ph}[{\gamma }_b]$$

$$H_I[{\Gamma }_{e,b}^{(3/2)}]=$$

$$=$$

$$E_0=$$

$$=$$

$$0\le$$

$$\sum _{i=1}^{\infty }{n}_i^{e/b}=N_{e/b}$$

$${\hat{\Pi }}_{\alpha }(p_{\alpha ,2},...,p_{\alpha ,N})=\sum _{i=2}^N\left(-\frac{1}{2}\frac{{\partial }^2}{\partial p_{\alpha ,i}^2}+\frac{{\omega }_{\alpha }^2}{2}{p}_{\alpha ,i}^2\right)$$

$$p_{\alpha }$$

$$\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)$$

$${\hat{H}}^{\prime }=$$

$$=$$

$$+\sum _{\alpha =1}^M\left(-\frac{1}{2}\frac{{\partial }^2}{\partial p_{\alpha }^2}+\frac{{\omega }_{\alpha }^2}{2}{p}_{\alpha }^2\right)+\sum _{\alpha =1}^M{\hat{\Pi }}_{\alpha }(p_{\alpha ,2},...,p_{\alpha ,N})$$

$$\stackrel{\text{eqref}CenterOfMass}{=}$$

$$+\frac{1}{2}\sum _{k\ne l}\left[w({\mathbf{r}}_k,{\mathbf{r}}_l)+\sum _{\alpha =1}^M\left(-\frac{{\omega }_{\alpha }}{\sqrt{N}}{q}_{\alpha ,k}{\mathbf{\lambda }}_{\alpha }\cdot {\mathbf{r}}_l-\frac{{\omega }_{\alpha }}{\sqrt{N}}{q}_{\alpha ,l}{\mathbf{\lambda }}_{\alpha }\cdot {\mathbf{r}}_k+{\mathbf{\lambda }}_{\alpha }\cdot {\mathbf{r}}_k{\mathbf{\lambda }}_{\alpha }\cdot {\mathbf{r}}_l\right)\right],$$

$${\hat{H}}^{\prime }$$

$$={\hat{T}}^{\prime }+{\hat{V}}^{\prime }+{\hat{W}}^{\prime },$$

$${\Psi }^{\prime }({\mathbf{r}}_1{\sigma }_1,...,p_{M,N})=\Psi ({\mathbf{r}}_1{\sigma }_1,...,{\mathbf{r}}_N{\sigma }_N,p_1,..,p_M)\chi (p_{1,2},...,p_{M,N}).$$

$${\Psi }^{\prime }({\mathbf{r}}_1{\sigma }_1,...,{\mathbf{r}}_N{\sigma }_N,q_{1,1},...,q_{M,N})={\Psi }^{\prime }({\mathbf{z}}_1{\sigma }_1,...,{\mathbf{z}}_N{\sigma }_N),$$

$${{\rm Y}}^{\prime }({\mathbf{z}}_1{\sigma }_1,...{\mathbf{z}}_N{\sigma }_N)\equiv {{\rm Y}}^{\prime }({\mathbf{r}}_1{\sigma }_1,...,p_M;p_{1,2},...,p_{M,N})$$

$$\underset{{{\rm Y}}^{\prime }}{\inf }⟨{{\rm Y}}^{\prime }⟩{\hat{H}}^{\prime }{{\rm Y}}^{\prime }$$

$$=⟨{\Psi }^{\prime }⟩{\hat{H}}^{\prime }{\Psi }^{\prime }.$$

$${\hat{H}}_{ph}={\hat{H}}_{ph}^{\prime }-\sum _{\alpha =1}^M{\hat{\Pi }}_{\alpha }.$$

Full text

\altaffiliation

Center for Computational Quantum Physics (CCQ), Flatiron Institute, 162 Fifth Avenue, New York NY 10010, USA

Experiments performed in the last decades (see, e.g., Refs. 1, 2) have made accessible the strong and ultra-strong111We follow the definition of the light-matter coupling regimes of Ref. 3. interaction regime between matter degrees of freedom and the quantized modes of optical cavities, which allows for the study of many new phenomena including modification of chemical reaction rates 4, 5, interacting photons in quantum non-linear media 6 or super-radiance of atoms in a photonic trap 7. At the same time it creates opportunities such as the modification of energy-transfer pathways within photosynthetic organisms 8, an increase of conductivity in organic semiconductors hybridized with the vacuum field 9 or the generation of long-distance molecular interactions that, for example, allow for energy transfer way beyond the short-range dipole-dipole mediated transfer (Förster theory) 10. All these phenomena are related by the emergence of hybrid light-matter quasi-particle states, called polaritons, that determine the properties of the respective coupled electron-photon system. The physics of these exotic states can be understood impressively well by model systems of Dicke-type 11, meaning several few-level systems coupled to some photon modes 3. Many experimentally found features of the (ultra-)strong coupling regime could be described 12, 13, 14, 15, 16 and much exciting new physics was predicted 17, 18, 19, 20, 21 using such models. This article focuses on the influence of (ultra-)strong coupling on the ground state of light-matter systems, a topic on which considerably less literature exists. Only recently, polaritonic ground states, which are believed to be fundamental for the understanding of polaritonic chemistry 1, have been started to be investigated in more detail 22, 3. However, limits and difficulties of few level-approximations have been pointed out 23, 24, 25, 26, 27, 28 and recently, new models have been used to investigate polaritonic chemistry 29, 15, 30, 31, 32, 33, 34, 35, 28, 36, 37, 38, 39, 40, 41. Still, many questions remain open, especially whether the collective (ultra-)strong coupling, predicted by the Dicke-model can actually modify ground state properties of single molecules17, 42. Another example is the ongoing discussion on the theoretical understanding of super-radiance 43, 24. These debates suggest that there is a need for new theoretical tools that treat matter and photons at the same level of theory 35, 27, 28. Up to now, standard theoretical modelling treats in detail either the photons or the matter, which becomes insufficient when the matter and photon degrees of freedom are equally important 27. We present in this paper further evidence that the impact of the light-matter interaction on the matter is far from trivial and can change from system to system.

However, the full quantum-mechanical description of just the electronic degrees of freedom is already computationally very challenging due to the exponential scaling of the wave function.222Imagine for instance a 4-electron system in one spatial dimension, coupled to one cavity mode (which corresponds to the Beryllium example, presented in Sec. Numerical Results.) Every eigenstate of the corresponding Hamiltonian would be a function that depends on 5 variables. Thus, if we wanted to calculate the exact ground state, assuming 100 grid points per coordinate (corresponding to e.g. a box size of 20 with spacing 0.2) and 8 byte per function value (double precision,) we would need $100^{5}\cdot 8\text{ Byte}\approx 75$ Gigabyte of working memory, which is close to the edge of current high performance technology. Instead, reformulations of the many-body problem in terms of reduced quantities like the electron density 44, 45 or the Green’s function 46, 47, 48 have been shown to provide accurate results for relatively low computational costs. Thus, working with reduced quantities seems to be a natural choice also for coupled light-matter systems49. Recently, quantum-electrodynamical density-functional theory (QEDFT) was introduced as an extension to pure electronic density-functional theory (DFT) 50, 51, 52, 53. First calculations showed the feasibility of QEDFT 54, 55 leading to the possibility to perform full first-principles calculations of real molecules coupled to cavity modes 56, 57. However, standard approximations in DFT (usually based on a non-interacting auxiliary Kohn-Sham system) become inaccurate if applied to strongly-correlated systems. This is very well studied in terms of the strong-correlation regime in electronic systems 58 and also observed for the existing QEDFT functionals that approximate the electron-photon interaction 56. Consequently, to study novel effects arising in the ultra-strong 59 or deep-strong coupling regime 60 of light and matter from first principles, i.e., without resorting to simplified few-level systems, one needs to develop new functionals for the combined matter-photon systems or explore alternative many-body methods.

Reduced density-matrix functional theory (RDMFT) is such a method. RDMFT is based on the electronic one-body reduced density matrix (1RDM) instead of the electronic density (as in DFT) as its basic functional variable 61. Similarly to DFT, approximations are necessary in RDMFT as it is not known how to express explicitly all expectation values of operators in terms of the 1RDM. Simple approximations within this approach 62 have proven very efficient in dealing with difficult electronic structure problems like the correct qualitative description of a dissociating molecule 63 or the prediction of the Mott-insulating phase of certain strongly-correlated solids 64. Thus, it seems worth to explore how RDMFT performs in describing also the strong interaction between molecular systems and cavity modes.

Specifically, we will discuss in this work the properties of coupled light-matter systems in a setting that resembles typical cavity experiments (see Fig. 1.) It turns out that transferring RDMFT to such systems involves overcoming additional difficulties in contrast to the DFT framework. The reason behind this are the conditions under which the corresponding reduced density matrices (RDMs) (which will be purely electronic, purely photonic and coupled) connect to the original wave function, which is crucial to construct a well-defined reduced density-matrix framework. Already for the purely electronic system the conditions under which such a connection exists (known as N-representability conditions) are not trivial 65, 66, 67. For RDMs in matter-photon systems, these conditions are entirely unexplored. Nevertheless, we manage to overcome this difficulty by mapping the original system to a higher-dimensional auxiliary system that allows for an effective description of the problem by fermions. Hence, the corresponding RDMs connect to the auxiliary wave functions under N-representability conditions of fermionic systems. Despite being fermions, the newly introduced quasi-particles will depend on electronic and bosonic degrees of freedom. This construction was recently introduced by some of the authors in Ref. 68 and used to construct a DFT scheme specifically for the strong-coupling regime of light-matter systems. In the auxiliary system, the Hamiltonian consists of only 1- and 2-body terms for the new quasi-particles, thus it has the same structure as the conventional electronic Hamiltonian for molecular systems, which allows us to apply electronic RDMFT theory without major modifications. We will present some results for model systems with the simplest known RDMFT functional, the Müller functional 62, and show that this dressed RDMFT is accurate from the weak to the strong-coupling regime. Then, we will present two examples highlighting that how matter reacts to the interaction with photons depends strongly on the system. For instance, for the same coupling strength, the repulsion between the particles can be locally suppressed in one system but enhanced in another. We finish with commenting on some open issues and challenges for future applications.

This article is structured as follows. First, we explain the physical setting of an electronic system in a cavity in Sec. Physical Setting. In Sec. Reduced Density Matrices for coupled light-matter systems, we discuss the RDMs that appear in the ground-state energy expression of the coupled light-matter system. We will explain the difficulties that are introduced by the coupling between fermions and bosons and having non-particle-conserving terms. In Sec. The “fermionization” of matter-photon systems, we introduce the higher-dimensional auxiliary system which will allow us to avoid most of the aforementioned difficulties. We show how to construct a proper RDMFT framework in this system in Sec. Dressed Reduced Density-Matrix Functional Theory, explain our numerical implementation in Sec. Numerical Implementation, and present the numerical results in Sec. Numerical Results. We finish with discussing possibilities as well as challenges of dressed RDMFT in Sec. Conclusion.

Physical Setting

To describe weak and strong matter-photon interaction, it is necessary to go beyond typical quantum-chemistry or solid-state physics theories that describe electrons in a local potential interacting via Coulomb interaction and being perturbed by a classical external electromagnetic field. Instead, we need to treat explicitly the quantum nature of light and the back-reaction between electrons and electromagnetic field excitations (photons) 35. Therefore, the framework of quantum electrodynamics (QED) needs to be employed. However, we do not want to treat QED in its full complexity but will apply some well-established approximations (see Ref. 35 for a detailed discussion.) First, we apply the Born-Oppenheimer approximation and treat the nuclei as fixed classical particles.333The extension to also include the nuclei as quantum particles is in principle straightforward by following, e.g., similar strategies like discussed in Refs. 55, 27 Second, we work in the non-relativistic limit, which for the typical energy scales of molecules and their low-energy excitations is usually sufficient. Third, we assume that the wavelength $\lambda$ of the relevant electromagnetic modes is much larger than the spatial extension $d$ of the electronic system ($\lambda\gg d$) such that the dipole approximation (here in the Coulomb gauge) is valid 69, 70. In the case of the dipole approximation, where every photon mode couples to all Fourier components of the charge current of the electronic subsystem 52, 53, an effective description with only a few modes is usually sufficient. The continuum of modes is then effectively taken into account by using instead of the bare mass the physical mass of the electrons 57, 27. Since we focus on equilibrium situations, the openness of the cavity can be neglected.

Therefore, the basic Hamiltonian that we use to describe strongly-coupled light-matter systems reads (we use atomic units throughout)

[TABLE]

Here the first two sums on the right-hand side correspond to the usual electronic many-body Hamiltonian $\hat{H}_{e}=\hat{T}+\hat{V}+\hat{W}$, used to describe the uncoupled matter system consisting of $N$ electrons in an external potential $v(\mathbf{r})$ interacting via the Coulomb repulsion $w(\mathbf{r},\mathbf{r}^{\prime})$. The third sum describes $M$ photon modes, that are characterized by their elongation $p_{\alpha}$, frequency $\omega_{\alpha}$ and polarization vectors $\bm{\lambda}_{\alpha}$. The polarization vectors include already the effective coupling strength $g_{\alpha}=|\bm{\lambda}_{\alpha}|\sqrt{\frac{\omega_{\alpha}}{2}}$35 and couple to the total dipole $\hat{\mathbf{D}}=\sum_{k=1}^{N}\mathbf{r}_{k}$ of the electronic system. The sum can be decomposed in a purely photonic part $\hat{H}_{ph}=\sum_{\alpha=1}^{M}\left(-\tfrac{1}{2}\tfrac{\partial^{2}}{\partial p_{\alpha}^{2}}+\tfrac{\omega_{\alpha}^{2}}{2}p_{\alpha}^{2}\right)$, the dipole self-interaction444Note that this term is necessary for the existence of a ground state70. $\hat{H}_{d}=\sum_{\alpha=1}^{M}\tfrac{1}{2}\left(\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}}\right)^{2}$, which we split again for later convenience in its one-body $\hat{H}_{d}^{(1)}=\sum_{\alpha=1}^{M}\sum_{k=1}^{N}\tfrac{1}{2}(\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{k})^{2}$ and two-body part $\hat{H}_{d}^{(2)}=\sum_{\alpha=1}^{M}\sum_{k\neq l}\tfrac{1}{2}(\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{k})(\bm{\lambda}_{\alpha}\cdot\mathbf{r}_{l})$, and the bilinear interaction $\hat{H}_{I}=-\sum_{\alpha=1}^{M}\omega_{\alpha}p_{\alpha}\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}}$. Some comments on this Hamiltonian are appropriate. Since we work in Coulomb gauge, the dipole-self energy term arises only for the transversal but not for the longitudinal part of the field71. However, if we assume a cavity then the Coulomb interaction is modified72. We can easily incorporate this into our framework since $w(\mathbf{r},\mathbf{r}^{\prime})$ is completely at our disposal and we do not rely on any kind of Coulomb approximation. Thus, we can also treat the influence of e.g., a plasmonic environment73.

The ground-state wave function of Eq. (1) is a function of $4N+M$ coordinates

[TABLE]

where $\sigma_{k}$ are the electronic spin degrees of freedom. The wave function $\Psi$ is as usual anti-symmetric with respect to the exchange of any two electron coordinates $\mathbf{r}_{j}\sigma_{j}\leftrightarrow\mathbf{r}_{k}\sigma_{k}$, and also depends on $M$ photon-mode displacement coordinates $p_{\alpha}$.

At this point, we want to remind the reader that there is no fundamental symmetry with respect to the exchange of two displacement coordinates $p_{\alpha}$ with $p_{\beta}$. The bosonic symmetry instead refers to the exchange of mode excitations, which are interpreted as photons in the number-state representation. To see this, we first use the ladder operators $\hat{a}_{\alpha}=\sqrt{\tfrac{\omega_{\alpha}}{2}}(p_{\alpha}-\tfrac{1}{\omega_{\alpha}}\tfrac{\partial}{\partial p_{\alpha}})$ and $\hat{a}^{+}_{\alpha}=\sqrt{\tfrac{\omega_{\alpha}}{2}}(p_{\alpha}+\tfrac{1}{\omega_{\alpha}}\tfrac{\partial}{\partial p_{\alpha}})$ to represent the sum of harmonic Hamiltonians, i.e.,

[TABLE]

Eigenstates of the individual $\hat{a}_{\alpha}^{+}\hat{a}_{\alpha}$ in the above representation are given by multiple applications of creation operators to the vacuum state $\ket{0}$, i.e., $\ket{\varphi^{n}_{\alpha}}=(\hat{a}^{+}_{\alpha})^{n}\ket{0}$.555Note that following Refs. 74, 47, we chose here explicitly a non-normalized basis $\{\ket{\varphi^{n}_{\alpha}}\}$ of the n-photon sector, with $\braket{\varphi^{n}_{\alpha}}{\varphi^{n}_{\alpha}}=n!$, which allows for a simpler definition of the bosonic 1RDM. The missing normalization factor is shifted to the resolution of identity in this basis, i.e. $\mathbb{1}=\frac{1}{N_{b}!}\sum_{\alpha_{1},...,\alpha_{N_{b}}=1}^{M}\ket{\alpha_{1},...,\alpha_{N_{b}}}\bra{\alpha_{1},...,\alpha_{N_{b}}}$, where $\ket{\alpha_{1},...,\alpha_{N_{b}}}=\hat{a}^{+}_{\alpha_{1}}\cdots\hat{a}^{+}_{\alpha_{N_{b}}}\ket{0}$ as defined later in the text. These eigenstates are connected to the displacement representation by $\varphi^{n}_{\alpha}(p_{\alpha})=\braket{p_{\alpha}}{\varphi^{n}_{\alpha}}=\braket{p_{\alpha}}{(\hat{a}^{+}_{\alpha})^{n}}{0}$. We can then express an $M$-mode eigenfunction as $\phi_{n_{1},...,n_{M}}(p_{1},...,p_{M})=\braket{p_{1}...p_{M}}{\phi_{n_{1},...,n_{M}}}=\bra{p_{1}...p_{M}}\hat{a}^{+}_{1})^{n_{1}}...$ $(\hat{a}^{+}_{M})^{n_{M}}\ket{0}$. In this form it becomes clear that a multi-mode eigenstate $\ket{\phi_{n_{1},...,n_{M}}}$ can be considered to consist of $N_{b}=n_{1}+...+n_{M}$ photons (mode excitations.) We can associate every such multi-mode eigenstate with a specific photon-number sector, i.e., the zero-photon sector is merely one-dimensional and corresponds to $\ket{\phi_{0_{1},...,0_{M}}}\equiv\ket{0}$, the single-photon sector is $M$-dimensional and corresponds to the span of $\hat{a}^{+}_{\alpha}\ket{0}\equiv\ket{\alpha}$ for all $\alpha$ and so on. For the multi-photon sectors we see due to the commutation relations of the ladder operators the bosonic exchange symmetry appearing, e.g., $\hat{a}^{+}_{\alpha_{1}}\hat{a}^{+}_{\alpha_{2}}\ket{0}=\hat{a}^{+}_{\alpha_{2}}\hat{a}^{+}_{\alpha_{1}}\ket{0}\equiv\ket{\alpha_{1},\alpha_{2}}$ for $\alpha_{1},\alpha_{2}\in\{1,...,M\}$. It is no accident that the bosonic symmetry becomes explicit in this representation since the different modes $\alpha$ determine how the photon wave functions looks like in real space (see also the discussion in App. LABEL:app:symmetry_photon.) A general photon state can therefore be represented by a sum over all photon-number sectors as $\ket{\Phi}=\sum_{n=0}^{\infty}\left(\tfrac{1}{\sqrt{n!}}\sum_{\alpha_{1},...,\alpha_{n}=1}^{M}\tilde{\Phi}(\alpha_{1},...,\alpha_{n})\ket{\alpha_{1},...,\alpha_{n}}\right)$.

Reduced Density Matrices for coupled light-matter systems

Having introduced our system of interest, we now want to discuss how to find its ground state. A ground state (if it exists) is defined as the state (possibly degenerate) that has the lowest energy expectation value

[TABLE]

This is the classical formulation of the variational principle due to Ritz and is well-defined for every Hamiltonian that is bounded from below. While well-known, Eq. (4) has the disadvantage that in practice the minimization has to be performed over an enormous configuration space that is spanned by all possible many-body wave functions. A possible reduction of computational complexity presents itself by the fact that the full wave function is usually not necessary to compute the energy expectation value but typically only reduced quantities are sufficient. Varying instead over the space of reduced objects makes the minimization simpler. For instance, in the case of an electronic many-body state $\psi_{e}(\mathbf{r}_{1}\sigma_{1},...,\mathbf{r}_{N}\sigma_{N})$ of $N$ electrons (see Tab. 1), the expectation value of a general (non-local) $q$-body operator $\hat{O}(\mathbf{r}_{1},...\mathbf{r}_{q};\mathbf{r}_{1}^{\prime},...\mathbf{r}_{q}^{\prime})$, which is given by $O=\braket{\psi_{e}}{\hat{O}\psi_{e}}$, can be determined via the electronic (spin-summed)666We define here only the spin-summed version of the qRDM, because in this work, we do not consider explicitly spin-dependent quantities. For instance, if we included magnetic fields in the Hamiltonian, the situation would change. $q$-body RDM ($q$RDM)777Note that there are different conventions in literature for the normalization of the $q$RDM. We followed Ref. 75.

[TABLE]

For instance, the well-known electronic 1RDM, that we denote in the following by $\gamma_{e}(\mathbf{r},\mathbf{r}^{\prime})=\Gamma_{e}^{(1)}(\mathbf{r};\mathbf{r}^{\prime})$, is sufficient to calculate all electronic single-particle observables such as the kinetic energy. A prominent example of a higher-order operator in the electronic case is the two-body Coulomb interaction among the $N$ electrons. To calculate its expectation value, we need to consider the diagonal of the 2RDM $\Gamma^{(2)}(\mathbf{r}_{1},\mathbf{r}_{2};\mathbf{r}_{1},\mathbf{r}_{2})$. In the chosen normalization, all RDMs satisfy the sum-rule $\int\mathrm{d}^{3q}r\Gamma^{(q)}(\mathbf{r}_{1},...,\mathbf{r}_{q};\mathbf{r}_{1},...,\mathbf{r}_{q})=\tfrac{N!}{(N-q)!}$. So for instance, the 1RDM integrates to the particle number $\int\mathrm{d}^{3}r\gamma_{e}(\mathbf{r},\mathbf{r})=N$, the 2RDM integrates to two times the number of pairs, etc. Additionally, higher and lower order RDMs are connected via $\Gamma^{(q)}_{e}(\mathbf{r}_{1},...,\mathbf{r}_{q};\mathbf{r}_{1}^{\prime},...,\mathbf{r}_{q}^{\prime})=\tfrac{1}{N-q}\int\mathrm{d}^{3}r\Gamma^{(q+1)}_{e}(\mathbf{r}_{1},...,\mathbf{r}_{q+1};\mathbf{r}_{1}^{\prime},...,\mathbf{r}_{q+1})$.

For bosons the same construction is possible. In our case, where we have a discrete set of possible single-boson states, a $N_{b}$ boson state in the (symmetrized) mode-representation $\psi_{b}(\alpha_{1},...,\alpha_{N_{b}})$ (see supporting information, Sec.  LABEL:app:symmetry_photon for more details) leads to the corresponding bosonic $q$RDM,

[TABLE]

Accordingly to the electronic case, we denote the 1RDM by $\gamma_{b}(\alpha,\beta)=\Gamma_{b}^{(1)}(\alpha;\beta)$. However, in the specific case of photons, where the number of particles is indetermined and we work with Fock-space wave functions $\ket{\Phi}$, we need to consider a Fock-space 1RDM of the form

[TABLE]

In an according manner one can define a bosonic Fock-space $q$RDM via $\Gamma^{(q)}_{b}(\alpha_{1},...,\alpha_{q};\alpha_{1}^{\prime},...,\alpha_{q}^{\prime})=\braket{\Phi}{\hat{a}^{+}_{\alpha_{1}^{\prime}}\cdots\hat{a}^{+}_{\alpha_{q}^{\prime}}\hat{a}_{\alpha_{q}}\cdots\hat{a}_{\alpha_{1}}\Phi}$.

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. (2) we can accordingly define $\Gamma_{e}^{(q)}\equiv\tfrac{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}^{+}_{\alpha_{1}}\cdots\hat{a}^{+}_{\alpha_{q}}\hat{a}_{\alpha_{q}^{\prime}}\cdots\hat{a}_{\alpha_{1}^{\prime}}\Psi}$.

In a next step, we see whether these standard ingredients of RDMFT are sufficient to express the energy expectation value of the coupled Hamiltonian of Eq. (1). For the purely electronic part, the different contributions can be expressed either explicitly by the electronic 1RDM or by the electronic 2RDM. The single-particle operators of $\hat{H}_{e}$ and the single-particle part of the dipole self-energy $\hat{H}_{d}^{(1)}$ are given in terms of the 1RDM by

[TABLE]

Here we have denoted on the left-hand side the explicit dependence of the expectation value on the 1RDM and 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 $\left(-\tfrac{1}{2}\nabla_{\mathbf{r}}^{2}+v(\mathbf{r})\right)$. 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

[TABLE]

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 find

[TABLE]

However, the bilinear coupling term is not given in a simple RDM form but becomes

[TABLE]

A new reduced quantity appears 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})$.888To see this in a simple manner 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 76. We re-express $\braket{\Psi}{\left[\left(\hat{a}^{+}_{\alpha}+\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}$. If we then define $\Gamma^{(3/2)}_{e,b}(\alpha;\mathbf{r},\mathbf{r}^{\prime})=\sum_{\sigma}\braket{\Psi}{\left[\left(\hat{a}^{+}_{\alpha}+\hat{a}_{\alpha}\right)\hat{\psi}_{e}^{\dagger}(\mathbf{r}\sigma)\hat{\psi}_{e}(\mathbf{r}^{\prime}\sigma)\right]\Psi}$ we can re-write $\braket{\Psi}{\left[\left(\hat{a}^{+}_{\alpha}+\hat{a}_{\alpha}\right)\bm{\lambda}_{\alpha}\cdot\hat{\bf D}\right]\Psi}=\int\mathrm{d}^{3}r\left(\bm{\lambda}_{\alpha}\cdot\mathbf{r}\right)\Gamma^{(3/2)}_{e,b}(\alpha;\mathbf{r},\mathbf{r})$. The bilinear interaction term therefore creates/annihilates bosons by interacting with the electronic subsystem. The $3/2$-body RDM has in general no simple connection to any $q$RDM, even if we extend the definitions to include combined matter-boson $q$RDMs.999Using the field-operator formulation, the usual $q$RDMs consist of strings of particle-number-conserving combinations of electron and boson operators. Integrating/summing out a number-non-conserving part of it does not lead to a simple relation to half-body RDMs in general. One obvious reason is that $q$RDMs conserve particle numbers, while half-body RDMs do not. Take, for instance, the Fock-space $\gamma_{b}(\alpha,\beta)=\braket{\Phi}{\hat{a}^{+}_{\beta}\hat{a}_{\alpha}\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}^{+}_{\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.

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

[TABLE]

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{r}_{1}\sigma_{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 RDMFT, and for electronic systems the properties of RDMs have been studied for more than 50 years 77. However, this seeming reduction of complexity is deceptive. 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 a physical wave function. This is indicated in Eq. (Reduced Density Matrices for coupled light-matter systems), 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. For systems with fixed particle numbers it is in principle known how to restrict the set of trial RDMs to physical ones. The corresponding restrictions are called N-representability conditions 65, 66, 67. However, only for the 1RDM of ensembles (fermionic or bosonic) the conditions are simple. In this case, by diagonalizing the 1RDM in its eigenbasis $\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, the conditions are

[TABLE]

for fermions and bosons, respectively. If the particle number $N_{e/b}$ of one species of the system is conserved, the respective sum-rule

[TABLE]

becomes a second part of the N-representability conditions. Consequently, to define a proper RDM framework for coupled electron-boson problems, one would need to know the corresponding constraints that connect the wave function with all the necessary RDMs. A glance at the history of an important example, the search for the N-representability conditions of the electronic 2RDM, suggests that finding similar conditions for the novel half-body RDMs together with connections between the fermionic and bosonic $q$RDMs is a very challenging task. The electronic-2RDM problem was proposed in 1960 78 and it took until 2012, to understand how to make the conditions explicit 67.

Although the connection of the different RDMs in coupled fermion-boson systems is a very interesting subject, and recent results for a grand-canonical formulation of fermions or bosons suggest that also a combined formulation is feasible 74, we will follow an alternative route in this work. We reformulate the problem in such a way that we “fermionize” the coupled fermion-boson problem, where we can apply the known conditions of the fermionic problem.

The “fermionization” of matter-photon systems

In this section, we explain in detail how a system described by the Hamiltonian (1), is mapped to an auxiliary space such that the coupled matter-light degrees of freedom can be modelled with new particles that are fermionic. We call them dressed or polaritonic particles, because they depend on electronic and photonic coordinates.101010Note that the use of a dressed particle picture allows to also describe Landau polaritons, as shown recently in Ref. 79. Thus, also such systems can in principle be considered with the presented approach. This “fermionization” procedure was introduced in a recent work by Nielsen et al. 68 and can be divided in three steps. First, we introduce for each mode auxiliary extra dimensions $(p_{\alpha,2},...,p_{\alpha,N})$, where the number of these extra dimensions depends on the number of electrons $N$. We therefore embed the physical configuration space in a higher-dimensional space, i.e., we now consider wave functions depending on $(\mathbf{r}_{1}\sigma_{1},...,\mathbf{r}_{N}\sigma_{N},p_{1},..,p_{M},p_{1,2},...,p_{1,N},...,p_{M,2},...,p_{M,N})$. Second, we add for every photon mode $\alpha$ a linear operator

[TABLE]

to the physical Hamiltonian of Eq. (1). This auxiliary Hamiltonian is a sum of quantum harmonic oscillators with respect to the auxiliary coordinates. The resulting Hamiltonian in the extended configuration space is $\hat{H}^{\prime}=\hat{H}+\sum_{\alpha=1}^{M}\hat{\Pi}_{\alpha}$ (we denote all quantities in the auxiliary space with a prime symbol.) Here we see that the auxiliary degrees of freedom do not mix with the physical ones. This will allow in a very simple manner to embed but also to reconstruct the physical wave function. In the third step, we perform an orthogonal variable transformation of the photonic plus auxiliary coordinates to new coordinates $(q_{\alpha,1},...,q_{\alpha,N})$ such that

[TABLE]

This whole procedure can be viewed as the inverse of a center-of-mass coordinate transformation 80. In total, we find the auxiliary Hamiltonian in the higher-dimensional configuration space given as

[TABLE]

where we inserted the definition of the dipole operator, $\hat{\mathbf{D}}=\sum_{k=1}^{N}\mathbf{r}_{k}$, 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 auxiliary photon coordinates $\mathbf{z}=(\mathbf{r},q_{1},...,q_{M})$, we can rewrite the above Hamiltonian as

[TABLE]

where we introduced the dressed Laplacian $\Delta^{\prime}=\sum_{i=1}^{3}\tfrac{\partial^{2}}{\partial r_{i}^{2}}+\sum_{\alpha=1}^{M}\tfrac{\partial^{2}}{\partial q_{\alpha}^{2}}$, the dressed local potential $v^{\prime}(\mathbf{z})=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}\right.$ $\left.+\tfrac{1}{2}(\bm{\lambda}_{\alpha}\!\!\cdot\!\mathbf{r})^{2}\right]$ and the dressed interaction kernel $w^{\prime}(\mathbf{z},\mathbf{z}^{\prime})=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}\right.$ $\left.-\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]$. Note that the choice of the linear auxiliary operator (Eq. (10)) and the coordinate transformation (Eq. (The “fermionization” of matter-photon systems)) have a certain freedom. The operator of Eq. (10) must not contain physical coordinates, such that the physical system cannot be influenced by this auxiliary operator. Further, it must together with the transformation give rise to polaritonic 1- and 2-body terms, as shown in Eq. (12). These requirements are met using an orthonormal transformation111111Note that there are many different orthonormal transformations, but the exact choice is not important for the formalism. It only needs to include the first line of Eq. (The “fermionization” of matter-photon systems). One specific example of such a transformation is68 $p_{\alpha,k}=\frac{1}{\sqrt{k^{2}-k}}(q_{1}+...+q_{k-1}-(k-1)q_{k})$ for $2\leq k\leq N$ alongside the first line of Eq. (The “fermionization” of matter-photon systems). together with the harmonic oscillators of Eq. (10). Since the operator of Eq. (10) only acts in the auxiliary coordinates, the normalized physical solution $\Psi$ of the original (time-independent) Schrödinger equation $E_{0}\Psi=\hat{H}\Psi$ in the standard configuration space is connected to a new physical solution of the auxiliary Hamiltonian $\hat{H}^{\prime}$ in a very simple manner, i.e.,

[TABLE]

Here the normalized solution $\Psi^{\prime}$ of the auxiliary Schrödinger equation $E_{0}^{\prime}\Psi^{\prime}=\hat{H}^{\prime}\Psi^{\prime}$ is found with $\chi$ being the ground state of $\sum_{\alpha}\hat{\Pi}_{\alpha}$. The ground state $\chi$ is merely a tensor product of individual harmonic-oscillator ground states and therefore exchanging $p_{\alpha,i}$ with $p_{\alpha,j}$ does not change the total wave function $\Psi^{\prime}$. If we rewrite this wave function in the new coordinates

[TABLE]

and due to the fact that we constructed $\chi$ to be symmetric with respect to the exchange of $q_{\alpha,k}$ and $q_{\alpha,l}$,121212Note that the auxiliary ground state of mode $\alpha$ is given as a Gaussian with respect to $\sum_{i=2}^{N}p_{\alpha,i}^{2}=\sum_{k=1}^{N}q_{\alpha,k}^{2}-\tfrac{1}{N}(\sum_{k=1}^{N}q_{\alpha,k})^{2}$, where we used that $p=\frac{1}{\sqrt{N}}\sum_{k=1}^{N}q_{\alpha,k}$. we realize that $\Psi^{\prime}$ is anti-symmetric with respect to the exchange of $(\mathbf{z}_{k}\sigma_{k})$ and $(\mathbf{z}_{l}\sigma_{l})$. Thus $\Psi^{\prime}$ is fermionic with respect to the polaritonic coordinates $(\mathbf{z}\sigma)$ and can be represented by a sum of Slater determinants of $(3+M)$-dimensional polaritonic orbitals $\varphi(\mathbf{z}\sigma)$. This makes the application of usual fermionic many-body methods possible and we can rely on fermionic $N$-representability conditions. However, besides the extra dimensions and the new fermionic exchange symmetry, the physical wave functions of the dressed auxiliary space also have a further $q_{\alpha,k}\leftrightarrow q_{\alpha,l}$ exchange symmetry. Simple approximations based on single polaritonic Slater determinants will violate this extra symmetry. We will remark on such violations when we introduce dressed RDMFT in the next section. To further see that indeed the constructed $\Psi^{\prime}$ is a minimal-energy state in the extended space with the appropriate symmetries, we first point out that for any trial wave function

[TABLE]

the $q$-exchange symmetry implies that the fermionic symmetry is in the $(\mathbf{r}\sigma)$ coordinates. Then it holds since $\hat{H}$ only acts on $(\mathbf{r}_{1},...,p_{M})$ and $\sum_{\alpha=1}^{M}\hat{\Pi}_{\alpha}$ only acts on $p_{1,2},...,p_{M,N}$ that

[TABLE]

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}$131313Note that $\hat{H}^{\prime}$ also has many eigenstates that are not of the form $\Psi^{\prime}=\Psi\chi$, with $\Psi$ anti-symmetric under exchange of $\mathbf{r}_{k}\sigma_{k}\leftrightarrow\mathbf{r}_{l}\sigma_{l}$ and $\chi$ symmetric under exchange of $q_{\alpha,k}\leftrightarrow q_{\alpha,l}$. In general, we thus have to enforce these properties to only retain the eigenstates of this form. However, for the ground state it is sufficient to enforce the symmetry under exchange of $q_{\alpha,k}\leftrightarrow q_{\alpha,l}$ together with the anti-symmetry under exchange of $\mathbf{z}_{k}\sigma_{k}\leftrightarrow\mathbf{z}_{l}\sigma_{l}$. 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 (The “fermionization” of matter-photon systems) 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 non-trivial 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 elongation of $p_{\alpha}=\frac{1}{\sqrt{N}}(q_{\alpha,1}+...+q_{\alpha,N})$ and its conjugate $\frac{\partial}{\partial p_{\alpha}}=\frac{1}{\sqrt{N}}(\frac{\partial}{\partial q_{\alpha,1}}+...+\frac{\partial}{\partial q_{\alpha,N}})$, 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 (The “fermionization” of matter-photon systems),

[TABLE]

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 $\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}$. 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{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

[TABLE]

Dressed Reduced Density-Matrix Functional Theory

From the observations of the previous chapter, it is clear that a simple dressed Kohn-Sham approximation can capture matter-photon correlations that are hard to capture with standard approximations of Kohn-Sham QEDFT. But from the experience with purely electronic DFT, we expect that for ultra- and deep-strong coupling situations the simple dressed approximations also become less reliable. Instead of developing more advanced approximations for a dressed Kohn-Sham approach, we propose to follow another route in this paper. Similar to electronic-structure theory, where RDMFT becomes a reasonable alternative to DFT methods when strong correlations become important 63, 64, 81, we present a dressed RDMFT approach to capture ultra- and deep-strong electron-photon coupling.

Let us therefore analyze the structure of $\hat{H}^{\prime}$, given in Eq. (13). It consists of only polaritonic one-body terms $\hat{h}^{(1)}(\mathbf{z})=-\tfrac{1}{2}\Delta+v^{\prime}(\mathbf{z})$, and two-body terms $\hat{h}^{(2)}(\mathbf{z},\mathbf{z}^{\prime})=w^{\prime}(\mathbf{z},\mathbf{z}^{\prime})$. It commutes with the polaritonic particle-number operator $\hat{N}^{\prime}=\int\mathrm{d}^{3+M}z\,\hat{n}(\mathbf{z})$, where we used the definition of the polaritonic local density operator $\hat{n}(\mathbf{z})=\sum_{i=1}^{N}\delta^{3+M}(\mathbf{z}-\mathbf{z}_{i})$. This means that the auxiliary system has a constant polaritonic particle number $N$. Additionally, the physical wave function of the dressed system $\Psi^{\prime}(\mathbf{z}_{1}\sigma_{1},\dots,\mathbf{z}_{N}\sigma_{N})$ is per construction anti-symmetric. This allows for the definition of a dressed (spin-summed) 1RDM

[TABLE]

in accordance to Eq. (5). By construction, this auxiliary density matrix reduces to the physical electron density matrix via $\gamma_{e}(\mathbf{r},\mathbf{r}^{\prime})=\int\mathrm{d}^{M}q\,\gamma(\mathbf{r},q_{1},..,q_{M};\mathbf{r}^{\prime},q_{1},..,q_{M})$. Furthermore, we introduce the (spin-summed) dressed 2RDM $\Gamma^{(2)}(\mathbf{z}_{1},\mathbf{z}_{2};\mathbf{z}_{1}^{\prime},\mathbf{z}_{2}^{\prime})=N(N-1)\sum_{\sigma_{1},...,\sigma_{N}}\int\mathrm{d}^{(3+M)(N-2)}z$ ${\Psi^{\prime}}^{*}(\mathbf{z}_{1}^{\prime}\sigma_{1},$ $\mathbf{z}_{2}^{\prime}\sigma_{2},\mathbf{z}_{3}\sigma_{3},..,\mathbf{z}_{N}\sigma_{N})$ ${\Psi^{\prime}}(\mathbf{z}_{1}\sigma_{1},\mathbf{z}_{2}\sigma_{2},$ $\mathbf{z}_{3}\sigma_{3},..,\mathbf{z}_{N}\sigma_{N})$. These dressed RDMs allow for expressing the energy expectation value of the dressed system by

[TABLE]

Thus, we can define the variational principle for the ground state only with respect to well-defined reduced quantities,

[TABLE]

To perform this minimization, we need to constrain the configuration space to the physical dressed RDMs that connect to an anti-symmetric wave function with the extra $q$-exchange symmetry by testing the appropriate N-representability conditions of the dressed 2RDM and the dressed 1RDM. Besides the by now well-known conditions for the fermionic 2RDM 67 and the fermionic 1RDM 65 we would in principle get further conditions to ensure the extra exchange symmetry. However, already for the usual electronic 2RDM the number of conditions grows exponentially with the number of particles, and it is out of the scope of this work to discuss possible approximations. The interested reader is referred to, e.g., Ref. 82. Instead, we want to stick to the dressed 1RDM $\gamma$ and approximate the 2-body part as a functional of the $\gamma$. The mathematical justification of RDMFT is given by Gilbert’s theorem 61, which is a generalization of the Hohenberg-Kohn theorem of DFT 83. More specifically, Gilbert proves that the ground state energy of any Hamiltonian with only 1-body and 2-body terms is a unique functional of its 1RDM. Following this idea, we will express the ground-state energy of the dressed system as a partly unknown functional $F^{\prime}$ of only the system’s dressed 1RDM

[TABLE]

For this minimization, we need a functional of the diagonal of the dressed 2RDM in terms of the dressed 1RDM as well as adhering to the corresponding $N$-representability conditions when varying over $\gamma$. We also see now the advantage of the dressed RDMFT approach, which avoids the original variation over all wave functions as well as a variation over many different RDMs as shown in Eq. $\eqref{eq:var_parinciple_RDM}$. Instead, we only need the dressed 1RDM, which has a comparatively simple connection to fermionic wave functions (at least when we vary over ensembles, i.e., Eq.s (Reduced Density Matrices for coupled light-matter systems,9).) The price we pay for this is twofold: First, we have a new symmetry that will most likely lead to extra $N$-representability conditions. Second, we need to increase the dimension of the natural orbitals by one for every photon mode. However, to capture the main physics of usual cavity experiments, often one effective mode is enough. Computations with four-dimensional dressed orbitals are numerically feasible. It is specifically such settings, where we envision a dressed RDMFT to be a reasonable alternative to other ab-initio approaches to cavity QED68, 52, 55, 40, 30, 41.

Another advantage of RDMFT in general is the direct access to all one-body observables. This transfers also to dressed RDMFT. The calculation of expectation values of purely electronic one-body observables is trivial with the knowledge of the dressed 1RDM, but also photonic one-body (and half-body) observables can be calculated, using the connection formula shown in the end of Sec. The “fermionization” of matter-photon systems. Thus, we are able to calculate very interesting properties of the cavity photons like the mode occupation or quantum fluctuations of the electric and magnetic field.

To see whether our approach is practical and accurate, we perform first simple calculations for coupled matter-photon systems. We will make the following pragmatic approximations: We only enforce the fermionic ensemble $N$-representability conditions141414This approximation is similarly employed in electronic RDMFT. For details we refer to, e.g., Ref. 84. and ignore presently the extra exchange symmetry between the q-coordinates,151515In all the numerical examples that we studied, this “fermion polariton approximation” was very accurate. One can find more details in the supporting information. and we employ simple approximations to the unknown part $W^{\prime}[\gamma]$ that have been developed for the electronic case. To do so, we further, similarly to the electronic case, decompose $W^{\prime}[\gamma]=E_{H}[\gamma]+E_{xc}[\gamma]$ into a classical Hartree part $E_{H}[\gamma]=\frac{1}{2}\int\int\mathrm{d}^{3+M}z\mathrm{d}^{3+M}z^{\prime}\gamma(\mathbf{z},\mathbf{z})\gamma(\mathbf{z}^{\prime},\mathbf{z}^{\prime})w^{\prime}(\mathbf{z},\mathbf{z}^{\prime})$ and an unknown exchange-correlation part $E_{xc}[\gamma]$. Almost all known functionals $E_{xc}[\gamma]$ are expressed in terms of the eigenbasis and eigenvalues of the 1RDM. In our case the dressed natural orbitals $\phi_{i}(\mathbf{z})$ and occupation numbers $n_{i}$ are found by solving $\int\gamma(\mathbf{z},\mathbf{z}^{\prime})\phi_{i}(\mathbf{z}^{\prime})\mathrm{d}^{3+M}z^{\prime}=n_{i}\phi_{i}(\mathbf{z})$. The simplest approximation is to only retain the fermionic exchange symmetry and employ the Hartree-Fock (HF) functional

[TABLE]

As the HF functional depends linearly on the natural occupation numbers, any kind of minimization will lead to the single-Slater-determinant HF ground state (which corresponds to occupations of 1 and 0.) 85 We call this approximation dressed Hartree-Fock (dressed HF.) We can go beyond the single Slater determinant in dressed RDMFT, if we employ a non-linear occupation-number dependence in the exchange-correlation functional. We here employ the Müller functional 62

[TABLE]

which has been re-derived by Bjuise and Baerends 86. The Müller functional has been studied for many physical systems 63, 86 and gives a qualitatively reasonable description of electronic ground states. Additionally, it has many advantageous mathematical properties 62, 87. A thorough discussion of different functionals goes beyond the scope of this work, and we only want to remark that a variety of functionals were proposed after $E_{M}[\gamma]$ and it is likely to have even better agreement with the exact solution by choosing more elaborate functionals.

Numerical Implementation

Besides the fact that we can re-use many ideas from electronic RDMFT, a further advantage of the dressed reformulation is that we can also re-use most of the numerical techniques developed for quantum chemistry and materials science. For instance, to determine the dressed orbitals we merely need to be able to solve higher-dimensional static Schrödinger-type equations. We only have to change the usual electronic potential $v$ to its dressed counter part $v^{\prime}$. This, together with a change of the electronic Coulomb interaction $w$ to its dressed counter part $w^{\prime}$, already allows to perform dressed HF calculations, at least under the approximation of violating the additional symmetry constraint, discussed in the previous section. If the code one uses is furthermore able to perform RDMFT minimizations, it is straightforward to extend the implementation to also solve coupled electron-photon problems via dressed RDMFT from first principles. We have done so with the electronic-structure code Octopus 88 and the implementation will be made available with the upcoming release.

Specifically, we rewrite the approximated energy functional in the natural orbital basis as

[TABLE]

We use this form to minimize the energy functional by varying the natural orbitals as well as the natural occupation numbers. To impose fermionic ensemble $N$-representability, we first represent the occupation numbers as the squared sine of auxiliary angles, i.e. $0\leq n_{i}=2\sin^{2}(\alpha_{i})\leq 2$, to satisfy Eq. (Reduced Density Matrices for coupled light-matter systems).161616Note that the $n_{i}$ are bounded by 2 because we employed a spin-summed formulation. If we considered natural spin-orbitals instead, the upper bound would be 1. The second part of the conditions (Eq. (9)), i.e., $\sum_{i=1}n_{i}=N$, as well as the orthonormality of the dressed natural orbitals, i.e., $\int\mathrm{d}^{3+M}z\phi^{*}_{i}(\mathbf{z})\phi_{j}(\mathbf{z})=\delta_{ij}$, are imposed via Lagrange multipliers as, e.g., explained in Ref. 88. We have available two different orbital-optimization methods, a conjugate-gradient algorithm171717We used this method only for some benchmark calcualtions as it is not yet optimized for the needs of RDMFT. and an alternative method that was introduced by Piris et al. in Ref. 89. The latter expresses the $\phi_{i}$ in a basis set and can use this representation to considerably speed up calculations in comparison to the conjugate-gradient algorithm. It was used for all results presented in this paper. However, it is not trivial to converge such calculations in practice and we developed a protocol to obtain properly converged results. The interested reader is referred to Sec. LABEL:APP:protocol of the supporting information.

Numerical Results

In the following, we present some examples of few electron systems in one spatial dimension. In the first part of this section, we validate our method by comparing to exact solutions of simple atomic and molecular systems. Then we show that our method also provides reasonable results for a more complex system. We finish the section with two examples that illustrate how our method can describe and uncover non-trivial changes of the matter due to its coupling to photons. We first present the dissociation of a molecule as an example for a chemical reaction. Despite the dipole approximated coupling the cavity photons affect the ground state locally differently. The changes also have a non-trivial dependence on the inter-atomic distance. As this system can be solved exactly, we can again validate that dressed RDMFT reproduces these intricate effects accurately. Finally, we show also that the ground state modifications of atomic systems, that have a very similar density profile outside of the cavity, are localized and depend strongly on the detailed electronic structure. These results highlight how cavity photons can at once locally enhance and suppress electronic repulsion and modify the electronic structure considerably.

The different systems are described by a local potential $v(x)$ and coupled to one photon mode. We transfer the systems in the dressed basis, that leads to a dressed local potential $v^{\prime}(x,q)=v(x)+\tfrac{1}{2}(\lambda x)^{2}+\tfrac{\omega^{2}}{2}q^{2}-\tfrac{\omega}{\sqrt{2}}q(\lambda x)$. Specifically, we consider a one-dimensional model of a helium atom ($He$), i.e., $v_{He}(x)=-\frac{2}{\sqrt{x^{2}+1}}$, a one-dimensional model of a hydrogen molecule ($H_{2}$), i.e., $v_{H_{2}}(x)=-\frac{1}{\sqrt{(x-d)^{2}+1}}-\frac{1}{\sqrt{(x+d)^{2}+1}}$, first at its equilibrium position $d=d_{eq}=1.628$ a.u., later with varying $d$, and a one-dimensional model of a Beryllium atom $Be$, i.e., $v_{Be}(x)=-\frac{4}{\sqrt{x^{2}+\epsilon^{2}}}$. For the latter, we consider a smaller softening parameter $\epsilon=0.5$ to make sure that all electrons are properly bound. We use the soft Coulomb interaction $w(x,x^{\prime})=1/\sqrt{|x-x^{\prime}|^{2}+1}$ 90, 91 for all test-systems. For the two-electron examples, we set the photon frequency in resonance with the lowest excitations of the respective “bare” systems, so outside of the cavity. For that we calculate the ground and first excited state of each system with the exact solver and find the corresponding excitation frequencies $\omega_{He}=0.5535$ a.u. and $\omega_{H2}=0.4194$ a.u. For $Be$ instead, we choose $\omega_{Be}=3.0$, which is not a resonance of the $Be$ atom, for rather numerical than physical reasons. As resonance is not an important feature for ground state calculations, different choices of $\omega$ do not crucially change the physics of the investigated system. This is in contrast to the excited states and the ensuing Rabi splitting 27. Instead, the chosen $\omega_{Be}$ considerably enhances the numerical stability of the calculations, which has the following reason: At the current state of our implementation, we need to make use of a basis set that we generate by a preliminary calculation. To generate a basis that captures electronic and photonic parts of the system equally well, we need to make sure that the energy scales of both degrees of freedom are similar. This can be controlled easiest by varying $\omega$. We want to stress that this basis-set issue is not a fundamental problem of the dressed orbital approach. On the one hand, we plan to control the photonic basis directly and on the other hand, we are working on optimizing an alternative conjugate-gradient routine that does not use a specific basis. Details can be found in Sec. LABEL:APP:dHF_basis_conv of the supporting information.

We start with the discussion of the $(N\!=\!2)$-particles examples. In this setting the dressed auxiliary system is 4-dimensional (2 particles with 2 coordinates each), which is still small enough to be solved exactly in a 4-dimensional discretized simulation box, so that we can compare dressed RDMFT (with the Müller functional of Eq. (19)), dressed HF (see Eq. (18)) and the exact solutions. We used the box lengths of $L_{x}=L_{q}=16$ a.u. and spacings of $dx=dq=0.14$ a.u. to model the electronic x and photonic coordinates q of the two dressed particles in the exact routine.181818We want to mention that the box length is not entirely converged with these parameters. In a (numerically very expensive) benchmark calculation, we observed a further decrease of energy with larger boxes (the calculations with respect to the spacing are converged,) but the changes in energy and density are only of the order of $10^{-5}$ or less. All the following results require a maximal precision of the order of $10^{-2}$ in energy as well as in the density and thus we can safely use the given parameters. Details can be found in Sec. LABEL:APP:ValidationMB of the supporting information. For dressed RDMFT and dressed HF instead, we needed to consider 2-dimensional simulation boxes for every dressed orbital and we set $L_{x}=L_{q}=20$ a.u. and $dx=dq=0.1$ a.u. We obtained converged results for $\mathcal{M}=41$ ($\mathcal{M}=71$) natural orbitals for $He$($H_{2}$). Details about how we determined these parameters can be found in the supporting information.

We first show (see Fig. 2) the deviations of the ground state energies for dressed RMDFT and dressed HF from the exact dressed calculation as a function of the dimensionless relation between effective coupling strength and photon frequency $g/\omega$191919This quantity is typically used as a measure for the strength of the light-matter interaction, see e.g. Ref. 28. for $He$ and $H_{2}$, respectively. We thereby go from weak to deep-strong coupling with $g/\omega=1$3. The deep-strong coupling regime has been reached in different systems like for instance for Landau polaritons60. For (organic) molecules the highest reported coupling strengths are in the ultra-strong regime of $g/\omega\approx 0.4$3. We see that while dressed HF deviates strongly for large couplings, dressed RDMFT remains very accurate over the whole range of coupling strength. Still, a more severe test of the accuracy of our method is if instead of merely energies, we compare spatially resolved quantities like the ground-state density $\rho(x,q)\equiv\gamma(x,q;x,q)$. To simplify this discussion, we separate the electronic and photonic parts of the two-dimensional density by integration, i.e., $\rho(x)=\int\mathrm{d}q\,\rho(x,q)$ and $\rho(q)=\int\mathrm{d}x\,\rho(x,q)$. The exact reference solutions show that with increasing $g/\omega$ the electronic part of the density becomes more localized, while the photonic part becomes broadened (for details of the effects of matter-photon coupling see the discussion following Fig. 10.) This behavior is captured qualitatively with dressed HF as well as with dressed RDMFT. The latter performs for the electronic density considerably better over the whole range of coupling strength, whereas for the photonic densities both levels of theory deviate in a similar way from the exact result. This is shown for $g/\omega=0.1$ in Fig. 3 for both test systems. Looking at the electronic densities, we can observe a feature that the ground state energy does not reveal. In some cases the effects of the two approximations are contrary to each other as we can see in the $He$ case. Here, the dressed RDMFT electronic density is more localized around the center of charge than the exact reference and the electronic density of dressed HF less. In other cases instead, both theories over-localize $\rho(x)$ (here visible for $H_{2}$.)

An even more stringent test of the accuracy of the dressed RDMFT approach is to compare the dressed 1RDMs. The essential ingredients of the dressed 1RDMs are their natural orbitals $\phi_{i}(x,q)$. Again, we separate electronic and photonic contributions and show their reduced electronic density $\rho_{i}(x)=\int\mathrm{d}q\,|\phi_{i}(x,q)|^{2}$. Fig. 4 depicts the first three dressed natural orbital densities of dressed RDMFT in comparison with the exact ones for both test systems. While it holds that for both systems, the lowest natural orbital density of the dressed RDMFT approximation is almost the same as the exact one, and the second natural orbital densities are only slightly different, the third natural orbital densities of $H_{2}$ differ even qualitatively. For $He$, similar strong deviations are visible for the fourth natural orbital. However, as long as such strong deviations only occur for natural orbitals with small natural occupation numbers, like in these cases ($H_{2}:n_{1}=1.878,n_{2}=0.102,n_{3}=0.015$, He: $n_{1}=1.978,n_{2}=0.020,n_{3}=0.001$), their (inaccurate) contribution to the density and total energy remains small.

To complete the picture, we also look at the photonic natural orbital densities, $\rho_{i}(q)=\int\mathrm{d}x\,|\phi_{i}(x,q)|^{2}$, the first 3 of which are plotted in Fig. 5, for $He$ and $H_{2}$. Here, the dressed RDMFT results even agree better with the exact solution than their electronic counterparts. Apparently, dressed RDMFT captures the photonic properties of the tested systems very accurately for the ultra-strong coupling regime. The accuracy drops with increasing $g/\omega$.

As an example for a photonic observable, we show in Fig. 6 the mode occupation $N_{ph}$ as a function of the coupling strength $g/\omega$ that we calculated by using Eq. (14), i.e. $N_{ph}=\frac{E_{ph}}{\omega}-\frac{N}{2}$, with the photon mode energy $E_{ph}=\sum_{i=1}^{\mathcal{M}}n_{i}\int\mathrm{d}x\mathrm{d}q\,\phi^{*}_{i}(x,q)\left(-\frac{1}{2}\frac{\mathrm{d}^{2}}{\mathrm{d}q^{2}}+\frac{w^{2}}{2}q^{2}\right)\phi_{i}(x,q)$. From weak to the beginning of the ultra-strong coupling regime ($g/\omega\approx 0.1$,) both dressed HF and dressed RDMFT capture $N_{ph}$ well. For very large coupling strengths, the deviations to the exact mode occupation becomes sizeable. This might sound counter-intuitive, as the photonic density is described comparatively well. The reason is that the photon occupation, in contrast to the density, is mainly determined by the second and third natural orbital, because the first natural orbital resembles a photonic ground state with occupation number zero in the studied cases. Dressed HF does not consider a second orbital (the first instead is doubly occupied) and thus cannot capture the effect. And for dressed RDMFT, the error in the second and third natural orbital is much larger than in the first (see Fig. 5.) However, it is probable that this can be improved by better functionals.

By comparing to the exact solution, we showed that the dressed-orbital construction seems to be a reasonable starting point for an approximate description of both the electronic and the photonic part of coupled matter-photon systems. Thus, we can now go one step further and present results for a many-body systems that cannot easily be solved exactly: the one-dimensional $Be$-atom in a cavity. In Fig. 7, we see the total energy as a functional of the coupling strength $g/\omega$ for dressed HF and dressed RDMFT, respectively. Like in the two-electron systems, the deviation between both curves increases for larger $g/\omega$ and as expected the dressed RDMFT energies are lower than the dressed HF results. Analyzing the ground-state densities, we see a similar trend as in the 2-particle systems. With increasing $g/\omega$, the electronic (photonic) part of the density becomes more (less) localized, though the details differ as we show in the last part of this section (see Fig. 11 and the corresponding part in the main text.) Comparing dressed RDMFT with dressed HF, we observe that the variation of the electronic (photonic) density with increasing coupling strength is less (more) prounounced for dressed RDMFT, as Fig. 8 shows. We conclude the survey of $Be$ with the mode occupation under variation of the coupling strength (see Fig. 9.) We see that the value of $g/\omega\approx 0.5$ separates two regions. For $g/\omega<0.5$ dressed RDMFT finds a larger mode occupation than dressed HF and for $g/\omega>0.5$ instead the dressed HF mode is stronger occupied. We found similar behavior also for the 2-particle systems, although the boundary between the 2 regions was considerably different there ($He:g/\omega\approx 0.8,H_{2}:g/\omega\approx 0.1$, see Fig. 6.)

We conclude this section with two examples, for which the light-matter interaction changes the bare systems non-trivially, depending not only on the coupling strength but also on the details of the electronic structure. We start with the dissociation of $H_{2}$ as an example for a chemical reaction, where we use $v_{H_{2}}(x)$ with different $d$. In Fig. 10, we see the density of two H-atoms under variation of the distance $d$ with and without the (strong) coupling to the cavity. We see that the influence of the cavity mode strongly depends on the exact electronic structure. The interaction with the cavity mode can locally reduce or enhance the electronic repulsion due to the Coulomb interaction, where the exact interplay between both effects depends on the interatomic distance. Thus, we can observe a number of different effects like pure localization of the density towards the center of charge ($d=1$) or localization combined with a local enhancement of repulsion such that the density deviations exhibit a double peak structure ($d=2$.) The local enhancement of electronic repulsion can grow so strong that the density at the center of charge is reduced but at the same time the density maxima shift closer to each other, which is an effective suppression of electronic repulsion ($d=3$.) This interplay is reflected in the natural orbitals and occupation numbers. The coupling shifts a considerable amount of occupation from the first natural orbital to the second and third one. The contribution to the total density of the former (latter) has the character of enhanced (suppressed) electron repulsion. To show the potential of these effects, we present calculations in the deep-strong coupling regime with $g/\omega=1.0$, where the effects reach the order of 10 % of the unperturbed density, which is enormous. For smaller coupling strengths of the order of $g/\omega=0.1$, these effects are as diverse, but naturally smaller with density deformations of the order of $10^{-3}$. However, as every observable depends on the density, such deviations are significant. Remarkably, dressed RDMFT reproduces the effects accurately.

In the second example, we compare the behavior of the $He$ and $Be$ atoms under the influence of the cavity. Though the shapes of the electronic density of the two bare systems are very similar (see insets in Fig. 11,) they behave very differently under the influence of the cavity, which can be seen in Fig. 11. The electronic density of $He$ is pushed towards its center of charge with increasing coupling strength, which can be understood as a suppression of the electronic repulsion induced by the Coulomb interaction. As $He$ can be understood very well with only one orbital this is to be expected.202020For $g/\omega=0.8$ we still observe $n_{1}=1.85$. However, it should be noted that the good agreement of dressed RDMFT with the exact calculation in comparison to dressed HF is exactly because of the contribution of the second natural orbital, that is (still considerably) occupied with $n_{2}=0.14$. Things change for $Be$, where we have several dominant orbitals. With increasing coupling strength, we see like in the dissociation example a subtle interplay between suppression and local enhancement of the electronic repulsion, that depends on the coupling strength. Thus for the same coupling strength, we can observe opposite ($g/\omega=0.1$ and $0.4$ in the plot) but also similar effects ($g/\omega=0.8$ in the plot) in two systems that have almost the same “bare” density shape. Like in the dissociation example, this intricate behavior can be understood by the interplay of the different natural orbitals contributing to the electronic density. In this particular case, the main physics happens in the second and third natural orbital, where the former (with a double-peak structure) loses a considerable amount of occupation to the latter (with a triple peak structure) with increasing coupling strength.

These (seemingly simple) examples show how subtle details of the electronic structure influence the changes induced by the coupling to photons. We see a non-trivial interplay between local suppression and enhancement of the Coulomb induced repulsion between the particles. This is reflected in the natural orbitals and occupation numbers of the light-matter system and thus influences all possible observables. Such small changes have been shown to theoretically affect chemical properties and reactions strongly, which are determined by an intricate interplay between Coulomb and photon induced correlations 28. Whether these modifications of the underlying electronic structure are indeed a major player in the changes of chemical and physical properties still needs to be seen. However, to capture such modifications in the first place (and study their influence) clearly needs an ab-initio theory that is able to treat both types of (strong) correlations accurately and is predictive inside as well as outside of a cavity. We have shown here that dressed RDMFT is a viable option to predict and analyze these intricate structural changes.

Conclusion

In this work we presented a RDMFT formalism for coupled matter-photon systems. This formalism is capable to account for the full quantum-mechanical degrees of freedom of the coupled fermion-boson problem. We discussed that extending the standard formulation of electronic RDMFT to systems with coupled fermionic and bosonic degrees of freedom is not straightforward. Then, we presented an alternative approach which overcomes most of the intricate representability issues by embedding the coupled matter-photon system in a higher-dimensional auxiliary space. Specifically, we introduced for a problem with $N$ electrons coupled to $M$ photon modes, $(N-1)M$ auxiliary coordinates, which allowed us to “fermionize” the coupled problem with respect to new polaritonic coordinates. The resulting dressed fermionic particles are governed by a Hamiltonian with only 1-body and 2-body terms and thus, one can apply any standard electronic-structure method. The extension is constructed in such a way that the auxiliary dimensions do not modify the original physical system and the physical observables are easy to recover. Notably, besides the possibility to study modifications of electronic systems due to a cavity mode, dressed RDMFT offers also the possibility to calculate purely photonic observables like the mode occupation or fluctuations of the electric and magnetic field. We used this framework to develop and implement dressed RDMFT in the electronic-structure code Octopus 88, and tested it with the Hartree-Fock and Müller functional. For simple one-dimensional models of atoms and molecules the obtained approximate results were in good agreement with the exact results from the weak to the deep-strong coupling regime. We then used our method to show that the modifications due to strong matter-photon coupling are far from trivial and depend on the detailed electronic structure. For a molecular as well as an atomic system we showed that strong coupling can locally enhance and suppress the Coulomb induced repulsion between electrons. This behaviour does not only depend on the strength of the matter-photon coupling but also on the details of the matter subsystem (e.g. the interatomic distance of the atoms of a molecule.) We showed that our method allows to predict the structures accurately inside and outside of the cavity and furthermore extends the well-established tools of natural orbitals to analyze coupled light-matter systems.

Although the presented method is practical only for a few photon modes, since the number of photon modes determines the dimension of the involved dressed orbitals, it is exactly these cases that are the most relevant in cavity and circuit QED experiments. Since dressed RDMFT is non-perturbative and seems to be accurate over a wide range of couplings, it is a promising tool to investigate long-standing problems of quantum optics, such as the quest for a super-radiant phase in the ground state of strongly-coupled matter-photon systems 24. Moreover, it is a very promising tool to investigate changes in the ground-state due to matter-photon coupling that can possibly modify chemical reactions 5. Recently, it has been shown that charge-transfer processes can be considerably modified due to strong coupling to a cavity mode28. And although the presented results were for reduced dimensionality, an extension to three spatial dimensions is straightforward. We can rely here again on an already existing implementation for RDMFT in Octopus. Work along these lines is in progress. Besides such interesting applications and fundamental questions of light matter-interactions, there are many open questions to answer also in the presented theory itself. For instance, how strong is the influence of the hitherto negelected $q$-exchange symmetry? First calculations for many particles indicate that it will become important to enforce this extra symmetry to stay accurate when going from the weak to the deep-strong coupling regime. Furthermore, it might become beneficial to avoid the “fermionization” that we employed and then very interesting mathematical questions about $N$-representability for coupled fermion-boson systems need to be addressed. Here, the understanding how to enforce the $q$-exchange symmetries in the dressed formulation could be very useful.

{acknowledgement}

F.B. would like to thank Nicole Helbig, Klaas Giesbertz, Micael Oliveira, and Christian Schäfer for stimulating and useful discussions. We acknowledge financial support from the European Research Council (ERC-2015-AdG-694097).

{suppinfo}

Survey on the bosonic symmetry of the photon wave function. Details about the convergence study of the numerical examples shown in the paper. Protocol for the convergence of a dressed HF/RDMFT calculation.

{tocentry}