Light-matter hybrid-orbital-based first-principles methods: the influence of the polariton statistics

A detailed understanding of strong matter-photon interactions requires first-principle methods that can solve the fundamental Pauli-Fierz Hamiltonian of non-relativistic quantum electrodynamics efficiently. A possible way to extend well-established electronic-structure methods to this situation is to embed the Pauli-Fierz Hamiltonian in a higher-dimensional light-matter hybrid auxiliary configuration space. In this work we show the importance of the resulting hybrid Fermi-Bose statistics of the polaritons, which are the new fundamental particles of the ``photon-dressed'' Pauli-Fierz Hamiltonian for systems in cavities. We show that violations of these statistics can lead to unphysical results. We present an efficient way to ensure the proper symmetry of the underlying wave functions by enforcing representability conditions on the dressed one-body reduced density matrix. We further present a general prescription how to extend a given first-principles approach to polaritons and as an example introduce polaritonic Hartree-Fock theory. While being a single-reference method in polariton space, polaritonic Hartree-Fock is a multi-reference method in the electronic space, i.e. it describes electronic correlations. We also discuss possible applications to polaritonic QEDFT. We apply this theory to a lattice model and find that the more delocalized the bound-state wave function of the particles is, the stronger it reacts to photons. The main reason is that within a small energy range many states with different electronic configurations are available as opposed to a strongly bound (and hence energetically separated) ground-state wave function. This indicates that under certain conditions coupling to the quantum vacuum of a cavity can indeed modify ground state properties.

results. We present an efficient way to ensure the proper symmetry of the underlying wave functions by enforcing representability conditions on the dressed one-body reduced density matrix. We further present a general prescription how to extend a given first-principles approach to polaritons and as an example introduce polaritonic Hartree-Fock theory. While being a single-reference method in polariton space, polaritonic Hartree-Fock is a multi-reference method in the electronic space, i.e. it describes electronic correlations. We also discuss possible applications to polaritonic QEDFT. We apply this theory to a lattice model and find that the more delocalized the bound-state wave function of the particles is, the stronger it reacts to photons. The main reason is that within a small energy range many states with different electronic configurations are available as opposed to a strongly bound (and hence energetically separated) ground-state wave function. This indicates that under certain conditions coupling to the quantum vacuum of a cavity can indeed modify ground state properties.
A plethora of experiments of atoms, molecules and solids embedded in quantum cavities 1-9 that were performed in the last two decades, have demonstrated the possibility to change the properties of matter by coupling it strongly to the modes of an optical cavity.
In the strong-coupling regime, 1 matter degrees of freedom strongly mix with a few effective photon modes such that hybrid light-matter states, called polaritons, emerge. The combined light-matter system can exhibit significantly different properties than the separate subsystems even at ambient conditions, which suggests various interesting applications in chemistry and material science. [1][2][3] Examples include the possibility of building polariton lasers, 5 the modification of chemical landscapes, 6 selective manipulation of electronic excitations, 7 the control of long-range energy transfer between different matter systems 8 or the emergence of new states of matter. [13][14][15] The variety of this (by far not exhaustive) list of effects caused by the emergence of polaritons reveals the complexity that arises when light and matter mix strongly. The theoretical description of these effects is far from trivial. Currently, many fundamental 1 For a definition and detailed discussions of light-matter coupling regimes, see Refs. [10][11][12] questions in the field of polaritonic chemistry are yet to be answered. Among others 16,17 are the questions whether collective strong-coupling can modify chemistry without influencing the local electronic structure 18,19 and whether strong coupling to only the vacuum field can modify chemical reactions. 20 Letting aside the cavity setting for a moment, in modern chemistry many fundamental questions about the behavior of matter have been answered by first-principle calculations.
These well-established methods routinely predict properties of matter, which might interact with "classical" time and spatially dependent electromagnetic fields. For example, in order to answer whether a specific chemical reaction happens or not, one calculates the potential-energy surfaces, which allows for the estimation of activation barriers. However, such complex first-principle methods (like Density Functional Theory (DFT) 21 or computational quantum-chemistry methods 22 ) are usually geared towards the precise quantummechanical description of the electrons (treating ions classically). Nevertheless, many phenomena require a quantum mechanical description of the ion dynamics. In a similar way, materials in cavities call for a quantum treatment of the electromagnetic field. [23][24][25][26] In the cavity setting (see Fig. 1 for a sketch), where besides the electrons and nuclei also the photons play a decisive role, first-principle calculations for only one species of particles are already routinely employed. Although only to provide parameters (such as dipole moments and excitation energies) which serve as essential input in cavity Quantum Electrodynamics (QED) models such as the Rabi and the Jaynes-Cummings model. 2 Photon properties and absorption spectra can be described well within these models in combination with, e.g.
input-output theories. 27,28 However, when it comes to polaritonic chemistry, first-principle theories for light and matter can provide different (complementary) information, 13 such as modified Cavity Born-Oppenheimer surfaces 23 or the calculation of the most important states in a chemical reaction. Thus, the profit of extending first-principle methods to the cavity setting is twofold: on the one hand they can describe coupled systems parameter-free, Figure 1: Sketch of our cavity setup. In our explicit example (see Sec. 4) the electrons of our matter subsystem are allowed to move parallel to the electric field E of the cavity mode yet restricted in the perpendicular directions, i.e. we consider a one-dimensional discretrized matter subsystem. Since the extension in the perpendicular directions is small compared to the wave length of the dominant cavity mode, the coupling is mediated via the total dipole d of the electrons. If the mode volume (distance between the mirrors) is small or the number of particles increased, new hybrid light-matter quasi-particles, i.e. polaritons, emerge.
i.e. less biased, and on the other hand they can provide parameters for cavity QED models 29 or even motivate new purpose-build models. This dual role of first-principle theories is well established, for example, in the context of solid-state physics. In this case, first-principle methods are used either for a full description of the solid or to provide parameters for, e.g.
the Hubbard models for more complex (strongly-correlated) solids.
The standard electronic-structure methods are geared towards a specific type of particle, i.e., electron or nucleus/ion. And depending on the situation it is usually one of these two particle species that dominates a physical or chemical property. In the case of strong lightmatter coupling, where hybrid light-matter states can change these properties, also at least a second particle type, the photon, becomes important as well. Thus one needs to develop specific first-principle methods that can treat several quantized particle species at the same time, such as the Kohn-Sham approach to Quantum-Electrodynamical Density Functional Theory (QEDFT). 25,[30][31][32] However, it has been recently shown that first-principles methods, such as QEDFT, can be based on the emerging hybrid particle, the polariton, directly.
For that the usual dipole-coupled Hamiltonian, i.e. the Pauli-Fierz Hamiltonian in the long-wavelength approximation, 33 is embedded into a higher-dimensional polaritonic Hilbert space in an exact way. 34,35 The resulting polaritonic or "dressed" Hamiltonian has the same structure as the usual matter-only quantum-mechanical Hamiltonian, i.e., one-body kinetic and external potential part and a two-body interaction term with respect to the (now higher dimensional) polaritonic coordinates. This allows one to employ well-established many-body methods, using polaritonic orbitals. This has been demonstrated for time-independent and time-dependent QEDFT 34 as well as for Hartree-Fock (HF) and Reduced-Density-Matrix Theory (RDMFT). 35 Such an approach has several key advantages: One does not need to invent new first-principle methods for the QED setting from scratch, but can use the machinery of the well-established electronic-structure methods. Further, one works directly in terms of the fundamental quasi-particle that determines the properties of the coupled system. And finally, simple wave functions in terms of polaritonic orbitals correspond to correlated (multi-determinant) wave functions in physical space.
The above presented approach is exact, since it maps the electron-photon (and more generally the electron-ion-photon) configuration space in a simple and explicit manner to a higher-dimensional polaritonic configuration space. Yet the resulting polaritonic wave functions have now hybrid Fermi-Bose statistics. This has interesting consequences for the firstprinciple approaches, which are commonly based on Slater determinants of single-particle orbitals. Yet a Slater determinant enforces purely fermionic symmetry. Thus, for the polaritonic extension of first-principle methods we need a construction that enforces this new type of combined fermion-boson symmetry in an easy and efficient manner. While in previous works this has only been approximately enforced and still allowed good results under specific conditions and for model systems, 34,35 it is to be expected that in general an error in the quasi-particle's statistics can lead to unphysical results.
In this work we show that indeed approximating the new statistics of the polariton wave function can result in unphysical predictions. Specifically, ignoring the hybrid Fermi-Bose character leads to violations of the Pauli principle, i.e., several electrons can occupy the same quantum state. By then analyzing in detail the hybrid character of the polaritonic wave functions we find necessary conditions that any physical polaritonic wave function needs to obey. In the case of ground-state wave functions, these conditions are sufficient to guarantee the physicality of the wave function. We here show how these conditions can be enforced on the level of reduced density matrices for common first-principle methods with the help of exact inequality constraints on the polaritonic Lagrangian. This result makes such methods a straightforward and reliable tool for the prediction of changes due to strong light-matter interactions. As an explicit example we consider how HF theory can be generalized to polaritonic problems and how the inequality constraints lead to new terms in the HF equations that enforce the combined fermion-boson character. We furthermore present numerical results for simple model systems that are, however, already challenging for a straightforward exact numerical calculation. We observe that how multi-electron systems react to strong coupling depends on structural details of the uncoupled system. Most importantly, we find that the more spatially extended the electronic wave function is, the stronger the electrons react to the photons, and the more different the coupled light-matter ground state is when compared to the uncoupled one. This goes together with an increase of the electronic correlations captured by polaritonic HF theory, which is multi-determinantal in the electronic subsystem. These results suggest that coupling to only the vacuum of a cavity can have a strong impact, specifically when the uncoupled electronic wave function is spatially delocalized. This also indicates that when a collective ensemble of emitters is treated from first-principles, there could be strong local modifications, since the collective wave function of all these particles would be highly extended.
The article is structured as follows. After briefly introducing the standard Hamiltonian of coupled electron-photon systems in Sec. 1, we discuss in detail the polariton statistics of the dressed system and elucidate how and why unphysical solutions may occur without the correct statistics in Sec. 2. We provide instructions on how to generalize a given electronicstructure theory to polaritons with the correct statistics in Sec. 3. We conclude this section by explicitly applying our instructions to one of the most ubiquitous approaches, namely HF theory. We then present corresponding numerical results in Sec. 4 and discuss physical implications. We conclude the article with a summary and outlook in Sec. 5.
1 Electrons dipole-coupled to cavity modes As indicated in Fig. 1, we consider an atomic or molecular system inside a cavity or, more generally, a nanophotonic environment. Since we assume that the extension of the matter subsystem is small compared to the wavelength of the cavity (perpendicular to the polarization direction of the dominating modes), we can simplify the full minimal-coupling Hamiltonian by adopting the dipole approximation. 32,33,36 If we furthermore assume the nuclei/ions clamped, i.e. we work in the Born-Oppenheimer approximation, the problem reduces to N electrons coupled to the cavity modes. While in principle we could also include the nuclei/ions in our description, since the structure of the Hamiltonian (which is the crucial ingredient for the dressed formulation, see Sec. 2) would not change, 37 we refrain from this more complex situation and concentrate on the electronic structure. In length form and atomic units the Hamiltonian in dipole-approximation for a molecular system in a quantum cavity then reads 31,36 the kinetic and external one-body parts,T [t] andV [v], respectively, and the two-body interaction termŴ [w]. Here the kinetic term is the usual Laplacian t(r) = − 1 2 ∇ 2 r , the external potential v(r) is due to the attractive nuclei/ions and w(r, r ) is the electron-electron repulsion. Usually this is just taken as the free-space Coulomb interaction w(r, r ) = 1/|r − r |, but in a cavity the interaction can be modified. 38 Since we do not rely on the specific form of the electron-electron repulsion we can, for instance, also easily model the changed interactions due to a plasmonic environment. 39 The fourth termĤ ph is the free field-energy of M effective modes of the cavity. The effective modes are characterized by their displacement coordinate p α , frequency ω α and polarization vectors λ α . The latter include already the effective coupling strength g α = |λ α | ωα 2 ∝ 1/ √ V 10,13 that is proportional to the inverse square-root of the cavity mode volume V. In the dipole approximation, the coupling between light and matter is described by the bilinear termĤ I together with the dipole self-energy termĤ d . Here the dipole operator is defined byD = N k=1 r k . We note that the dipole self-energy termĤ d is of utmost importance, especially if equilibrium properties are to be considered (as is the case in this work). That is, because this term is responsible for the stability of matter coupled to photon modes 40 and also guarantees many further fundamental properties of the coupled light-matter system. 41 As a consequence, only by including this (in cavity QED models often discarded) term we can have a well-defined ground state wave function in the basis-set limit.
The corresponding ground-state wave function of Eq. (1) depends then on 4N + M coordinates Ψ(r 1 σ 1 , . . . , r N σ N ; p 1 , . . . , p M ), where σ k are the electronic spin degrees of freedom. The wave function Ψ is anti-symmetric with respect to the exchange of any two electron coordinates r j σ j ↔ r k σ k , and also depends on M photon-mode displacement coordinates p α . The anti-symmetry in rσ enforces the Fermi statistics and thus the Pauli exclusion principle, i.e. two electrons cannot occupy the same quantum state. It is important to note that there is no fundamental exchange symmetry between different displacement coordinates. This is due to the fact that the space of each mode is a single-state bosonic Fock space, i.e. it counts how many excitations are in that mode, and hence the bosonic symmetry is in the associated creation and annihilation operators (see, e.g., Sec. 2 of Buchholz et al. 35 ).
Finally, a comment on the basic physical entities in length gauge is in order. Since the length gauge mixes matter and photon degrees of freedom, the displacement coordinates p α do not correspond to the electric field but rather (as the name already indicates) to the auxiliary displacement field. 36,40,41 As a consequence, in contrast to the unitarily equivalent velocity form,Ĥ ph is not directly proportional to the number of photons but contains matter contributions. Indeed, in length gauge the photon-number operator for mode α is 41 which contains the interaction terms. This highlights that in length gauge we already implicitly work with light-matter quasi-particles. This is, however, not yet sufficient to see that polaritons emerge as the fundamental quasi-particles of the light-matter theory. How this can be made explicit is discussed in the next section.

Dressed electron-photon system
Before we discuss the construction that makes the polariton explicitly the fundamental quasiparticle of the dipole-coupled Hamiltonian of Eq. (1), we first want to make some general comments about polaritons and the multi-reference character of coupled electron-photon systems.
As already pointed out before, the hallmark of strong electron-photon coupling is the emergence of light-matter hybrid states or polaritons. The basic models 10 to describe such hybrid states consider two relevant electronic states, labelled by |g ("ground") and |e ("excited" state), and two photonic states, the vacuum |0 and one-photon state |1 . The resulting polaritons are then denoted as with coefficients α, β that depend on details of these models. From a quantum-chemical perspective we therefore see that polaritons correspond to multi-reference wave functions, even if we described the electronic states as single Slater determinants. If we go beyond these most simple of models the description of the coupled light-matter states require even more terms for an accurate description, especially for strongly-coupled ground states. 42 Systems that require multi-reference states are well known in electronic-structure theory and their accurate description is among the hardest challenges in the field. An example is bond stretching in the hydrogen molecule, where the wave function has a multi-reference character.
This prototypical system is commonly used as a challenging test case for first-principle methods. [43][44][45][46] However, the term multi-reference depends crucially on the basic entities, i.e.
the "single references", that are used to build the "multi-reference" state. Transferred to the problem of polaritonic physics, we face the fundamental problem that separate matter and photon wave functions 3 (as in the simple example above) can be very inefficient in describing polaritonic states of the coupled system. Instead, we want to describe such systems directly by polariton degrees of freedom that depend on both, photonic and electronic variables.
Specifically, we want to define a one-particle (orbital) basis of polaritons, on which we want to base our considerations. The corresponding polaritonic wave function will differ (however in an explicit and trivial manner, see Eq. (12)) from the physical wave function, which we usually construct from separate electronic and photon wave functions. The main technical advantage of this dressed construction, discussed in the following, is that already a singlereference polariton wave function corresponds to a multi-reference coupled electron-photon wave function and hence polaritonic (or dressed) orbitals are potentially more efficient than using the two separate ones. Figure 2: Sketch of the auxiliary construction for an example wave function Ψ(x, y, p) of two one-dimensional electrons (x, y) coupled to one photon mode with displacement coordinate p. The coupling is indicted by the double arrows x ↔ p, y ↔ p. The electronic orbital wave functions are symbolized by the ground and first excited state of a box with zero-boundary conditions and the photon mode by a wiggly line. The corresponding dressed wave function Ψ (x, q 1 , y, q 2 ) has instead two photon-coordinates q 1 , q 2 that are related to the physical coordinate p by p = 1/ √ 2(q 1 + q 2 ). On the wave function level, this connection can be utilized to introduce two polariton orbitals with coordinates (x, q 1 ) and (y, q 2 ), respectively, that are interacting. This new interaction is indicated by a double arrow between the two orbitals.
While we focus in this work on the correct hybrid statistics of the polaritons in the dressed construction, 34,35 we nevertheless want to give a short recapitulation. To turn the coupled electron-photon problem of Sec. 1 into an equivalent and exact dressed problem we will follow three steps (for a simplified sketch of the construction, see Here and in the following, we will denote all quantities in the auxiliary configuration space with a prime.
2. We next construct an auxiliary Hamiltonian in the extended configuration space of the depends only on these new auxiliary coordinates. This construction guarantees that the auxiliary degrees of freedom do not mix with the physical ones, which will ensure a simple and explicit connection between the physical and auxiliary system.
3. Finally we perform an orthogonal coordinate transformation of the physical and auxiliary photon coordinates (p 1 , ..., p M,N ) → (q 1,1 , ..., q M,N ) such that Note that the second line is automatically satisfied for any orthogonal transformation and the first line defines p α as the "center-of-mass" of all the q α,i with uniform relative In total, we then find the auxiliary Hamiltonian in the higher-dimensional configuration space given aŝ where we inserted the definition of the total dipole operator and reordered the expressions, such that the terms with only one index and the terms with two different indices are grouped together. Introducing then a (3+M )-dimensional polaritonic vector of space and transformed photon coordinates z = rq with q ≡ (q 1 , . . . , q M ), we can rewrite the above Hamiltonian aŝ where we introduced the dressed one-body terms and the dressed two-body interaction term We see here that only the conditions (6), but not the details of the coordinate transformation of step 3 are important for our construction. 35 We want to stress that the crucial part of this coordinate transformation is the replacement of p α in the interaction terms p α λ α ·D. Instead of p α only, now all q α,i couple to the dipole of the matter system just with a rescaled couplingstrength by the factor 1/ √ N (see Appendix A for an explicit example of such a coordinate transformation). Further, as will be discussed in more detail in Sec. 3, we are in practice not interested in an exact solution of the dressed problem. Since the auxiliary configuration space is much larger than the original configuration space, we made an (beyond simple systems) unfeasible numerical problem even more unfeasible. Yet, the dressed formulation allows for relatively simple approximation schemes, e.g. HF theory in terms of a single polaritonic Slater determinant (see Sec. 3).
To conclude this brief summary of the dressed construction, we want to mention that while the matter observables stay unchanged, the photon observables like the photon energŷ H ph or the displacement coordinate of the photon modes p α have different representations in the auxiliary configuration space. We can, however, explicitly derive the new representations as discussed in detail in Buchholz et al. 35 The most important examples that we will also use in this work are the photon energy that enters the total energy Let us next discuss the wave function Ψ in the auxiliary configuration space. The wave function Ψ in the usual configuration space is a (normalized) solution of the (timeindependent) Schrödinger equation E 0 Ψ =ĤΨ. SinceĤ =Ĥ + M α=1Π α andΠ α acts only on the auxiliary coordinates, we can simply construct with χ being the (normalized) ground state of M α=1Π α , which is a product of individual harmonic-oscillator ground states. Clearly, Ψ is a normalized solution of the auxiliary Schrödinger equation E 0 Ψ =Ĥ Ψ . In principle any combination of eigenstates of the auxiliary harmonic oscillators would lead to a new eigenfunction forĤ but since we here focus on the ground state the natural choice is the lowest-energy solution. Rewriting this wave function in the new coordinates and employing the polaritonic coordinates zσ ≡ rqσ, we arrive at This polaritonic wave function as the ground state of (7) is the reformulation of the original electron-photon problem of (1) we were looking for. Since all the new photonic coordinates belong to harmonic oscillator ground states, exchanging p α,i with p α,j does not change the total wave function Ψ and this property transfers to the exchange of any coordinate q α,i and q α,j . Hence we have now a bosonic symmetry with respect to the q coordinates. Since the electronic part of the auxiliary system is not affected by the coordinate transformation, the electronic symmetries are the same in the physical and auxiliary system, i.e. we have a fermionic symmetry with respect to rσ. Together these two fundamental symmetries imply that the polaritonic coordinates zσ have fermionic character. The symmetries of the polaritonic wave function Ψ can be summarized as This means that though the dressed wave function has fermionic statistics (15) in terms of the polaritonic coordinates zσ, due to the constitutive relations (14) it actually consists of two types of particles: one with fermionic character and another with bosonic character.
Consequently, the polariton wave function Ψ has a hybrid Fermi-Bose statistics. As consequences of these symmetries we find the Pauli exclusion principle for the electrons, yet for the auxiliary photon coordinates we find that many photonic auxiliary entities can occupy the same quantum state.
Let us briefly highlight that the constitutive relations and the hybrid fermion-boson symmetry are fundamental to the polaritonic wave function. The fermionic statistics of Eq. (15), which are merely a consequence of the constitutive relations, are in general not sufficient to guarantee a physically reasonable result. Let us illustrate this with a simple two-electron-one-mode example with no electron-electron and electron-photon interactions.
The exact physical ground state is just a product of a Slater determinant consisting of the two lowest eigenfunctions ψ 1 , ψ 2 of and χ 0 (p) = ( ω π ) 1/4 e −ωp 2 /4 , which is the ground state of a harmonic oscillator with frequency ω. We obtain the dressed version of Ψ by multiplication with another oscillator ground state χ 0 (p 2 ) with the same frequency. Performing the coordinate transformation (6) in this case simply replaces the coordinates (p, p 2 ) with (q 1 , q 2 ) since they are orthogonal. 4 Note that for a coupled problem the transformation becomes more involved (see Appendix A 4 Specifically, we need to perform the transformation for the term χ 0 (p)χ 0 (p 2 ) ∝ e −ωp 2 /2 e −ωp 2 2 /2 . Since an orthogonal transformation for a set of variables (x 1 , x 2 , ...) → (x 1 , x 2 , ...) leaves any expression of the form for details). The correct auxiliary ground state reads where in the last line we subsumed the electronic and photonic orbitals with the same coordinate index to a polariton orbital, i.e. ψ i (r j σ j )χ n (q j ) ≡ φ in (z j σ j ). This wave function is obviously antisymmetric with respect to the exchange of z 1 σ 1 and z 2 σ 2 but it also obeys the constitutive relations that enforce the hybrid symmetry in rσ and q, respectively. Now let us consider a possible wave function that only obeys the overall fermionic symmetry, e.g. Only , where 1 and 2 are the eigenenergies corresponding to ψ 1 and ψ 2 , a simple minimization would yield the right symmetry solution. If, however, Indeed, if we enforce the constitutive relations also onΨ by adding two extra terms, we see Thus we obtain the desired result of ruling out the solution that violates the Pauli principle.
It can be shown 35 that enforcing the conditions of Eq. (14) is sufficient for the dressed system to attain as ground state Ψ = Ψχ, where Ψ is the ground state of the original Hamiltonian of Eq. (1) and χ the product of ground-state harmonic oscillators. For excited states, the constitutive relations are necessary but not sufficient to single out the eigenfunctions ofĤ that correspond to the original Hamiltonian in terms of simple products. 34 So, how can we enforce the right symmetry and with this the right hybrid statistics in practice? We would like to use the polaritonic orbitals φ i (zσ) as our basic physical entity and this makes the constitutive relations quite special. This is due to the fact that they concern always just a subset of the coordinates of one orbital. Contrary to the usual approach in many-body physics, where a single-orbital basis is chosen and the fermionic/bosonic symmetry is enforced by building Slater determinants/permanents. One could try to transfer this approach to the polariton case and construct a many-body basis with elements that satisfy Eq. (14) and thus are Slater determinants not only with respect to the polaritonic, but also with respect to the electronic coordinates (or equivalently permanents with respect to the photonic coordinates). This would lead to basis elements with "mixed-index" orbitals, i.e. orbitals φ(r i , q j , σ i ) that depend on coordinates with different indices i = j. If we want to calculate the expectation value of an observable and integrate over all coordinates, we see that coordinate i and j are coupled alone due to the "mixed-index" orbital. This has severe consequences. For example, we cannot make use of the orthonormality relation to set certain terms to zero when we calculate matrix elements and one-body terms become "two-body like." The number of such anomalous terms grows factorially with the particle number, which makes such a type of ansatz infeasible in practice. For details the reader is referred to Appendix B.
As an alternative to looking at the polaritonic wave function directly, the physical conditions (14) are visible in the dressed one-body reduced density matrix (1RDM), 35 which is given explicitly by The Fermi statistics of the wave function Ψ with respect to the polaritonic coordinates zσ (15) is also apparent in γ[Ψ ] in the form of so-called N -representability conditions. 47 By using the natural orbitals φ i and the natural occupation numbers n i , which are defined by the eigenvalue equation n i φ i =γ ψ i , we represent γ in its diagonal form The fermionic N -representability conditions become especially simple in this representation and are given by In general, it can be proven 47 that any matrix that fulfils the conditions (18) is connected to an ensemble of N -body states that are fermionic with respect to the exchange of their coordinates.
From the dressed 1RDM we can define the electronic 1RDM and the auxiliary photonic 1RDM Again, we can define the according natural orbitals ψ and The Fermi statistics with regard to only the electronic coordinates rσ thus becomes apparent by considering the electronic natural occupation numbers n e i n e i ≥ 0, ∀i where we split the conditions in three parts for later convenience. The equivalent bosonic symmetry of the auxiliary photonic coordinates leads instead to the conditions Note that the normalization of γ e/b to the electron number N is a direct consequence of the auxiliary construction that considers exactly N polaritons for a system with N electrons.
This becomes explicitly visible in the fact that the normalization of γ by definition transfers ). Additionally, the lower bounds of γ e/b , c.f. Eqs. (23a) and (24a), transfer from γ, because the partial trace operation is a completely positive map. 48 We can conclude that if (18)  This now shows explicitly also for an interacting wave function that at most one electron can occupy a specific quantum state, while many auxiliary photon quantities can occupy a single quantum state. Further, the dressed 1RDM γ[Ψ ] itself has only natural occupation numbers between zero and one and is therefore fermionic, yet it contains a fermionic and a bosonic subsystem. It is important to note that the original wave function Ψ did not have this simple hybrid statistics but only fermionic symmetry, since the physical p α did not follow any specific statistics. Further, that we genuinely have formulated the coupled electron-photon problem in terms of hybrid quasi-particles becomes most evident by actually using single-particle (polariton) orbitals φ i (zσ) to expand the dressed 1RDM of Ψ . Why we consider the 1RDM and how we can do this efficiently will be discussed in the following.
Let us illustrate this with the example from before. The correct auxiliary ground state Ψ satisfies the conditions of Eq. (23), since The two electronic orbitals are the eigenfunctions (natural orbitals) of γ e [Ψ ] with natural occupation numbers n 1 = n 2 = 1. If we do the same calculation withΨ , we get instead which violates the N -representability conditions (23)  Nevertheless, the conditions (23) are sufficient to guarantee that an ensemble (mixed state) Γ = j w j |Ψ j Ψ j | with j w j = 1 exists, which yields this 1RDM 47 . The 1RDM is then said to be ensemble N -representable. Thus, though (23) does not guarantee that Ψ satisfies (14), it does guarantee that it satisfies the Pauli exclusion principle and additionally that an ensemble constructed from |Ψ j which do satisfy (14) exists.
To obtain a computationally tractable procedure, we therefore use the construction presented in the next section to ensure the polariton statistics implied by Eq. (14), instead of the factorially growing number of "mixed-index" orbitals. We will consider all fermionic density matrices in the auxiliary configuration space, which we characterize by the conditions of Eq. (18) in terms of polaritonic orbitals φ i (zσ). We then constrain this space by enforcing the N -representability conditions of Eq. (23) for the 1RDM of the electronic subsystem.
Since this guarantees that only (ensembles of) fermionic wave functions are allowed, also the minimal energy solution has fermionic symmetry with respect to rσ. This together with the zσ antisymmetry implies that q has bosonic symmetry. We call this construction the polariton ansatz for strong light-matter interaction. In the next section, we will based on the polariton ansatz provide a detailed prescription to generalize a given electronic-structure theory to treat ground states of coupled electron-photon systems from first principles.
3 First-principle theories: from electronic to polari-

tonic bases
In this section, we lay out in detail how one can transform a given electronic-structure theory that meets some minimal requirements into its polaritonic version. The goal of such a "polaritonic-structure theory" is to find the ground state of the Hamiltonian of Eq. (1) by considering the ground state of the auxiliary Hamiltonian of Eq. (7). We define the according variational principle for the ground-state energy E 0 as where P = {Ψ : Ψ ↔ (14)} is the set of all normalized many-polariton wave functions that obey the constitutive relations of Eq. (14). For our purposes, as explained in Sec. 2, we will instead consider the larger set M of all (mixed-state) density matrices Γ = j w j |Ψ j Ψ j | with j w j = 1, that obey the hybrid Fermi-Bose statistics. The minimal energy also in this more general set corresponds to the pure state of Eq. (25), i.e.
The main trick now is in how we construct this set. We do so by first considering the yet larger setM = {Γ :Ψ j = C j,K Φ K }, i.e. density matrices made of superpositions of Slater determinants Φ K = det(φ K,1 · · · φ K,N )/ √ N ! of polariton orbitals φ K,i . This guarantees the overall Fermi statistics in terms of the polaritonic coordinates zσ. We then constrain this larger set to where  18)). As we explained in Sec. 2, this is sufficient to ensure the constitutive relation (14). However, this does not automatically imply that the wave functionsΨ j that are used to construct the constrained density matrices obey the basic symmetries. Rather, they are also auxiliary quantities and it is only the density matrices that are the physical objects. Due to the ensemble conditions there are, however, wave functions with the exact hybrid symmetries associated. We thus avoid the direct construction of the exponentially growing correlated electron-photon states. 5 Most importantly, the polariton picture gives any coupled problem of the form of Eq. (1) the same structure as a purely electronic problem with two-body interactions. Consequently, we can transfer every type of electronic-structure theory to the coupled electron-photon problem, if the theory provides an expression for the 1RDM (since we need the 1RDM to 5 Strictly speaking, we should construct the ensembles Γ = j w j |Ψ j Ψ j | which generates γ[Ψ ] to evaluate the energy as Tr{ΓĤ} to remain variational (26). But instead we evaluate the energy directly from Ψ which only satisfies (23) as Ψ |Ĥ Ψ . However, because γ(zσ, z σ ) is correct, the error is only in the correlation part of the two-body part of the energy. This means that for electronic structure theories only based on the 1RDM, e.g. Kohn-Sham DFT, HF and RDMFT, it is exact, since only the 1RDM is relevant in those theories. test the N-representability constraints). The main steps how to do so are depicted in Fig. 3.

Hamiltonian
Basis States Lagrangian of the EST More specifically, the connection betweenĤ m and E m is given by the particle number N Typically, one needs to impose some constraints on the parametrization of the wave function Ψ to make it physical c k [Ψ] = 0, e.g. orthonormality of orbitals or the norm of the CI coefficients. The generic electronic-structure minimization problem is formulated as subject to c k [Ψ] = 0.
We can solve Eq. (28) by, e.g. minimizing the Lagrangian where is a Lagrange-multiplier term. Instead of minimizing E m directly, one minimizes L m with respect to the orbitals and the Lagrange-multipliers k .
Today, a plethora of standard electronic-structure codes exist that solve (28) very efficiently for many different theory levels and thus allow for a highly accurate description of electronic structure.
If we consider the coupled electron-photon Hamiltonian of Eq. (1) instead of the purely electron HamiltonianĤ m , we find that we need to build new approximation strategies and implementations to deal with the coupled electron-photon Hamiltonian directly. However, by transforming the problem into its dressed counterpart, i.e. we consider (7), we can utilize the full existing machinery for the electronic case. In particular, this means that we have now Since there are many possible strategies to solve the minimization problem (31) and a good choice might depend on the specific electronic-structure theory that is considered, we finish the general discussion here. Instead, we conclude the section with a concrete example. We will apply the above rules to HF theory, which leads to polaritonic HF theory. This means that we approximate the density matrix of the exact dressed wave function of Eq. (13) by the density matrix of a single Slater determinant with orbitals φ 1 , . . . , φ N , i.e.
where P N denotes the permutation group on N elements and the index j is chosen such that it is even (odd) for an even (odd) permutation π i ∈ S N . Further, we consider a spinrestricted formalism, i.e. we assume that the number of electrons N is even and define φ 2k−1 (zσ) = φ k (z)α(σ), φ 2k (zσ) = φ k (z)β(σ) for k = 1, . . . , N/2, where α, β are the usual spin-orbitals. We again note that we do not necessarily enforce with our constraints that the auxiliary Slater determinant has the right symmetry but rather its 1RDM. In this regard polaritonic HF becomes actually a 1RDM functional theory for polaritonic problems rather than a wave-function based method. 35 With this ansatz, we calculate the energy expectation value for the Hamiltonian of Eq. (7), which reads where we introduced the "dressed" Coulomb-operatorĴ i which acts aŝ and the "dressed" exchange-operatorK i that acts as 50 The polaritonic one-and two-body terms are given by (8), (9) and (10), respectively. With this we find that E HF = E HF (t , v , w , {φ k }). Consequently, we also find structurally the same derivative (that for HF theory is called the Fock-matrix), which reads Since we consider only one Slater determinant, the orbitals φ i are also the eigenfunctions of the system's dressed 1RDM γ. Because of the spin-restriction, it suffices also to consider the spin-summed version γ(z, z ) = 2 N/2 i=1 φ i * (z )φ i (z), which we denote with the same symbol. We see that γ[Φ ] has occupations (eigenvalues) of 2 instead of 1 because of the spin-summation. This transfers to the natural occupation numbers of the electronic 1RDM γ e (r, r ) = dqγ(rq, r q). Now, we have defined all the terms that enter the minimization problem (31) and we can proceed with demonstrating how to solve it. Algorithmic-wise, we are confronted with enforcing the additional inequality constraints (30) in extension to the original (HF) minimization problem in (28). We start by noting that the constraint functions g i depend on γ{φ k }, which can be directly calculated from the polariton orbitals, via the eigendecomposition of γ e . Since the diagonalization of γ e is a non-trivial step for large systems (or in real-space) and thus can be a bottleneck of the minimization, it is helpful to consider natural and dressed orbitals as independent variables of the minimization and enforce their connection as an additional constraint. 7 We thus define g i = g i [γ e {φ k , ψ e i }] and include the necessary orthonormality of the ψ e i by a third set of conditions that we include in the minimization by a third Lagrange-multiplier term − ijθ ij f ij . Note that this construction automatically linearizes the constraints (30) during one minimization step, where the ψ e i are fixed.
To enforce now these inequality constraints, we use an augmented Lagrangian algorithm, following the book of Nocedal and Wright 51 , Ch. 17.3. 8 We chose this algorithm, since it simply extends a given Lagrangian with penalty terms. Hence, we can make use of any existing implementation that solves the minimization problem of Eq. (28) (30) quadratically, but has no effect in the so-called feasible region of configuration space, where the conditions (30) are satisfied.
Specifically for our example, the extra Lagrangian term of the translation rules depicted in Fig. 3 is given by The full Lagrangian for the polaritonic HF minimization problem reads then and the corresponding first order conditions for a minimum (stationary point) of L HF are where we considered φ k and φ k * as independent and definedĜ i φ k (rq) = n e k dr ψ e i * (r )φ k (r q) ψ e i (rσ). Additionally, we can diagonalize the Lagrange-multiplier matrices¯ ij = δ ij j and  (32) and (34). We did this for a model Hamiltonian and will present the results in the next section.
Before we come to the results, we want to conclude this section with a brief remark on how to generalize the just presented example of polaritonic HF theory. Since we consider in HF theory only one Slater determinant, the equations become especially simple. However, including (polaritonic) correlation beyond this approximation is easily possible by, for instance, polaritonic QEDFT. To do so, we can leave the structure of the implementation exactly as it is, but just exchange the operatorsK i in (32) and (34)   However, for the exemplification, we consider a one-dimensional lattice system that couples to one photon mode in dipole approximation with frequency ω (a sketch of the setup has been depicted in Fig. 1). We want to remark that including more modes is in principle possible (see Sec. 3), but the one-mode case is the standard. The Hamiltonian is of the form of Eq. (1) and readŝ with hopping t = 1 2∆x 2 corresponding to a second-order finite difference approximation for a grid with spacing ∆x, where we choose t = 0.5 for all calculations, which corresponds to a spacing between neighbouring sites of ∆x = 1 bohr and a local scalar potential with value v i on site i. We have set the Coulomb repulsion to zero in this example to highlight the influence of the matter-photon coupling and how well the polaritonic HF approach can capture it. Nevertheless, due to the dipole self-energy termĤ d we have a mode-induced dipole-dipole interaction among the electrons. This type of interaction is important in many fundamental quantum-optical questions, such as the quest for a super-radiant phase in the strong-coupling case. [54][55][56] Further, the electron basis B m is determined by the number of sites, x i = i − x 0 is the position with respect to the middle of our lattice x 0 ,ĉ ( †) i,σ are 9 A similar consideration holds also for a real-space description, where the Fock-matrix cannot be explicitly constructed, but is diagonalized iteratively by, e.g., a conjugate-gradient algorithm like the one of our implementation. 52 The scaling of such a method in the electronic case is of the order O(B m ln B m N 2 ) and can even be reduced with state-of-the-art algorithms to O(B m ln B m N ), 53 where N is the number of electrons. Such methods thus scale better than a direct diagonalization of the Fock-matrix, but at the same time the underlying basis size B m describes the number of grid points and thus is typically much larger than B m in orbital-based codes. However, for very large systems, both descriptions can have comparable B m . A real-space like description of the displacement coordinates of the photon modes would in a similar manner be quite inefficient for a few photons, but might become even advantageous for large photon numbers. the fermionic creation (annihilation) operators that satisfy the anticommutation relation [ĉ i,σ ,ĉ † j,σ ] + = δ ij δ σσ , andn i =ĉ † i,↑ĉ i,↑ +ĉ † i,↓ĉ i,↓ is the density operator. We have chosen a simple lattice model, since this allows us to have still exact numerical reference data to compare to. In contrast to standard electronic-structure problems, there are are no (numerically exact) references solutions currently available for real three-dimensional systems in the continuum. To the best of our knowledge there are only QEDFT simulations (at several levels of approximations) available for realistic three-dimensional systems. 57,58 For the implementation and to go to the polariton picture, we express the matterĤ m plus dipole partĤ d of our Hamiltonian in matrix form by using the basis states |ψ i,σ =ĉ † i,σ |0 . As a basis for the photon subsystem, we utilize the eigenstates χ i of the photon energy operator, i.e.Ĥ ph χ α = (α + 1/2)χ α , which are photon number states. To calculate the coupling term of the energy expression H I , c.f. Eq. (32), we express the displacement operator However, the interpretation of the operators is different from before, because we have to apply them to polaritonic basis states. Since we consider the spin-restricted formalism introduced in the end of Sec. 3, we neglect the spin-dependency of the electronic part of the basisψ i,σ → ψ i and define |φ iα = |ψ i χ α . We can then derive the kernel expression For the two-body term, a definition analogously to (40) is more difficult, since we have to differentiate the two polaritonic coordinates. For the sake of the analogy, we formally writê where the upper indices differentiate the two polaritonic orbitals that both have an electronic and photonic part. This is to be understood in the following sense: To define the the operators only act on the basis elements with the same indices. The kernel (w self ) i 1 α 1 i 2 α 2 j 1 β 1 j 2 β 2 for the self-interaction partŴ self = λ 2 Bm i 1 =j 2 =1n i 1n j 2 x i 1 x j 2 reads for example With these definitions, we can calculate the polaritonic HF energy expression, c.f. (32) and the polaritonic HF Fock-matrix, c.f. (34). Then, we employ the augmented Lagrangian algorithm as discussed in Sec. 3 to find the polaritonic HF ground state of the model system.
As a first example, we illustrate the violation of the Pauli principle if we do not enforce the right symmetries, as also discussed in Sec. 2 for a simple uncoupled problem. To this end, we compare ground-state energies, electronic 1RDMs and the photon number of a small 4-electron system obtained with the two different HF ground states, i.e. polaritonic HF using density matrices with the exact symmetry, c.f. (14), and polaritonic HF with only fermionic symmetry (which we call in this section fermionic HF). We can expect deviations between both polaritonic-HF theory levels for systems that contain more than one orbital. In our spin-restricted case this corresponds to more than two electrons and that is why we chose here N = 4. Further we set the external potential to zero, i.e. v i = 0 ∀i. Since we need to calculate the exact coupled electron-photon many-body ground state from a configuration space that grows exponentially fast with the size of the basis sets and the electron number, we choose a small box of length L = 5 bohr. This corresponds to a matter basis of B m = 6 spatial sites times two spin states for each electron. For the photon-subsystem, we consider B ph = 5 photon number states for which all relevant quantities are sufficiently converged. 10 Despite the small basis sets and electron number employed, the many-body configuration space has the considerable size of (2B m ) N * B ph ≈ 10 5 , which is already at the edge of standard exact diagonalization solvers: matrices of this size can still be diagonalized without special efforts like parallelization. Since we only aim for a benchmark study here, this limitation is not problematic, but it shows how expensive the exact solutions of coupled electronphoton systems computationally are. The need for numerically manageable approximations is evident here.
We first compare the electronic 1RDMs γ e , c.f. Eq. (19), and the photon numbers N ph = N ph , c.f. Eq. (3) using the connection formula of Eq. (11), for varying coupling strengths g/ω = λ/ √ 2ω and ω = 0.4. We choose these quantities because they provide us with a consistent measure of how well the electronic part and the photonic part of the system are approximated. The electronic 1RDM determines all electronic one-body quantities and is therefore a very precise measure of the quality of the approximation in the electronic sector.
Indeed, the photonic 1RDM in the single mode case corresponds to the photon number N ph , which is an equivalently good measure for the quality in the photonic sector. In Fig. 4 (a), we display the difference of the exact electronic 1RDM from the one of the polaritonic HF and fermionic HF approximations, measured by the Frobenius norm A 2 = ij A 2 ij for a matrix A ij . We see that for all coupling strengths the polaritonic HF 1RDM (dasheddotted orange line), which enforces the right hybrid statistics, remains very close to the exact solution indicating that the electronic subsystem is captured very well within this approximation. The fermionic HF (solid blue line) approximation, however, deviates strongly due to its wrong purely fermionic character. In both cases, fermionic HF deviates much stronger from the exact reference than polaritonic HF due to the wrong symmetry.
The same behavior is also encountered in the photonic subsector, where in Fig. 4  We therefore find, similarly to the simple uncoupled problem in Sec. 2 (for g/ω = 0 we recover this case exactly), that for the same energy expression using the wrong statistics leads to sizeable errors.
The same problem is encountered also in Fig. 5 (a), where we display the total energy E = Ĥ of the coupled system as a function of the coupling strength g/ω. While the polaritonic HF (dashed-dotted orange line) is variational, i.e. due to the right statistics we are always equal or above the exact energy (dashed green line), the fermionic HF (blue solid line) breaks the proper symmetry and thus can reach energies below the physically accessible ones. However, again in close analogy to the uncoupled example in Sec. 2, if we increase the frequency of the photon field such that it is much more costly to excite photons than electrons, the minimal-energy conditions can single out the correct statistics, as displayed in Fig. 5(b). That is, for ω large enough the constraints g i [γ{φ k , ψ e i }] ≥ 0 of Eq. (31) are trivially fulfilled.
(a) (b) Figure 5: Total energy for the 4-electron system with ω = 0.4 hartree (a) and ω = 0.8 hartree (b) for varying coupling strength g/ω. While in (a) the fermionic HF approximation (solid blue line) can achieve unphysically low energies when compared to the exact solution (dashed green line) due to the wrong statistics, in (b) the minimal-energy condition singles out the right statistics without further constraints. The polaritonic HF (dashed-dotted orange line) by construction always has the right hybrid statistics and thus is variational, i.e. the energy is always above the exact energy.
Having illustrated the need to enforce the correct symmetry for the use of the polaritonic basis, we employ polaritonic HF to study the effect of electron-photon coupling versus electron localization. We consider a matter system with a local potential v(x) = N/ √ which represents a potential well that is deep (shallow) for small (large) . The parameter thus represents the level of confinement of the potenial v(x) which is depicted in green (for various values of ) in Fig. 6. To reduce boundary effects (which also represent a form of confinement), we consider a box length of 30 bohr (B m = 30), corresponding to a real-space grid from x = −14.5 bohr to x = 14.5 bohr. We consider the case of a 2-electron and a 4-electron system, respectively, and set ω = 0.1 hartree, which is far away from the regime where the fermionic HF approximation is valid. Thus the right hybrid statistics of the polaritons are crucial. We consider again B ph = 5, for which all the results are converged. We want to stress here that the corresponding many-body space for 4 particles has a dimension of (2B s ) N * B ph = 64.8 · 10 6 and thus is practically inaccessible by exact diagonalization.
Let us first consider how the electronic ground-state density changes when coupling and the localization are varied. To facilitate the comparison between the N = 2 and N = 4 case we plot in Fig. 6 the normalized electronic ground-state density ρ(x)/N , where ρ = γ e (x, x) is the diagonal of the electronic 1RDM (in blue) and the normalized confinement v(x)/N (in green). In Fig. 6 (a) we show the uncoupled 2-particle case and in (b) we use g/ω = 0.2 for the 2-particle case for varying . In Fig. 6 (c) and (d) we show the same plots for the 4-particle case. In both cases we see that for strongly-confined electrons, i.e. for small values of , the influence of the strong light-matter coupling on the density is negligible. This is in agreement with the usual assumption underlying, e.g. the Jaynes-Cummings model, that the ground state for atomic systems is only slightly affected by coupling to the photons of a cavity mode. Much higher coupling strengths would need to be employed in order to see a sizeable effect for strong localization. In contrast, once we lift the confinement and the electrons get delocalized, the influence of the light-matter coupling becomes appreciable.
The induced changes are not uniform but depend on the details of the electronic structure, i.e. in the N = 2 case we have a clear localization effect (Fig. 6 (b)) while for N = 4 we have an enhancement or even emergence of the double-peak structure (Fig. 6 (d)). That the changes in the electronic structure due to strong electron-photon coupling are non-trivial are in agreement with previous studies. 35,57 (a) Figure 6: Every plot depicts the normalized electron density ρ(x)/N (blue) with corresponding local potentials v(x)/N = 1/ √ x 2 − 2 (green, rescaled by a factor of 0.25 for better visibility) for a series of softening parameters of the 2-(upper row) and 4-electron (lower row) system and for coupling strength g/ω = 0 (left column) and g/ω = 0.2 (right column). We see how different both systems respond to the coupling depending on the degree of confinement, measured by . Note that the legends only depict three example lines (the smallest, largest and middle values of ), respectively.
To make these observations more quantitative we display in Fig. 7 the normalized changes in the electronic 1RDMs depending on the confinement and the coupling strength, i.e. ∆γ e = γ e g/ω, − γ e g/ω=0, 2 /N , in panel (a) and (d) for the 2-particle case and the 4-particle case, respectively. Also, we show the photon number N ph in the ground state in dependence of confinement and coupling strength in (b) and (e) for the 2-particle case and the 4-particle case, respectively. As a third quantity we consider ∆n e = i ||n e i,g/ω, − n e i,g/ω=0.0, || 2 , where n i are the natural occupation numbers. For the zero-coupling case they are all either zero or one, which corresponds to a single Slater determinant in the electronic subspace. If they are between zero and one they indicate a correlated (multi-determinantal) electronic state. Therefore ∆n e measures the photon-induced correlations and also highlights that although polaritonic HF is a single-determinant method in the polaritonic space, for the electronic system it is a correlated (multi-determinantal) method. For both, the 2-and the 4-particle case we find consistently that the more delocalized the uncoupled matter system is, the stronger the coupling modifies the ground state. Although this effect depends on the details of the electronic structure as we saw in Fig. 6, the plots of Fig. 7 indicate that this behavior is quite generic. The reason is that within a small energy range many states with different electronic configurations are available as opposed to a strongly bound (and hence energetically separated) ground-state wave function (see Fig. 8). A glance on the correlation measure ∆n e in panel (c) and (f) of Fig. 7 strengthens this explanation: For large and g/ω the electronic correlation is strongest and thus many electronic configurations contribute to the states of this parameter regime. This indicates also that the effective one-body description of the ground state of many cavity-QED models might be inaccurate in this regime. Additionally, we observe that the maximal values of ∆γ e in the 2-particle case are slightly larger than in the 4-particle case, which again is related to the different electronic structures of the two systems.
For the photon numbers, however, which are depicted in panel (b) and (d) of Fig. 7, we see that in general the number of photons is larger in the 4-particle case. This is due to the simple reason that the more charge we have the more photons are created. Nevertheless, the amount Figure 7: The plots show key quantities of the model system as a function of the coupling strength g/ω and the localization parameter for the 2-electron (upper row) and 4-electron (lower row) case. In the first column, we depict the normalized deviation ∆γ e = ||γ e g/ω, − γ e g/ω=0.0, || 2 /N of the electronic 1RDM to a reference for the same and g/ω = 0 (blue). In the second column, we show the total photon number N ph (green) and in the third column, the total deviation of the electronic natural occupation numbers ∆n e = i ||n e i,g/ω, −n e i,g/ω=0.0, || 2 is displayed. This is a measure of the induced electron-photon correlation. We observe in all the cases that when the bare matter wave function becomes more delocalized ( larger), the modifications due to the matter-photon coupling become stronger. of photons does not just double (as expected from a simple linear relation) but is almost three times higher. This highlights the non-linear regime of electron-photon coupling that we consider here. Again the number of photons increase also with the delocalization and hence the parameter is a very decisive quantity. All these results point towards an interesting parameter in the context of strong light-matter coupling: the localization of the matter wave function. In agreement with a recent case study for simple 2-particle problems, 41 systems that are less confined react much stronger to a cavity mode. This does not only hold for the ground state, 41 but suggests that if we want to observe genuine modifications of the ground state due to strong light-matter coupling, we should consider matter systems that have a spatially extended wave function. One way would be large molecular or solid-state systems, the other way would be an ensemble of emitters. In the latter case the strong influence would lead to local changes, in contrast to the Dicke-like description of collective strong coupling, where we have N independent replicas of the same, perfectly localized system. In both cases, it seems plausible that there are strong modifications of the electronic structure even if there is only coupling to the vacuum of a cavity. Modifications of chemistry by merely the vacuum seem therefore to be feasible. At the same time an interesting perspective with respect to collective strong-coupling arises: Maybe it is not the simple Dicke-type collectivity, 59 where an excitation is delocalized over many replica, that drives the changes in chemistry, but rather a genuine cavity-photon mediated spatial delocalization of the ensemble wave function. To answer this question, further investigations especially for realistic three-dimensional systems and ensembles including the Coulomb interaction have to be conducted. And polaritonic first-principle methods as introduced in this work seem specifically well-suited to answer these interesting and fundamental questions. 2 ∂ 2 x − v(r)]ψ i = e i ψ i for N = 2 are shown as a function of the confinement parameter . We see that the energies approach each other with increasing and thus decreasing confinement. This means that for a fixed coupling strength, the less confined the system is the more states are available in a small energy range for photon-induced modifications of the ground state.

Conclusion and Outlook
In this article, we have highlighted the influence of the hybrid Fermi-Bose statistics of the polariton wave function in the dressed construction. This dressed construction allows to use (matter-only) first-principle methods to investigate strong light-matter coupled systems.
We have provided a simple and general prescription on how to turn a given electronicstructure method into a polaritonic-structure method by introducing conditions in terms of the 1RDM to ensure the right hybrid statistics. We have given several examples where ignoring the correct Fermi-Bose statistics leads to unphysical results. As a specific example of a polaritonic-structure method we have transformed electronic HF theory to polaritonic HF theory and have demonstrated numerically that this approach stays accurate even for very strong coupling between light and matter. Finally, we have applied polaritonic HF to investigate the influence of electron localization in the context of strong electron-photon coupling and found for non-trivial examples that the more delocalized the uncoupled matter wave function is, the stronger it reacts to the modes of a cavity, which comes along with an increase of electronic correlation. This result indicates that there might be strong modifications of the ground state for spatially extended systems even by only coupling to the vacuum of a cavity. Specifically this raises the question whether an ensemble of emitters with a spatially extended wave function might show collective strong-coupling effects that also modify the local electronic structure, in contrast to the Dicke-type collective coupling picture that does not consider electronic correlation.
The presented results provide a relatively straightforward way to perform simulations for strong light-matter coupling based on well-established theoretical and numerical methods.
The main two practical bottlenecks are the increase of the dimension of the basic orbitals of the theory and the inclusion of the new constraints that enforce the right statistics. The former bottleneck will in most cases mean that one has instead of 3-dimensional orbitals now 4-dimensional ones, since it is usually a single effective mode that is considered, and is therefore not overly numerically expensive. It is crucial to realize here that, besides the necessarily larger orbital bases, the scaling of a method is not effected by the inclusion of the photon field 11 . The second bottleneck is instead challenging algorithmic-wise. As implemented for the test examples, one has to construct the full dressed 1RDM and then determine the electronic 1RDM as well as diagonalize it in each iteration step of the selfconsistency cycle. Further, for each orbital one needs an extra Lagrange multiplier in the augmented-Lagrangian method. These two things taken together as well as the non-linearity of the constraint lead to a quite slow convergence. Yet, there are different methods to enforce the new constraints and also the existing implementation can be made much more efficient.
To extend, for example, the existing implementation of the dressed construction in the real-space electronic-structure code Octopus 60 to enforce the hybrid statistics, one would not explicitly diagonalize the electronic 1RDM, which is numerically prohibitively expensive 11 Note that we do not want to underestimate the influence of an entirely new dimension in the problem. The given statement depends of course on the necessary size of the photon-basis. However, in the cases that we studied so far, we observed the converged results with photon bases that where of the order of the number of particles. but resort to semi-definite programming methods 61 that exploit the reformulation of the constraints in terms of positivity-conditions (see Sec. 3). 12 . In orbital-based codes, our in Sec. 3 presented algorithm could be well applicable, if not too large systems are considered and a direct diagonalization of the electronic 1RDM is feasible. Note that in contrast to the matter-description, the extra photon orbitals do not necessarily have to be constructed in the code. The analytically known structure of the photon-number states allows to calculate all necessary matrix elements analytically. We have exploited this also in our implementation (see Sec. 4).
We believe that such extensions, though they require a certain amount of work, are indeed worthwhile. The presented results for a simple yet still numerically challenging model system highlight that much is still to be understood and discovered in the context of polaritonic chemistry and physics. Despite the great success that cavity QED models have in explaining certain aspects of this emerging field, similar to condensed-matter physics, genuine firstprinciple methods seem to be necessary to get a detailed understanding. Many of the exciting experimental results, like the changes in chemical reactions, still are poorly understood and the basic mechanisms are yet to be discovered. To find explanations that are as unbiased as possible, as well as to make quantitative predictions calls for the availability of genuine firstprinciple methods for coupled matter-photon systems. At this point, we want to stress that we considered polaritonic HF just as a test-case to get a basic understanding of the hybrid statistics from a physical but especially from a numerical perspective. This understanding is the basis to implement the whole machinery in a standard electronic-structure code like Octopus. Once, this is accomplished, we will be able to perform polaritonic RDMFT and QEDFT calculations. Previous works 34,35 showed that both levels of theory seem very promising for an accurate description of the strong coupling between molecular systems and a cavity mode. The main open gap of these methods was a numerically feasible inclusion of 12 Note that independently of the precise algorithm, the bottleneck in a full real-space description of polaritonic orbitals is the extra dimension. However, with modern high-performance clusters, systems that are not too large should be manageable even in 4D. the hybrid statistics, which we closed with this work.
Let us finally also comment on possible applications of the above construction for other multi-species situations. It seems possible, provided we can define species with the same number of particles, that one can instead of working with complicated multi-species wave functions, work with the combined density matrices and enforce the ensemble representability conditions on the subsystems. This does not necessarily need any dressed construction. For instance, think about the Schrödinger equation for electrons and nuclei/ions. Assuming that we have one kind of nuclei/ions we could express the combined density matrix in terms of electron-nuclei/ion pairs. It seems interesting to investigate the above procedure also in the context of such cases. A The exact mapping between the physical and auxiliary wave function for two electrons and one mode The ground state of Hamiltonian (1) for two electrons and one photon mode is of the form Ψ(x 1 , x 2 , p) = i,j,α C α ij ψ i (x 1 )ψ j (x 2 )χ α (p) for a given electronic spin-spatial one-particle basis {ψ i (x)} (x = (rσ)) and a photonic basis {χ α (p)}, that we specifically choose as eigenfunctions of the photon HamiltonianĤ ph , defined in Sec. 1. We map Ψ to the auxiliary ground state Ψ (x 1 , x 2 , p, p 2 ) = i,j,α C α ij ψ i (x 1 )ψ j (x 2 )χ α (p) χ 0 (p 2 ), where χ 0 (p 2 ) is the vacuum state of the auxiliary harmonic oscillator, which by construction has the same frequency and thus is the lowest eigenstate ofĤ ph . Then, we perform an orthogonal coordinate transformation, c.f. Eq. (6), which for the two-particle case reads Ψ (x 1 , x 2 , p, p 2 ) −→ Ψ (x 1 , x 2 , p = 1/ √ 2(q 1 + q 2 ), p 2 = 1/ √ 2(q 1 − q 2 )).
In this expression, only the first line is what would appear in a standard Slater determinant, i.e. overlaps of orthonormal orbitals with a constant norm for all φ a/b . However, with the terms that stem from the additional symmetry requirements, many new "mixed-index" terms arise, whenever one or more of the according orbitals have coordinates with different indices. All these 8 terms are two-body like integrals, which is in stark contrast to the calculation of the norm for a normal Slater determinant that involves only one-body terms.
Their computation is non-trivial, which also is contrast the normal terms of the first line, which are analytically given by the orthonormalization condition. The occurrence of such terms also means that the norm of Ψ ab depends on the specific form of φ a/b and thus has be calculated explicitly to normalize Ψ ab . The number of such terms for an N-body generalized determinant is given by the possible permutations of polaritonic and electronic coordinates N ! 2 minus the N ! "ordinary" terms and grows factorial with the number of particles, i.e.
N ! 2 − N !. We see that the normalization and in the same way also the calculation of expectation values of a wave function that explicitly exhibits the symmetry (14) requires the numerical calculation of (over)exponentially many non-trivial terms. This explicit ansatz is thus infeasible in practice and instead one should use more efficient ways to enforce the polariton symmetry as discussed in section 2.

C Numerical Details
For the numerical results of Sec. 4, we wrote a Python code relying mainly on the routines of the package NumPy. 13 The code specifically constructs the one-body Hamiltonian, defined in Eq. (34) allowing in principle for an arbitrary lattice basis for the matter system (in this publication we always considered a one-dimensional real-space grid, but one could also consider, e.g. atomic-orbital and implement the respective matrix elements) and Focknumber states for the one photon-mode. The minimization routine of the code is based on the conjugate-gradient algorithm described in Payne et al. 52 , Chap. V, replacing the Lagrangian and the gradient expression by (37)  calculations the overall convergence criterium of max(||(38a)||, |E m − E m−1 |) < 10 −4 , where E m is the total energy of the m-th iteration of the outer loop.