Reduced Density-Matrix Approach to Strong Matter-Photon Interaction.

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 nonperturbative 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.

Experiments performed in the last decades (see, e.g., refs (1 and 2)) have made accessible the strong and ultrastrong (we 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 nonlinear 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−16] and much exciting new physics was predicted[17−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.[3,22] However, limits and difficulties of few-level approximations have been pointed out,[23−28] and recently, new models have been used to investigate polaritonic chemistry.[15,28−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 molecules.[17,42] Another example is the ongoing discussion on the theoretical understanding of super-radiance.[24,43] These debates suggest that there is a need for new theoretical tools that treat matter and photons at the same level of theory.[27,28,35] Up to now, standard theoretical modeling 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. (Imagine, for instance, a four-electron system in one spatial dimension, coupled to one cavity mode (which corresponds to the beryllium example, presented in Numerical Results). Every eigenstate of the corresponding Hamiltonian would be a function that depends on five 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 of 0.2) and 8 byte per function value (double precision), we would need 1005 × 8 Byte ≈ 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−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 systems.[49] Recently, quantum-electrodynamical density-functional theory (QEDFT) was introduced as an extension to pure electronic density-functional theory (DFT).[50−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 noninteracting 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 ultrastrong[59] or deep-strong coupling regime[60] of light and matter from first-principles, that is, 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] Similar 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 exploring 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 Figure ). It turns out that transferring RDMFT to such systems involves overcoming additional difficulties in contrast to the DFT framework. The reason behind this is 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−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 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.

Figure 1

Typical setting of a cavity experiment. A matter system (here represented by a diatomic molecule) is put inside an optical cavity that enhances specific modes of the electromagnetic field (here represented by the lowest cavity mode, but in principle many modes can become important). By that, the coupling between the matter system and the light modes can be considerably enhanced with respect to the free space. The dipole of the molecule should be aligned with the polarization of the enhanced mode and its position is assumed at the field maximum. Note that, in principle, also higher multipole moments can become important.

This article is structured as follows. First, we explain the physical setting of an electronic system in a cavity in Physical Setting. In 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 “Fermionization” of Matter–Photon Systems, we introduce the auxiliary system that will allow us to avoid most of the aforementioned difficulties. We show how to construct a proper RDMFT framework in this system in Dressed Reduced Density-Matrix Functional Theory, explain our numerical implementation in Numerical Implementation, and present the numerical results in Numerical Results. We finish with discussing possibilities as well as challenges of dressed RDMFT in 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 (the extension to also include the nuclei as quantum particles is, in principle, straightforward by following, e.g., similar strategies like discussed in refs (27 and 55)). Second, we work in the nonrelativistic limit, which for the typical energy scales of molecules and their low-energy excitations is usually sufficient. Third, we assume that the wavelength λ of the relevant electromagnetic modes is much larger than the spatial extension d of the electronic system (λ ≫ 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.[27,57] 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)Here the first two sums on the right-hand side correspond to the usual electronic many-body Hamiltonian Ĥ = T̂ + V̂ + Ŵ, used to describe the uncoupled matter system consisting of N electrons in an external potential v(r) interacting via the Coulomb repulsion w(r, r′). The third sum describes M photon modes that are characterized by their elongation pα, frequency ωα, and polarization vectors λα. The polarization vectors include already the effective coupling strength (35) and couple to the total dipole D̂ = ∑r of the electronic system. The sum can be decomposed in a purely photonic part , the dipole self-interaction (note that this term is necessary for the existence of a ground state[70]) Ĥ = ∑α=11/2(λα · D̂)2, which we split for later convenience in its one-body Ĥ(1) = ∑α=1∑1/2(λα · r)2 and two-body part Ĥ(2) = ∑α=1∑1/2(λα · r)(λα · r), and the bilinear interaction Ĥ = −∑α=1ωαpαλα · 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 field.[71] However, if we assume a cavity then the Coulomb interaction is modified.[72] We can easily incorporate this into our framework since w(r, r′) is completely at our disposal and we do not rely on any kind of Coulomb approximation. Thus, we can also treat the influence of, for example, a plasmonic environment.[73]

The ground-state wave function of eq is a function of 4N + M coordinateswhere σ is the electronic spin degree of freedom. The wave function Ψ is, as usual, antisymmetric with respect to the exchange of any two electron coordinates rσ ↔ rσ, and also depends on M photon-mode displacement coordinates pα.

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α with pβ. 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 and to represent the sum of harmonic Hamiltonians, that is,Eigenstates of the individual âα+âα in the above representation are given by multiple applications of creation operators to the vacuum state |0⟩, that is, |φα⟩ = (âα+)|0⟩. (Note that, as in refs (47 and 74), we chose here explicitly a non-normalized basis {|φα⟩} of the n-photon sector, with ⟨φα|φα⟩ = n!, which will be convenient for the discussion of bosonic reduced density matrices in the next section. The missing normalization factor is shifted to the resolution of identity, that is, = 1/Nb!|α1,...,α⟩⟨α1,...,α|, where |α1, ..., α⟩ = âα+···âα+|0⟩, as defined later in the text.) These eigenstates are connected to the displacement representation by φα(pα) = ⟨pα|φα⟩ = ⟨pα|(âα+)|0⟩. We can then express an M-mode eigenfunction as ϕ(p1, ..., p) = ⟨p1···p|ϕ⟩ = ⟨p1···p| (â1+)···(â+)|0⟩. In this form it becomes clear that a multimode eigenstate |ϕ⟩ can be considered to consist of N = n1 + ... + n photons (mode excitations). We can associate every such multimode eigenstate with a specific photon-number sector, that is, the zero-photon sector is merely one-dimensional and corresponds to |ϕ0⟩ ≡ |0⟩, the single-photon sector is M-dimensional and corresponds to the span of âα+|0⟩ ≡ |α⟩ for all α and so on. For the multiphoton sectors we see due to the commutation relations of the ladder operators the bosonic exchange symmetry appearing, for example, âα+âα+|0⟩ = âα+âα+|0⟩ ≡ |α1, α2⟩ for α1, α2 ∈ {1, ..., M}. It is no accident that the bosonic symmetry becomes explicit in this representation since the different modes α determine how the photon wave functions look in real space (see also the discussion in Supporting Information, section 5). A general photon state can therefore be represented by a sum over all photon-number sectors as where Φ̃(α1, ..., α) = ⟨α1, ..., α|Φ⟩.

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 valueThis is the classical formulation of the variational principle due to Ritz and is well-defined for every Hamiltonian that is bound from below. While well-known, eq 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 ψ(r1σ1, ..., rσ) of N electrons (see Table ), the expectation value of a general (nonlocal) q-body operator Ô(r1, ... r; r1′, ... r′), which is given by O = ⟨ψ|Ôψ⟩, can be determined via the electronic (spin-summed; we define here only the spin-summed version of the q-body RDM (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) qRDM (note that there are different conventions in the literature for the normalization of the qRDM; we followed ref (75))For instance, the well-known electronic 1RDM, that we denote in the following by γ(r, r′) = Γ(1)(r;r′), 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 Γ(2)(r1,r2;r1,r2). In the chosen normalization, all RDMs satisfy the sum-rule . So for instance, the 1RDM integrates to the particle number ∫d3rγe(r, r) = N, the 2RDM integrates to two times the number of pairs, and so on. Additionally, higher and lower order RDMs are connected via

Table 1

Physical Wave Functions and the Corresponding Symbols of This Section

Ψ(r1σ1, ..., rNσN;p1, ..., pM)	electron-photon many-body state
ψe(r1σ1, ..., rNσN)	purely electronic many-body state
ψb(α1, ..., αNb)	photonic many-body state in mode representation with fixed particle number
Φ(α1,α2, ...)	photonic many-body state in Fock space
ϕei (r)/ϕbi(α)	electronic/photonic natural orbital

For bosons, the same construction is possible. In our case, where we have a discrete set of possible single-boson states, a N boson state in the (symmetrized) mode-representation ψ(α1, ..., α) (see Supporting Information, section 5, for more details) leads to the corresponding bosonic qRDM,According to the electronic case, we denote the 1RDM by γ(α, β) = Γ(1)(α; β). However, in the specific case of photons, where the number of particles is undetermined and we work with Fock-space wave functions |Φ⟩, we need to consider a Fock-space 1RDM of the formIn an according manner one can define a bosonic Fock-space qRDM via Γ((α1, ..., α; α1′, ..., α′) = ⟨Φ|âα+···âα+âα···âα Φ⟩.

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 , we can accordingly define Γ( ≡ d, as well as Γ( ≡ ⟨Ψ|âα+···âα+âα···âαΨ⟩.

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 . 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 Ĥe and the single-particle part of the dipole self-energy Ĥd(1) are given in terms of the 1RDM byHere we have denoted on the left-hand side the explicit dependence of the expectation value on the 1RDM, and the subscript | indicates that r′ is set to r after the application of the semilocal single-particle operator (−1/2∇2 + v(r)). The expectation value of the electronic interaction energy Ŵ and the two-body part of the dipole self-energy Ĥ(2) are given in terms of the (diagonal) of the 2RDM byHence, 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 findHowever, the bilinear coupling term is not given in a simple RDM form but becomesA new reduced quantity appears that mixes light and matter degrees of freedom and can be interpreted as a 3/2-body operator Γ(3/2)(α; r, r′). (To see this in a simple manner, we also lift the continuous fermionic problem into its own Fock space and introduce genuine field operators ψ̂†(rσ) and ψ̂(rσ) with the usual anticommutation relations. Similar 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 ⟨Ψ|[(âα+ + âα)λα · D̂ ]Ψ⟩ = ∑σ∫d3r⟨Ψ|[(âα+ + âα)ψ̂†(rσ)ψ̂(rσ)(λα · r)]Ψ⟩. If we then define Γ(3/2)(α; r, r′) = ∑σ ⟨Ψ|[(âα+ + âα)ψ̂†(rσ)ψ̂(r′σ)]Ψ⟩, we can rewrite ⟨Ψ|[(âα+ + âα)λα · D̂]Ψ⟩ = ∫d3r(λα · r)Γ(3/2)(α; r, 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 qRDM, even if we extend the definitions to include combined matter–boson qRDMs. (Using the field-operator formulation, the usual qRDMs consist of strings of particle-number-conserving combinations of electron and boson operators. Integrating/summing out a number-nonconserving part of it does not lead to a simple relation to half-body RDMs, in general.) One obvious reason is that qRDMs conserve particle numbers, while half-body RDMs do not. Take, for instance, the Fock-space γ(α, β) = ⟨Φ|âβ+âαΦ⟩. In the special case that |Φ⟩ 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 γ(α, β) = dβ*dα, where dα is the total displacement of the coherent state of mode α. In this case, we also know ⟨Φ|âβ+Φ⟩ = dβ*. If we now assume all but one mode, say mode 1, having zero displacement, then we only know γ(1, 1) = |d1|2 from the bosonic 1RDM. We do, however, in general, not know what d1* is. For other states, such a connection is even less explicit.

Putting the inter-relations among the different RDMs aside for the moment, the minimization for the coupled matter–boson problem can be reformulated bySo, in principle, we could replace the variation over all wave functions, Ψ, by their respective set of RDMs needed to define the energy expectation values. Instead of varying over the full configuration space (r1σ1, ..., p), the above reformulation seems to indicate that we can replace this by varying over (r, r′) for the diagonal of Γ(2) and also for the 1RDM γ, together with a variation over (α, β) for γ and over (α,r) for Γ(3/2). 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 , where {γ, Γ(2), γ, Γ(3/2)} → Ψ 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−67] However, only for the 1RDM of ensembles (fermionic or bosonic) are the conditions simple. In this case, by diagonalizing the 1RDM in its eigenbasis γ = ∑n(ϕ)*ϕ, where the ϕ are called the natural orbitals and the n are the natural occupation numbers, the conditions arefor fermions and bosons, respectively. If the particle number N of one species of the system is conserved, the respective sum-rulebecomes 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 qRDMs 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 “fermionize” the coupled fermion-boson problem in such a way that can apply the known conditions of the fermionic problem.

“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 modeled with new particles that are fermions. We call them dressed or polaritonic particles, because they depend on electronic and photonic coordinates. (Note 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 into three steps. First, we introduce for each mode auxiliary extra dimensions (pα,2, ..., pα,), 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, that is, we now consider wave functions depending on (r1 σ1, ..., r σ, p1, ..., p, p1,2, ..., p1,, ..., p, ..., p). Second, we add for every photon mode α a linear operatorto the physical Hamiltonian of eq . 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 Ĥ′ = Ĥ + ∑α=1Π̂α (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α,1, ..., qα,) such thatThis 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 aswhere we inserted the definition of the dipole operator, D̂ = ∑r, 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 z = (r, q1, ..., q), we can rewrite the above Hamiltonian aswhere we introduced the dressed Laplacian Δ′ = ∑3 + ∑α=1, the dressed local potential v′(z) = v(r) + and the dressed interaction kernel w′(z,z′) = w(r,r′) + ∑α=1q · r′ − q′ · r + · rλ · r′]. Note that the choice of the linear auxiliary operator (eq ) and the coordinate transformation (eq ) have a certain freedom. The operator of eq 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 . These requirements are met using an orthonormal transformation together with the harmonic oscillators of eq . (Note 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 . One specific example of such a transformation is[68]p = (q1 + ... + q − (k − 1) q) for 2 ≤ k ≤ N, alongside the first line of eq .) Since the operator of eq only acts on the auxiliary coordinates, the normalized physical solution Ψ of the original (time-independent) Schrödinger equation E0Ψ = ĤΨ in the standard configuration space is connected to a new physical solution of the auxiliary Hamiltonian Ĥ′ in a very simple manner, that is,Here the normalized solution Ψ′ of the auxiliary Schrödinger equation E0′Ψ′ = Ĥ′Ψ′ is found with χ being the ground state of ∑αΠ̂α. The ground state χ is merely a tensor product of individual harmonic-oscillator ground states, and therefore, exchanging pα, with pα, does not change the total wave function Ψ′. If we rewrite this wave function in the new coordinatesand due to the fact that we constructed χ to be symmetric with respect to the exchange of qα, and qα,, we realize that Ψ′ is antisymmetric with respect to the exchange of (zσ) and (zσ). (Note that the auxiliary ground state of mode α is given as a Gaussian with respect to pα,2 = qα,2 – 1/N(qα,)2, where we used p = 1/√Nqα,.) Thus Ψ′ is fermionic with respect to the polaritonic coordinates (zσ) and can be represented by a sum of Slater determinants of (3 + M)-dimensional polaritonic orbitals φ(zσ). This makes the application of the 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α, ↔ qα, 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 Ψ′ is a minimal-energy state in the extended space with the appropriate symmetries, we first point out that for any trial wave functionthe q-exchange symmetry implies that the fermionic symmetry is in the (rσ) coordinates. Then it holds since Ĥ only acts on (r1, ..., p), and ∑α=1Π̂α only acts on p1,2, ..., p so thatAlthough we constructed the auxiliary space explicitly in a way that the physical wave function Ψ can be reconstructed exactly from its dressed counterpart Ψ′, by integration of all auxiliary coordinates, this does not hold for all types of operators. (Note that Ĥ′ also has many eigenstates that are not of the form Ψ′ = Ψχ, with Ψ antisymmetric under exchange of rσ ↔ rσ and χ symmetric under exchange of qα, ↔ qα,. 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α, ↔ qα,, together with the antisymmetry under exchange of zσ ↔ zσ.) For operators that depend only on electronic coordinates, there is no difference, and we have ⟨Ψ|ÔΨ⟩ = ⟨Ψ′|ÔΨ′⟩. This is not surprising because the coordinate transformation (eq ) 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, that is, any operator that depends only on the elongation of (qα,1 + ... + qα,) and its conjugate , 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 , that we can straightforwardly generalize in the auxiliary space to . The connection between Ĥph and Ĥph′ is given by the definition of the coordinate transformation (eq ),Since the expectation value of Π̂α is known analytically, , we have . This can be generalized to any operator that contains terms of the form and , where α and β denote any two modes. The reason is that the transformation (eq ) 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 α, . In the auxiliary system, we have

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 Ĥ′, given in eq . It consists of only polaritonic one-body terms ĥ(1)(z) = −1/2Δ′ + v′(z), and two-body terms ĥ(2)(z, z′) = w′(z, z′). It commutes with the polaritonic particle-number operator N̂′ = ∫d3+zn̂(z), where we used the definition of the polaritonic local density operator n̂(z) = ∑δ3+(z – z). This means that the auxiliary system has a constant polaritonic particle number N. Additionally, the physical wave function of the dressed system Ψ′(z1 σ1, ..., z σ) is per construction antisymmetric. This allows for the definition of a dressed (spin-summed) 1RDMin accordance to eq . By construction, this auxiliary density matrix reduces to the physical electron density matrix via γ(r, r′) = ∫dqγ(r,q1, ..., q;r′,q1, ..., q). Furthermore, we introduce the (spin-summed) dressed 2RDM Γ(2)(z1,z2;z1′,z2′) = N(N – 1)∑σ∫d(3+zΨ′*(z1′σ1, z2′σ2, z3σ3, ..., zσ)Ψ′(z1σ1, z2σ2, z3σ3, ..., zσ). These dressed RDMs allow for expressing the energy expectation value of the dressed system byThus, we can define the variational principle for the ground state only with respect to well-defined reduced quantities,To perform this minimization, we need to constrain the configuration space to the physical dressed RDMs that connect to an antisymmetric 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, for example, ref (82). Instead, we want to stick to the dressed 1RDM γ and approximate the 2-body part as a functional of γ. 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′ of only the system’s dressed 1RDMFor 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 γ. 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 . 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., eqs and 9). The price we pay for this is 2-fold: 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 QED.[30,40,41,52,55,68]

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 last paragraph of “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 conditions (this 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 (in 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′[γ] that have been developed for the electronic case. To do so, we further, similarly to the electronic case, decompose W′[γ] = E[γ] + E[γ] into a classical Hartree part E[γ] = 1/2∫∫d3+zd3+z′γ(z, z)γ(z′,z′)w′(z, z′) and an unknown exchange-correlation part E[γ]. Almost all known functionals E[γ] are expressed in terms of the eigenbasis and eigenvalues of the 1RDM. In our case, the dressed natural orbitals ϕ(z) and occupation numbers n are found by solving ∫d3+z′γ(z,z′)ϕ(z′) = nϕ(z). The simplest approximation is to only retain the fermionic exchange symmetry and employ the Hartree–Fock (HF) functionalAs 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 nonlinear occupation-number dependence in the exchange-correlation functional. We here consider the Müller functional[62]which has been rederived 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, thus, we only want to remark that a variety of functionals were proposed after E[γ] 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 reuse many ideas from electronic RDMFT, a further advantage of the dressed reformulation is that we can also reuse 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 counterpart v′. This, together with a change of the electronic Coulomb interaction w to its dressed counterpart w′ 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 asWe 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, that is, 0 ≤ n = 2 sin2(θ) ≤ 2, to satisfy eq . (Note that the n is bounded by 2 because we employed a spin-restricted formulation. If we considered natural spin–orbitals instead, the upper bound would be 1). The second part of the conditions (eq ), that is, ∑n = N, as well as the orthonormality of the dressed natural orbitals, that is, ∫d3+zϕ*(z)ϕ(z) = δ, are imposed via Lagrange multipliers as, for example, explained in ref (88). We have available two different orbital-optimization methods, a conjugate-gradient algorithm (we used this method only for some benchmark calculations, 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 ϕ 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 the Supporting Information, section 4.

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 nontrivial 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 nontrivial dependence on the interatomic 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 . Specifically, we consider a one-dimensional model of a helium atom (He), that is, , a one-dimensional model of a hydrogen molecule (H2), that is, vH(x) = − , first at its equilibrium position d = deq = 1.628 au, later with varying d, and a one-dimensional model of a Beryllium atom Be, that is, . For the latter, we consider a smaller softening parameter ϵ = 0.5 to make sure that all electrons are properly bound. We use the soft Coulomb interaction (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 ωHe = 0.5535 au and ωH2 = 0.4194 au. For Be instead, we choose ωBe = 3.0 au, 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 ω 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 ω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 ω. 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 the Supporting Information, section 2.3.1.

We start with the discussion of the (N = 2)-particles examples. In this setting, the dressed auxiliary system is four-dimensional (two particles with two coordinates each), which is still small enough to be solved exactly in a four-dimensional discretized simulation box, so that we can compare dressed RDMFT (with the Müller functional of eq ), dressed HF (see eq ) and the exact solutions. We used the box lengths of L = L = 16 au and spacings of dx = dq = 0.14 au to model the electronic x and photonic coordinates q of the two dressed particles in the exact routine. (We 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 the Supporting Information, section 1). For dressed RDMFT and dressed HF instead, we considered two-dimensional simulation boxes for every dressed orbital and we set L = L = 20 au and dx = dq = 0.1 au. We obtained converged results for () natural orbitals for He (H2). Details about how we determined these parameters can be found in the Supporting Information.

We first show (see Figure ) 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/ω (this quantity is typically used as a measure for the strength of the light–matter interaction, see, e.g., ref (28)) for He and H2, respectively. We thereby go from weak to deep-strong coupling with g/ω = 1.[3] The deep-strong coupling regime has been reached in different systems like for instance for Landau polaritons.[60] For (organic) molecules the highest reported coupling strengths are in the ultrastrong regime of g/ω ≈ 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 ρ(x, q) ≡ γ(x, q; x, q). To simplify this discussion, we separate the electronic and photonic parts of the two-dimensional density by integration, i.e., ρ(x) = ∫dqρ(x, q) and ρ(q) = ∫dxρ(x, q). The exact reference solutions show that with increasing g/ω, the electronic part of the density becomes more localized, while the photonic part becomes broadened (the reader is referred to the last two paragraphs of this section and Figures and 11 for details about the effects of matter−photon coupling). 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/ω = 0.1 in Figure 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 overlocalize ρ(x) (here visible for H2).

Figure 2

Differences of dressed HF (dHF) and dressed RDMFT (dRDMFT) from the exact ground state energies (in Hartree) as a function of the coupling g/ω for the (one-dimensional) He atom (left) and (one-dimensional) H2 molecule (right) in the dressed orbital description. Dressed RDMFT improves considerably upon dressed HF. For both systems, the energy of dressed RDMFT remains close to the exact one, the error of dressed HF instead increases with the coupling strength.

Figure 10

We show the differences in the electronic density of the H2 molecule for three different bond lengths d (as examples of the dissociation) for g/ω = 1.0 compared to g/ω = 0, calculated exactly (ρex(x), left) and with dressed RDMFT (ρdR(x), right). We see that for small d, the cavity mode reduces the electronic repulsion and localizes the charges at the bond center (d = 1 < deq = 1.628) in comparison to the free molecule (insets). For larger d, the electronic repulsion is locally enhanced such that the charge deviations are separated in two peaks (d = 2). For very large d, this interplay between local suppresion and enhancement of repulsion becomes more pronounced (d = 3). The dressed RDMFT calculations capture the behavior very well.

Figure 11

We show the differences in the electronic density (ρ(x)) of He (left) and Be (right) for three different coupling strengths compared to the atoms outside the cavity (insets), calculated with dressed RDMFT. We see that the effect of the cavity is very different for both systems: The strong localization of the electronic density for He indicates the suppression of electronic repulsion for all coupling strengths. For Be instead, we see additionally local enhancement of the repulsion. The interplay of enhancement and suppression changes with increasing coupling strength.

Figure 3

Deviations of dressed HF (dHF) and dressed RDMFT (dR) ground state densities from the exact solution (ρex, depicted in the insets) for the He atom (top) and the H2 molecule (bottom) with coupling g/ω = 0.1. We separate the electronic (x, left) and photonic (q, right) coordinates as explained in the text. For both systems, dressed RDMFT finds a considerably better electronic density than dressed HF, which is consistent with the better result in energy (see Figure ). The photonic densities are reproduced almost exactly for both levels of theory.

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 ϕ(x, q). Again, we separate electronic and photonic contributions and show their reduced electronic density ρ(x) = ∫dq |ϕ(x, q)|2. Figure 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 density is only slightly different, the third natural orbital density of H2 differs 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 (H2: n1 = 1.878, n2 = 0.102, n3 = 0.015; He: n1 = 1.978, n2 = 0.020, n3 = 0.001), their (inaccurate) contribution to the density and total energy remains small.

Figure 4

First three natural orbital densities ρex/dR((x) of the exact (ex) and dressed RDMFT (dR) calculations for the He atom (top) and the H2 molecule (bottom) with coupling g/ω = 0.1. We see in both cases that ρex(1)(x) is almost exactly reproduced by dressed RDMFT, but ρdR(2)(x) deviates already visibly from ρex(2)(x) (left). However, it is in both cases qualitatively correct. This changes for ρdR(3)(x) of H2, which has one node more than ρex(3)(x). For He instead, ρdR(3)(x) is reproduced correctly (right).

To complete the picture, we also look at the photonic natural orbital densities, ρ(q) = ∫dx|ϕ(x, q)|2, the first three of which are plotted in Figure , for He and H2. 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 ultrastrong coupling regime. The accuracy drops with increasing g/ω.

Figure 5

We show the differences Δρ( = ρdR((q) – ρex((q) between the dressed RDMFT (dR) and the exact (ex) photonic natural orbital densities ρex/dR(q) for the three highest occupied natural orbitals for the He atom (left) and the H2 molecule (right) for coupling strength g/ω = 0.1. For both systems, the exact ρex((q) have a similar shape as the density (see inset). We see in both cases that dressed RDMFT captures the exact solution very well.

As an example for a photonic observable, we show in Figure the mode occupation Nph as a function of the coupling strength g/ω that we calculated by using eq , that is, Nph = − , with the photon mode energy From weak to the beginning of the ultrastrong coupling regime (g/ω ≈ 0.1), both dressed HF and dressed RDMFT capture Nph well. For very large coupling strengths, the deviations to the exact mode occupation becomes sizable. This might sound counterintuitive, 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 Figure ). However, it is probable that this can be improved by better functionals.

Figure 6

Total mode occupation Nph, calculated from the exact, dressed HF and dressed RDMFT solutions for He (left) and H2 (right). We see that both dressed RDMFT and dressed HF underestimate Nph. In the ultrastrong coupling regime for g/ω > 0.3 both dressed HF and dressed RDMFT (with the Müller functional) deviate strongly from the exact solution.

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 system that cannot easily be solved exactly: the one-dimensional Be-atom in a cavity. In Figure , we see the total energy as a function of the coupling strength g/ω for dressed HF and dressed RDMFT, respectively. Like in the two-electron systems, the deviation between both curves increases for larger g/ω 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 two-particle systems. With increasing g/ω, 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 Figure 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 Figure shows. We conclude the survey of Be with the mode occupation under variation of the coupling strength (see Figure ). We see that the value of g/ω ≈ 0.5 separates two regions. For g/ω < 0.5, dressed RDMFT finds a larger mode occupation than dressed HF, and for g/ω > 0.5 instead, the dressed HF mode is more strongly occupied. We found similar behavior also for the two-particle systems, although the boundary between the two regions was considerably different there (He: g/ω ≈ 0.8, H2: g/ω ≈ 0.1, see Figure ).

Figure 7

Total energy of the dressed HF and dressed RDMFT calculations of Be for increasing g/ω. We observe the same trend as for the two-electron systems: for both levels of theory, the energy grows with increasing g/ω, though for dressed HF faster than for dressed RDMFT.

Figure 8

Shown are the electronic (ρdHF/dR(x), left) and photonic (ρdHF/dR(q), right) densities of Be for dressed HF (dHF) and dressed RDMFT (dR) for two different coupling strengths subtracted from their counterparts in the no-coupling limit (ρdHF/dR(x/q)). We see in the electronic (photonic) case that the dressed RDMFT deviations are less (more) pronounced than for dressed HF.

Figure 9

Total mode occupation Nph of Be for dressed HF and dressed RDMFT. We see that dressed RDMFT exhibits larger Nph until a coupling strength of g/ω ≈ 0.5. For larger coupling the dressed HF mode occupation becomes higher.

We conclude this section with two examples, for which the light-matter interaction changes the bare systems nontrivially, depending not only on the coupling strength but also on the details of the electronic structure. We start with the dissociation of H2 as an example for a chemical reaction, where we use vH(x) with different d. In Figure , 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 toward 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/ω = 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/ω = 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 Figure ), they behave very differently under the influence of the cavity, which can be seen in Figure . The electronic density of He is pushed toward 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. (For g/ω = 0.8, we still observe n1 = 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 n2 = 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/ω = 0.1 and 0.4 in the plot) but also similar effects (g/ω = 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 nontrivial 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. It has been shown that such small changes are capable to strongly affect chemical properties and reactions, 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 an 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 one-body and two-body terms and, thus, can be applied to 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 behavior 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 nonperturbative 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 mode.[28] 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 of 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.