Ab-initio non-relativistic quantum electrodynamics: Bridging quantum chemistry and quantum optics from weak to strong coupling

By applying the Born-Huang expansion, originally developed for coupled nucleus-electron systems, to the full nucleus-electron-photon Hamiltonian of non-relativistic quantum electrodynamics (QED) in the long-wavelength approximation, we deduce an exact set of coupled equations for electrons on photonic energy surfaces and the nuclei on the resulting polaritonic energy surfaces. This theory describes seamlessly many-body interactions between nuclei, electrons and photons including the quantum fluctuation of the electromagnetic field and provides a proper first-principle framework to describe QED-chemistry phenomena. Since the photonic surfaces and the corresponding non-adiabatic coupling elements can be solved analytically, the resulting expansion can be brought into a compact form which allows us to analyze aspects of coupled nucleus-electron-photon systems in a simple and intuitive manner. Furthermore, we discuss structural differences between the exact quantum treatment and Floquet theory and show how existing implementations of Floquet theory can be adjusted to adhere to QED. From this generalized Born-Huang expansion an adapted Born-Oppenheimer approximation for nuclei on polaritonic surfaces can be deduced. This form allows a direct application of first-principle methods of quantum chemistry such as coupled-cluster or configuration interaction approaches to QED chemistry. By restricting the basis set of this generalized Born-Oppenheimer approximation we furthermore bridge quantum chemistry and quantum optics by recovering simple models of coupled matter-photon systems employed in quantum optics and polaritonic chemistry. We finally highlight numerically that simple few-level models can lead to physically wrong predictions, even in weak-coupling regimes, and show how the presented derivations from first principles help to check and derive physically reliable simplified models.

By applying the Born-Huang expansion, originally developed for coupled nucleus-electron systems, to the full nucleus-electron-photon Hamiltonian of non-relativistic quantum electrodynamics (QED) in the long-wavelength approximation, we deduce an exact set of coupled equations for electrons on photonic energy surfaces and the nuclei on the resulting polaritonic energy surfaces. Since the photonic surfaces and the corresponding non-adiabatic coupling elements can be solved analytically, the resulting expansion can be brought into a compact form which allows us to analyze certain aspects of coupled nucleus-electron-photon systems in a simple and intuitive manner. Furthermore, we discuss structural differences between the exact quantum treatment and Floquet theory, that along this line existing implementations can be adjusted to incorporate QED and how standard drawbacks of Floquet theory can be overcome. We then highlight, by assuming that the relevant photonic frequencies of a prototypical cavity QED experiment are in the energy range of the electrons, how from this generalized Born-Huang expansion an adopted Born-Oppenheimer approximation for nuclei on polaritonic surfaces can be deduced. By restricting the basis set of this generalized Born-Oppenheimer approximation we span the ark from quantum chemistry to quantum optics by recovering simple models of coupled matter-photon systems employed in quantum optics and polaritonic chemistry. We finally highlight numerically that simple few level models can lead to physically wrong predictions, even in the weak coupling regime, and show how the presented derivations from first principles help to check and derive physically grounded simplified models.

I. INTRODUCTION
In the last decade tremendous experimental advances have allowed to investigate and control complex manybody systems strongly coupled to photons [1][2][3]. In such situations of strong light-matter interactions novel physical effects can be observed such as symmetry protected collisions of strongly interacting photons [4], Bose-Einstein condensation [5] and room-temperature polariton lasing [6] of exciton-polaritons, or even the control of the energy-levels in living bacteria [7]. Such dramatic changes of physical and chemical properties of many-body systems can even be observed if no real photons are present and it is only the vacuum of, for example, an optical cavity that the matter couples to. Examples are a change of chemical reactivity under strong coupling to the vacuum electromagnetic field [8], different transition states in gas phase and cavity [9] and multiple-Rabi splittings under vibrational strong coupling [10]. Such experimental results highlight that disregarding the photonic degrees of freedom when calculating chemical and physical properties of many-body systems (as usually done in quantum chemistry and solid-state physics) can become inadequate when we change the bare electromagnetic vacuum, for example, by an optical cavity or nanoplasmonic devices, especially for strong light-matter coupling. In such cases we have to take into account electronic, nuclear and photonic degrees of freedom at the same time [11]. While the resulting equations that one would need to solve in principle are well-known [11][12][13], they are numerically unfeasible in terms of full many-body wavefunctions. One way to make these equations numerically tractable is by extending methods of many-body theory such as density/current-functional theory [14,15] or Green's function methods [16,17] to coupled matter-photon systems [13,[18][19][20][21][22]. While first ab-initio results for coupled matter-photon problems are available [23][24][25] so far most calculations have used simplified descriptions based on quantum-optical models like the Dicke or Tavis-Cummings model [26,27]. Although these models can be derived under certain assumptions from higherlying theories [28] and are able to reproduce well certain aspects of the experimental results [29,30], their limitations are as of yet not well explored. To investigate these possible limitations as well as to provide a consistent way of how to improve shortcomings is especially timely, since such models have been used to predict interesting new effects [31,32] in the context of the emerging field of polaritonic chemistry.
Furthermore, scrutinizing these quantum-optical models from a quantum-chemical perspective could also help to analyze long-standing problems of quantum optics such as the realizability of a super-radiant phase transition as predicted by the Dicke model [33,34]. In this work we investigate implications to electronic structure methods of working with the complete Hamiltonian of non-relativistic quantum-electrodynamics (QED) [11][12][13] in the long-wavelength approximation. Performing the Born-Huang expansion [35,36] of the coupled nucleus-electron-photon wavefunction allows us to exactly rewrite the eigenvalue problem as a set of arXiv:1804.00923v1 [quant-ph] 3 Apr 2018 coupled equations for nuclear, electronic and photonic degrees of freedom. In contrast to previous approaches we do not combine the photonic degrees of freedom with the nuclei [37] but group them with the electrons. While this gives rise to photonic potential-energy surface and non-adiabatic couplings for the electrons and thus also changes the usual electronic surfaces to polaritonic (electron-photon quasi-particles) surfaces for the nuclei, the photonic part can be solved analytically with the help of a shifted harmonic oscillator basis. This allows us to rewrite the Born-Huang expansion for coupled nucleus-electron-photon systems in a compact form. In this form the importance of the often ignored photon-mediated dipole self-energy term [38] becomes evident and we show how photonic observables can be constructed from the Born-Huang expansion. This expansion, although structurally reminisced of Floquet theory, illustratively presents the differences to the exact QED treatment. Nevertheless, we present how existing implementations can be adjusted to incorporate QED effects [39] and how standard drawbacks of Floquet theory can be overcome along this line. Then, by assuming that the kinetic-energy contributions of the nuclei are small and that the frequency of the relevant photon modes is in the energy range of the electronic excitations, for example, due to an optical high-Q cavity, we deduce an adopted Born-Oppenheimer approximation (the explicit polariton approximation) which is similar to the approach presented in Ref. [28]. By further restricting the basis expansion for the electronic and the photonic subsystem we then arrive at a simple few-level approximation that resembles often employed model Hamiltonians. We then scrutinize the resulting few-level approximations for a numerically exactly solvable electron-photon problem and highlight how in the weak-coupling regime some integrated quantities like the energies are well reproduced but photon occupation and real-space quantities like the density can be qualitatively wrong. By including only a few more states the results improve considerably provided one increases the electronic and photonic basis sets consistently. In the ultra-strong coupling regime, however, many more states need to be included for the integrated quantities to become accurate. This highlights that for strongly coupled problems the bare, that is, uncoupled, basis expansion is not very efficient and results based on a bare description, for example, using the standard Born-Oppenheimer states and dipole-coupling elements, become less trustworthy. This paper is structured as follows: In Sec. II we introduce the basic Hamiltonian, present the Born-Huang expansion and show the analytic solution of the photonic subsystem. In Sec III we then discuss the physical implications of the Born-Huang expansion and derive the explicit polariton approximation, which is a generalization of the quantum-chemical Born-Oppenheimer approximation. In Sec. IV we derive the small polariton approximation with the help of a bare basis expansion and connect to quantum-optical models. Next, in Sec. V we then give numerical details and present the accuracy of few-level approximations for the case of a GaAs quantum ring. Finally, we conclude and give an outlook in Sec. VI. In the appendix we further discuss how the coupling strength effectively changes in collective systems and show how the resulting framework can also be used to consider periodic or time-dependent systems.

II. THEORY
We start by presenting the basic QED non-relativistic Hamiltonian for nuclei, electrons and photons that we consider. Then we perform the Born-Huang expansion in terms of nuclear, conditional electronic and photonic wavefunctions. Since the Hamiltonian is spin independent (we ignore the fine-structure influence by spin-orbit and spin to magnetic field interactions) we can discard for notational simplicity and without loss of generality the spin degrees of freedom of the wavefunctions. By properly symmetrizing the wavefunctions at the end the physical eigenstates can be found. We further use atomic units throughout the work.

A. Pauli-Fierz Hamiltonian in the long-wavelength limit
Lets assume that for a coupled matter-photon system the relevant photon modes have wavelengths that are large compared to the typical size of a matter subsystem, for example, a molecule. This typically happens, for instance, if optical or microwave modes of the vacuum field are enhanced by enclosing the matter subsystem in a cavity. In this case an approximation of the full Pauli-Fierz Hamiltonian of non-relativistic QED in Coulomb gauge [11][12][13], where we neglect the spatial dependence of the photon modes in the length form [23,24,40] H =Ĥ n +Ĥ e +Ĥ ne +Ĥ p +Ĥ ep +Ĥ np .
is expected to be accurate. Here the nuclear Hamiltonian for N n nucleî includes similarly the corresponding sum over electronic kinetic components and the Coulomb electron-electron interaction. The nuclear-electron interaction is given accordingly byĤ Further, the photonic contribution for M p modes is then given bŷ which incorporates the total dipoleR = Ne j=1r j − Nn j=1 Z jRj of electrons and nuclei [41]. The photonmodes M p take into account the most important modes (see also Subsec. III D) but can in principle run from the fundamental resonator mode up to a maximum sensible frequency, for example an ultra-violet cut-off at rest mass energy of the electrons. The quantized oscillators representing the photonic system consists of the canonical coordinate corresponding to the displacement fieldq α = 1 √ 2ωα (â † α +â α ) and its conjugate momen- [23,24] with [q α ,p α ] = iδ αα . The fundamental coupling strength λ α = λ α e α describes the coupling between the total dipole and the photonic mode α with wavevector k α and transversal polarization vector e α . Here the coupling strength depends on the form of the mode functions S α (r), that is the geometry of the cavity and the chosen reference point for our matter subsystem [13,23,24]. This coupling is an inherent feature of the physical set-up, for example, the form and nature of the cavity, and cannot directly be amplified by the number of charged particles. We discuss this in some more detail in App. A and highlight how the number of emitters can, however, enhance the effective coupling strength when a reduced description is employed. The numerical investigations in Sec. V then refer to quite strong effective couplings, approachable for example in circuit QED [3,42].
The self-polarization part in the photonic Hamiltonian 1 2 (λ α ·R) 2 naturally arises in the length form to make the Hamiltonian bounded from below, which is a prerequisite to allow for ground states of an interacting cavity manybody system [38]. In the Subsec. III B we give a very intuitive physical picture for this abstract statement. The self-polarization further renders the Hamiltonian invariant under translations provided we have a charge neutral system, that is, Nn j=1 Z j = N e . In the case that the total system is not charge neutral, the center-of-charge couples to the photonic field and for eigenstates a translation byR , namelyR →R +R, leads to a trivial elongation of the photonic displacement byq α →q α − λα ωα ·R (see also Sec. II C). For the above nucleus-electron-photon Hamiltonian we then want to determine the eigenfunctions [43] HΨ i (R n , r, q) = E i Ψ i (R n , r, q) . (3) with the many-body energies E i , where R n , r and q are collective variables defined as R n = (R 1 n , R 2 n , ..., R Nn n ), r = (r 1 , r 2 , ..., r Ne ) and q = (q 1 , q 2 , ..., q Mp ). Next we perform a Born-Huang expansion where we expand into subsystem wavefunctions. This expansion can be performed in multiple different ways resulting in alternative physical interpretations and consequences for approximations. We will elaborate on their relevance and implications a little later. Here we use χ µ i (R n ) that represent nuclear wavefunctions andΨ µ ({R n }, r, q) the polaritonic components describing the correlated electronphoton system. We chooseΨ µ such that they form an orthonormal basis in the electron-photon subsystem and assume that they depend in an yet unspecified parametric way on the positions of the nuclei {R n }. Due to this parametric dependence we will later find equations that couple the different subsystem wavefunctions. In a second step, we further expand the polaritonic wavefunction into electronic ψ k ν (r, {R n }) and photonic Φ k (q, {R}) subsystem wavefunctions We use here that in the photonic subspace we employ an orthonormal basis of wavefunctions Φ k (q, {R}) that parametrically depend on the total dipole of the matter (nuclei and electrons) subsystem, that is, Φ k |Φ k p = dqΦ * k (q, {R})Φ k (q, {R}) = δ kk , as well as electronic wavefunctions parametrically dependent on the position of the nuclei such that This allows to re-express the normalization as Note that by the introduced expansion we find for each polaritonic eigenstateΨ µ electronic states ψ k µ associated to a photonic excitation Φ k . This shows that the electronic space is repeated with associated photonic excitations. We can see here a similarity to Floquet-theory [44] (see Sec. III D for details) where a matter system driven by a classical external field is considered. In this case, by assuming periodic driving, a time-dependent problem can be rewritten as an eigenvalue problem in a Hilbert space including time. The resulting eigenvalue equation is unbounded from below which expresses itself by the Floquet block index l ∈ Z that is usually interpreted (despite the absence of actual photons) as the number of photons involved in the process. The Floquet approach therefore has no well-defined ground state and allows for negative "photon excitations" which are often interpreted as the emission of photons. To identify the physically correct occupation of a "photon-dressed" system in Floquet theory is a very subtle issue of intense discussion in the community [45]. In contrast to Floquet-theory, for the fully coupled matter-photon problem the number of possible photonic excitations k are bounded from below by k = 0, which is the vacuum of the electromagnetic field.

B. Coupled equations in separated Hilbert spaces
Let us next derive, by applying the full Hamiltonian (1) on the discussed expansion, coupled equations for the parametrically dependent subsystem wavefunctions. We employ here that in configuration space (R n , r, q) multiplication operators like potentials and interactions commute with the wavefunctions. This allows us to formally treat the dependence on other subsystems parametrically, for example, for the electronic wavefunction the electronnucleus interaction becomes In a more physical interpretation, as long as the particles have no quantum-character, the coupling between different systems is purely recovered by their meanvalues. This simple factorization is no longer valid if the Hamiltonian includes derivatives and therefore assigns a quantum-character to the particles. The kinetic contributions act on all functions, that is, there are non-vanishing contributions such as ∇ Rj χ µ i ∇ Rj ψ k µ , which explicitly couple nuclear, electronic and photonic degrees of freedom. The initial eigenvalue equation then becomes where the photonic coupling elements are given by (see App. B for a detailed derivation) and In this form we can identify three equations which have to be solved in a self-consistent manner in order to satisfy the above combined nucleus-electron-photon problem. The first one is given by the photonic equation of (11). Since the total dipole R that shows up in the coupling to the photon subsystem wavefunction is given by merely the electronic and nuclear coordinates, for eigenstates with the term (11) simplifies to ∞ µ,l,k=0 The dependence of the photonic eigenvalue on the total dipole therefore leads to a photonic potential-energy surface ε k ({R}) ≡ ε k ({R n , r}). We denote the parametrically dependent photonic Hamiltonian as the photonic Born-Oppenheimer Hamiltonian We can then shift the photonic potential-energy surface into (10) and define a photon-adopted electronic Born-Oppenheimer Hamiltonian according tô With this definition we then solve for the electronic eigenfunctions of (10) including the photonic potential-energy surfaces aŝ The remaining nuclear equation is now the combination of (7), (8), (9) and the additional potential-energy sur- We finally end up with three equations (14), (18) and (19) that have to be solved self-consistently. Their physical interpretation is that electrons move adiabatically on photonic energy-surfaces ε l (r, {R n }) while excitations of the electronic system are coupled by photonic excitations l. The bilinear coupling (12) transfers electronic momentum between different electronic states, mediated by photonic excitations. However, within the long-wavelength approximation, the photon itself does not transfer momentum to the electron, see also App. C. The quadratic coupling (13) in contrast is an energetic shift between eigenstates. The equivalence between the presented Born-Huang expansion and the Power-Zienau-Woolley transformation [46] is elaborated in section III A. In combination, (14) and (18) constitute the polaritonic subsystem which is interacting with the nuclei. The nuclei move adiabatically on the polaritonic surfaces E ν (R n ) with additional couplings mediated by photonic and electronic excitations as well as mixed electron-photon excitations. Bare non-adiabatic coupling elements without photons have to be adjusted, as for example discussed in Sec. IV A, to account for novel contributions and the change of the electronic structure under photonic influence ψ l ν |∇ Rj |ψ l µ . The non-adiabatic couplings are often negligible as long as the Born-Oppenheimer surfaces E ν (R n ) (in our case these are the polaritonic potentialenergy surfaces) are energetically well separated but become relevant close to conical interactions [28,[47][48][49]. Finally, we point out that a similar construction of coupled equations can be deduced also for the time-dependent case (see App. D).
C. Analytic solution of the photonic subspace So far we have reformulated the high-dimensional problem of a correlated nucleus-electron-photon system into a set of three lower-dimensional, but coupled equations. The solution of these coupled equations is of course still as hard (maybe even harder) than the original problem. In order to make the problem numerically tractable we either need to introduce approximations or we provide analytical solutions to parts of these equations. In this subsection we will do the latter and bring the problem into a more compact form. We use that in the photonic equations (14) and (D2), respectively, the parametric dependence on the total dipole allows for an analytic solution. The dipole merely introduces a coherent shift in the photonic harmonic oscillator equations. Let us elaborate on this briefly.

Equation. (14) obeys the generic form
with a given initial state φ 0 (q) = φ(q, 0), q 0 = φ 0 |q|φ 0 andq 0 = φ 0 |p|φ 0 . As R(t) is a given external perturbation in this context, the solution to the above equation is The coherent shift operator, which is a combination of a time-dependent translation and boost (translation in momentum space), fulfills condition (21) up to a timedependent phase [50] where the classical trajectory of the displacement coor-dinateq is given by We therefore see explicitly that a coherently driven photon mode just follows exactly the classical trajectory and only other observables provide access to the "quantumness" of the photon mode. That this holds is due to the quantization procedure of the electromagnetic field, which makes sure that without coupling to the matter subsystem the expectation values of the field operators reproduces the Maxwell equations in vacuum. So mode by mode the quantum harmonic oscillators need to reproduce the classical oscillators as long as we only have external sources [13,51]. In the time-independent case, where the classical equation of motion merely reduce tȯ q = 0 and q(t) = q 0 we immediately arrive at the shifted eigenstates since it has to hold that q 0 = − λ ω · R, where R is the total static dipole. Consequently we have for the quantum stateŝ Since in the photonic subsystem we only have different shifted harmonic oscillators the resulting photonic wavefunction parametrically dependent on R becomes where we have defined q 0 α = − λα ωα · R and we have a multi-index k ≡ (k 0 , ..., k Mp ) that collects the individual mode excitations. With this we find the parametricallydependent photonic energy surfaces as As a consequence the photonic energy-surfaces in Eq. (17) are just constants that merely shift the total energy. Therefore, the photons do not affect the coupled equations directly, since they do not enact a force on the other particles, that is, ∇ R ε k = 0. Instead the photons affect the electrons and nuclei only via the non-adiabatic coupling elements. In contrast to the electron-nuclei coupling elements these photonic coupling elements are known analytically, that is, given in Eq. (12) and (13).
Thus the photons in the long-wavelength approximation do not introduce extra quantities that need to be determined numerically but rather change the usual Born-Huang expansion and lead to new but analytically known couplings.

III. IMPLICATIONS
After we have solved the photonic subsystem analytically in the presented Born-Huang expansion, let us restate the coupled problem that solves Eq. (3) in a more compact form. By solving and subsequently with the obtained polaritonic Born- we find that the exact eigenstate is recovered as . Let us consider in the following a few implications of this new form of the original problem before discussing some approximation strategies.

A. Relations to the Power-Zienau-Woolley transformation and the diabatic picture
Let us first compare the above form with the Power-Zienau-Woolley transformation [46] that allows to define a multipole form of the minimal coupling Hamiltonian. The length-form Hamiltonian that we use here can be obtained by approximating the quantized vector potential operator by its value at zero (or any other reference point)Â(r i ) →Â(0). This approximation leads to the momentum form of the minimal coupling Hamiltonian [13,20,38] (without loss of generality we just show how the momentum of the electrons are adopted) wherep α is as defined before. Then performing an operator-valued boost of the form exp[−iR ·

Mp α=1
λα ωαp α ] leads to the length-form of Eq. (1). In the photonic subsystem in the Born-Huang expansion this operatorvalued boost turns into a translation of the displacement coordinate for a given dipole moment R and therefore is equivalent toD † (q 0 ), which is just the tensor product of the individualD † (q 0 ) as defined in Eq. (22). Thus when we apply the coherent shift operator in the Born-Huang expansion we take the step back to the momentum form of the long-wavelength approximated minimal-coupling Hamiltonian. Formally we can connect the above operator-valued boost to a multipole expansion of the Power-Zienau-Woolley transformation e i drP(r)·A(r) of classical physics, where P(r) is a polarization field [46,52,53]. By then promoting the classical quantities to operators a multipole form of the minimal-coupling Hamiltonian can be defined [52,53]. In the above momentum form of Eq. (26) it becomes straightforward to generalize beyond the long-wavelength approximation by reverting our very initial assumption of spatially independent modes and keeping explicitly λ α (r) [13,38]. Keeping the spatial dependence will change our expansion considerably, as it is no longer only the total dipole the photons couple to. In the general case the photon subsystem wavefunction will depend parametrically on all individual coordinates of the other particles, while in a formal multipole expansion higher order terms will be introduced. This will lead, similar to the decoupling of electrons and nuclei, to non-constant photonic surfaces that depend on the spatial position of the matter subsystem. As a consequence the combined system will try to occupy a minimal potential point given by the involved modes of the photon field.
In the context of the momentum form of the nonrelativistic QED Hamiltonian we also briefly want to highlight the difference between an adiabatic and diabatic picture [54]. In the length form we have decided to use the adiabatic picture, that is, we assume that the nuclei move on electronic surfaces and the electrons move on photonic surfaces. The physical rationale is that we assume that the nuclei "move slower" than the electrons, which in turn "move slower" (be aware that displacement coordinates q are not real-space movements) than the photons. This has been expressed formally in our Born-Huang expansion by taking into account parametrically only the classical quantity R for the photon subsystem wavefunction. But in principle we can use other arrangements. For instance, in Ref. [37] the photons were grouped with the nuclei due to their obvious similarity to the simplest approximations of the quantized nuclei by (harmonic) phonons. The electrons are therefore considered "fast" with respect to the rest. This is a diabatic picture. In this case it is the photons and nuclei that move on the electronic potential-energy surfaces. Take, for example, the polaritonic wavefunctionΨ µ but now we expand in terms of electronic subsystem wavefunctions parametrically dependent on the momentum-form photon coordinates and the position of the nuclei as If we consider the momentum form Eq. (26) we can get rid off the parametric photon-coordinate dependence by a boost in real space of the form exp[−iR ·

Mp α=1
λα ωα p α ]. This boost is the transformation from momentum to length form but with p α as a number not as an operator. By this we have eliminated the dependence of the electronic wavefunction on the photon coordinate. Therefore, for what the photonic subsystem wavefunction Φ k µ (p, {R n }) is concerned, the electronic potential-energy surfaces become flat and it is only coupling elements that will change the photon field. In the diabatic case, however, there are usually no analytic solutions for the electrons available.
Let us at this point also highlight, that the length and momentum form of the non-relativistic QED Hamiltonian in the long-wavelength approximation, get their name from the above discussed translations and boosts. In the momentum form, a translation in real space has no effect on the photons at all, and a boost of the matter subsystem merely shifts the photon field coherently, leaving the eigenstates invariant. In the momentum form it is the variations/fluctuations in momentum space, for example, due to non-zero momentum matrix elements ∇ lk , that lead to physical effect on the eigenstates. In the length form (after we have performed already an operator-valued boost) we can translate the matter subsystem in real space and it only amounts to a simple coherent shift of the photonic subsystem, which again leaves the eigenstates invariant. It is the real-space variations/fluctuations, for example, due to non-zero dipole matrix elements, that lead to physical effects on the eigenstates (see Sec. IV B for an example).

B. The relevance of self-polarization
The physical and mathematical relevance of the selfpolarization, also called the dipole self-energy, α (λ α · R) 2 /2 has been discussed in multiple publications [23,24,38,55,56]. The here presented reformulation of the nonrelativistic QED problem in the long-wavelength approximation allows an alternative, very intuitive and physical explanation of the most important fact, that is, that the coupled nucleus-electron-photon system is not stable without the self-polarization [38]. This becomes evident if one considers the photonic surfaces ε k as defined in Eq. (23). In the case that the self-polarization term is discarded in the original Eq. (20), the energy expressions ε k get a dipole-dependent shift −(λ · R) 2 /2. Therefore the photonic potential-energy surfaces ε l in Eq. (24) would not be constant, instead they would be quadratically divergent. This leads to a linearly divergent force As a consequence Eq. (24) does not posses an equilibrium solution (if we consider infinite space), i.e, E µ ({R n }) → −∞, and with this also the original problem becomes unbounded from below, that is, If we restricted to a finite box then the ratio between coupling strength and box size (which in turn restricts the value of R) determines whether we get a stable system or whether the minimal energy is reached by tearing the system apart (nuclei on one side and electrons on the other to maximize R). With the self-polarization such a problem does not appear and a stable, box-size independent solution is well-defined. That this problem is not more often encountered has to do with the fact that mostly the coupled system is solved in a restricted state space. As long as our electronic expansion is small (see the discussion in Secs. IV A, IV B and V) this divergence does not emerge because the limited set of electronic states cannot describe such a divergent shift. Take, for instance, the minimal setting of only a single electronic excitation. The system is then represented by the Pauli-matrix algebra witĥ R →σ z and (λ α ·R) 2 → λ 2 α1 with the identity on the two dimensional vector space denoted by1. The self-polarization has no effect in this smallest electronic space besides a shift in energy. However, this is not the case for the exact solution and the missing contribution is essential to capture ground-state observables. This effect will be illustrated by a numerical example in Sec. V.
Let us finally remark that the self-polarization contribution is necessary to connect the momentum and the length form of the non-relativistic QED Hamiltonian in the long-wavelength approximation and that the selfpolarization term is not equivalent toÂ 2 (0). Further the influence of the self-polarization term for dynamical phenomena and within the weak-coupling regime is less evident. Only after sufficiently long propagation large deviations will be observable.

C. Observables
Although excitations of the photonic system are just implicitly addressed in the nuclear (25) and electronic (24) components, observables can be calculated as usual For instance, we can consider photonic observables, such as the mode-occupation. The physical interpretation of displacementq α and furthermore creation and destruction operators change between momentum and length form, that is, the occupation of mode α in momentum form is related to the electro-magnetic field energydensity [53]. We solve the photonic system as discussed in Sec. II C. The coherent translationD transfers the wavefunction Ψ i →Ψ i into the momentum frame while the operators remain formally identical. In our current formulation this leads to the following mode-resolved form In this way we can construct similar photonic observables such asq α and any other observable of the combined nucleus-electron-photon system. But of course, we can also access other observables, such as the electronic density or excitation energies (see Sec. V B for examples). As it will turn out, depending on which observable we want to consider, we will need different levels of approximations to the full Born-Huang expansion to get accurate results. Overly reduced models will not be able to capture -even qualitatively -some considered observables.

D. Born-Oppenheimer including photons: the explicit polariton approximation
The coupled equations of the Born-Huang expansion, that is, Eqs. (24) and (25) in its compact form, can be drastically simplified if certain coupling elements can be neglected, that is, for specific subsystem wavefunctions the derivative with respect to their parametrical dependence vanishes. This is equivalent to the fact that the induced energetic surfaces are well separated. That this holds true for all electronic subsystem wavefunctions is one of the basic assumptions of the traditional Born-Oppenheimer approximation of quantum chemistry. If we make the same assumption we still have the photonic coupling elements that connect the different nuclear subsystem wavefunctions. So in order to arrive at a similar simple form as the original Born-Oppenheimer approximation we have to also assume that these coupling elements are negligible. This will be the case if we -besides the usual physical rationale of "slowly moving" (almost classical) nuclei -also assume that the frequencies of the photonic modes is mainly in the electronic energy range. In this case it seems reasonable to assume that the coupling elements are not important, similar to the starting assumption of Ref. [28]. Under those assumptions, electron-photon and nuclear coordinates factorize, leading to The many-body wavefunction Ψ i is now approximated by a polaritonic excitation ν and uncorrelated nu- As a consequence of the Born-Oppenheimer approximation, the nuclear excitations are calculated on polaritonic quasiparticle energy-surfaces. The photonic contribution depends parametrically on the total dipole R. This leads to the following simplified form of the Born-Huang expansion where the electron-photon-solution feeds into the nuclear equation Clearly, other approximations of the fully coupled Born-Huang expansions are possible, for example, we could assume that the photonic modes are only within the nuclear-excitation energy range. To distinguish the above specific choice we call the coupled Eqs. (27) and (28) the explicit polariton approximation to highlight that we consider nuclei moving on polaritonic surfaces.
To solve the explicit polariton equations we restrict first our photonic excitation space to a maximum number of possible excitations l ∈ {0, 1, ..., l max }, that is, the sum ∞ k=0 → lmax k=0 is truncated at a finite photonic occupation. The cost we then pay is that we effectively end up with the same high-dimensional problem of the initial fully coupled electron-photon Hamiltonian in a restricted photon-space. For the electronic subsystem we will in the following, however, also use a suitable basis expansion which we will truncate in practice. How this is done and how the polaritonic problem of Eq. (27) is expressed in such a basis expansion we will discuss in Sec. IV A. The explicit polariton in a finite basis expansion (for photonic and electronic subsystem, respectively) can be efficiently used to solve coupled systems for weak and strong coupling for finite and also periodic systems (see App. C). The additional effort to include the photonic interaction into a common solution of the nuclear-electron problem is computationally negligible as long as we can treat the relevant set of the excited electronic states efficiently and the accuracy of the restricted electronic and photonic basis can be estimated based on physical intuition or arguments. That this is not always as intuitive as one might hope for and that there is a subtle interplay between number of electronic and photonic excitations that we take into account we will show explicitly in Sec. V. Further, while certain simplified quantities such as changes in excitation energies due to the emergence of polaritonic states might seem already converged only after taking into account a minimal number of photonic excitations and electronic basis states, the convergence in other quantities and especially of the underlying approximate wavefunction with respect to the photonic and electronic basis set might be much slower. We will investigate the restricted basis issue later and will as reference use the frequently employed restriction to merely just one photonic excitation per mode. We will name the resulting approximate solution the small polariton in analogy to Ref. [57], where the diabatic electron-nuclear problem is approached in a single-hopping approximation. Therefore the small polariton involves at maximum a single excitation per mode l α ∈ {0, 1}. In Sec. IV B we will use this approximation to connect the above quantum-chemical approach to quantum-optical models that are used to describe coupled matter-photon systems in a simplified way and we will discuss implications.
While the above explicit polariton approximation is expected to be a reliable approximation for the weak and strong coupling regime within a restricted set of excitations, in the domain of ultra-strong coupling, that is, λ α / √ 2ω α ∼ 1, other possible Born-Huang expansions, for example, such as the diabatic approach briefly discussed in Sec. III A, might be more appropriate. Let us also briefly remind the reader at this point about the connection of the explicit polariton approximation and the Floquet approach. As we have seen in equilibrium, the effect of the photons on the electronic structure is purely determined by fluctuations within the shifted harmonic oscillators which connect electronic excitations. The classical contribution, namely the shiftq α →q α − q 0 α , is without effect on the electronic structure. This clarifies that the eigenstates of Eq. (27) are polaritonic wavefunctions projected on the electronic subspace, that is, the solution of Eq. (27) can be constructed as a vectorial expansion for each mode in sub-blocks according to ψ ν = {ψ 0 ν , ψ 1 ν , ..., ψ lmax ν }. Diagonalizing this matrix is structurally very similar to solving the Floquet problem with the difference, that the energy is bounded from below and the second order contributions participating via ∆ lk involve additional couplings. Expanding the full basis explicitly corresponds to the photonic Fock-state expansion in momentum form.

IV. CONNECTION TO QUANTUM-OPTICAL MODELS
We have now approximated the full Born-Huang expansion by the explicit polariton approximation given in Eqs. (27) and (28). This has already reduced the complexity considerably. Still solving the explicit polariton approximation in full real space for the electrons is very challenging due to the involved Coulomb interaction among the particles. Thus for practical purposes also a restriction of the electronic basis set is desirable. Also in most models for coupled matter-photon systems the solution of the full many-electron problem is avoided by working in a restricted basis for the electrons. In its simplest form this restriction leads to a two-level approximation to the matter subsystem, which is a common approximation strategy encountered in quantum optics, for example, in the Jaynes-Cummings model or for many two-level systems in the Tavis-Cummings model. In order to connect to these well-known models and highlight possible problems in an overly simplified treatment as well as for practical purposes we will in the following also restrict the electronic basis set.

A. Basis-expansion in electronic eigenstates
To restrict also the electronic subspace we expand the polariton in a common electronic basis. Although this can be any complete basis, the exact electronic eigenfunctions will simplify the following steps. What will emerge is a way how to use the basis in such a way, that the photonic influence can be treated by a simple diagonalization in this basis. As a (not necessarily orthonormal) basis we choose some electronic wavefunctions ϕ n (r, {R n }) and then we expand the many-body electronic wavefunction This expansion applied to Eq. (27) leads to where S jn = ϕ j |ϕ n e H l BO,jn = ϕ j |Ĥ l BO |ϕ n e c ν ln S nn c ν n l .
Here we suppressed the parametric dependence on the nuclear configuration for brevity and we did not assume that the basis is orthonormal. We remind the reader that the analytically known coupling elements ∇ kl and ∆ kl are completely independent of nuclear and electronic coordinates.
Along this line it becomes evident how the non-adiabatic nuclear-polariton coupling elements ψ l ν |∇ Rj |ψ k µ e can be calculated from the knowledge of bare coupling elements ϕ n |∇ Rj |ϕ m e . We expand into (30) and observe that two contributions occur. The second line represents the super-position of bare non-adiabatic elements according to the participation of those states in the polariton solution. Furthermore, the linear coefficients themselves introduce an additional component ∇ Rj c µ mk ({R n }) as their value depends for example on the dipole-transition element ϕ n (r, {R n })|r j |ϕ m (r, {R n }) e which will change depending on the nuclear configuration as discussed in Sec. IV B.
Next we choose a basis and here we draw the connection to quantum chemistry by a selection of Slater determinants. For this, we multiply from the left with an electronic (excited) Slater determinant The opera-torT (n) excites the wavefunction n-times, that is, ϕ i (r, {R n })|T (j) ϕ 0 (r, {R n }) e = S ij . Further, within Hartree-Fock theory, the photonic coupling elements would open the usually closed Roothaan-Hall equation [58] into all possible excited components. Here, however, we do not restrict to Hartree-Fock theory but rather perform a configuration-interaction expansion of Eq. (29) within a restricted subspace of photonic and electronic excitations or with a variational minimization of the polaritonic energy according to In this form it becomes apparent that we can in principle use and extend common quantum chemical minimization procedures, such as Configuration-Interaction or Coupled-Cluster techniques [59], to solve the polaritonic problem.

B. A natural connection to model-systems
Let us now build the connection to a very common system applied in quantum optics. We assume a set of identical, independent molecules, the excited Slater-determinants are then exact solutions of the electronic system, coupled by one photonic mode α = 1. With close to degenerate energies, the first excited component of this small polariton subspace is represented by the unperturbed set n=1 ; E ν=1 = E 0 + ∆E e } with bare electronic excitation ∆E e and detuning of the photon mode δ is an anti-symmetric (Ŝ − ) many-body fermionic wavefunction with a symmetrized (Ŝ + ) excitation represented by the single excitation operatorT (1) . In addition to this fully symmtrized excitation, there exist N e − 1 excitations with mixed symmetry which are often referred to as dark states. They posses a vanishing transition dipole due to anti-symmetric components in the many-body excitation [28] and do not directly couple to the transversal mode by emission or absorption of a single photon. Nevertheless, higher order processes as the Lamb-shift still affect those states. As is usually done, we will disregard these dark states in the following. Next, instead of the electronic momentum elements ∇ jn e , in molecular systems the dipole moment is typically preferred in calculations . It is connected to the momentum element ∇ jn e by the operator identity In the approximation of electronically non-interacting molecules, that isŴ ee (r, r ) = 0, the electronic eigenvalues are isolated molecular excitations ∇ jn e = − j|[Ĥ,r] − |n e = (E n − E j ) j|r|n e = ∆E nj r jn and the energy-difference leads to a transformation of the characteristic bilinear coupling scale The coupling elements in length and momentum form are identical up to a factor (ω − δ)/ω [53]. The positionr is defined in relation to the center-of-charge. We now solve this ansatz in the smallest possible subspace for electronic j and photonic excitations l such that (j, l) ∈ {(0, 0), (0, 1), (1, 0), (1, 1)} were we choose the pure electronic Slater-determinants ϕ M B n to form an orthonormal basis S jn = δ jn .
The two-level small polariton approximation becomes extensively problematic as coupling or detuning increases and is not suited to describe the ground-state of a correlated system although it includes the anti-rotating contributions, that is, there is no rotating frame assumed. Especially the self-polarization contribution can become troublesome as discussed in Sec. III B. This is made more explicit and elaborated in Sec. V. The small polariton with a single electronic excitation can be represented by the 4-by-4 matrix which can be block diagonalized in a degenerate single excitation 2-by-2 matrix subspace whose solution resembles the upper and lower polariton with an energy-or Rabi-splitting of ∆E = δ2 +Ω 2 , the generalized collective detuning δ = δ + L and the collective Rabi-frequencỹ Hereby r {Rn} 01 depends parametrically on the nuclear configuration in terms of one individual molecular state but identical for all molecules. Assuming the molecular system consists of different species or molecular configurations, one can partially symmetrize those to recover a collective behavior for each set of identical configurations [28,60]. The obtained energy splitting fits the lowest order contribution, that is n = 0, of the Tavis-Cummings model [27] up to the factor (ω − δ)/ω with It becomes evident that within such a strongly simplified expansion, that is a single photonic and electronic excitation, the resulting solution will vary depending on the selected coupling form. The only exception occurs exactly on resonance or for a completely converged expansion (see Sec. V). The conserving fluctuations ∆ ll introduce a collective frequency or diamagnetic shift that detunes electronic and photonic excitations and can be interpreted to some extend as a collective renormalization of the electronic excitation or vice versa the photonic energy. In general, the coefficients c i=1,±

01
(δ({R n }),Ω({R n })) depend on the number of particles, not only via the Rabi-splitting but also due to the collective detuning. Furthermore, their behavior is essential close to conical intersections as they determine the non-adiabatic nuclear-polariton coupling elements ψ l ν |∇ Rj |ψ k µ e discussed in Sec. IV A (see Eq. (30)). Such a detuning as shown in Eq. (32) caused by the collectivity of a frequency-shift was measured [61,62] and extensively discussed in Ref. [63]. As a consequence of the collective detuning, the resonance does no longer coincide with zero detuning δ but appears for values δ = ω − ∆E e = −L. As a consequence of the quadratic coupling, the collective molecular participation or the coupling itself has to be significant otherwise the effect is negligible. However, within strong coupling, as for example ω α ∼ 0.1, √ N e λ ∼ √ 100 · 0.01 ∼ 0.1 we can see that L is on the order ofΩ and consequently non-negligible. In combination with the remaining block involving the counter-rotating contributions the collective many-body (cavity-matter) energetic levels are given by Here, E − 1 and E + 1 correspond to lower and upper polariton with conserving excitations, that is, the states are super-positions of excited matter or cavity but never both. In contrast, E 0 describes the renormalized groundstate of the correlated system while E 2 is connected to the coupled state where electronic and photonic system are excited simultaneously. If we drop the excitation non-conserving couplings, that is, the coupling from no excitation at all to both subsystems are excited ∇ 01 e · ∇ 10 and ∇ 10 e · ∇ 01 , E 0 and E 2 do not change beside energetic shifts via L. The small polariton does naturally include counter-rotating terms within this single excitation subspace as presented above. Whether the small polariton approximation is reasonable depends on the relevance of multi-photon processes. As long as we expect them to be negligible the small polariton approximation is expected to capture the essential physics. This gives a very intuitive criterion to judge the quality of this minimal model. But, as will be discussed in Sec. V, this intuition can sometimes be a little misleading as the quantitative agreement depends on the interplay of photonic and electronic excitations as well as the observables of choice.
Next we consider what the small polariton can teach us about the alleged Rabi-phase transitions in coupled matter-photon systems. From Eqs. (33) and (34), we can observe That is, what was the ground state at λ = 0 is connected to the ground state for all λ > 0. Only in the thermodynamic limit of N e → ∞, |Ω| → ∞, the groundstate becomes degenerate with the lower polariton. Performing the rotating-wave approximation and solving the following Tavis-Cummings model, namely shifting the first polariton by Eq. (31) will result in a crossing of lower polariton and the non-sensitive ground-state, that is, E − 1 − E 0 < 0 ∀λ > λ c with a critical coupling λ c . This type of a Rabi phase transition is not observed within the small polariton. In contrast to the out-ofequilibrium super-radiance transition which has been experimentally verified, a phase-transition in equilibrium is still highly debated [33,[64][65][66]. Also note that neglecting the counter-rotating contributions is only valid forΩ ∆ . Although the single-photon collectivity that is captured by the small polariton approximation is dominant, it only corresponds to the linear contribution for collective effects. Nonlinearities can be introduced by matterphoton interaction where also the self-polarization can be substantial and/or multi-photon participation which can lead indeed to drastic transitions as will be presented in Sec. V and has been observed also in Refs. [24,40]. Multi-photon participations are not represented in the above single-photon approach and a non-trivial collectivity leading to a ground-state transition can consequently not be ruled out from the above simplifications. As for example presented in Refs. [33,65,66], the phasetransition does drastically depend on the interplay between self-polarization, bilinear coupling and Coulomb interaction. An extended discussion of this question including an accurate investigation of the electron-photon coupled system by high-level renormalization-group techniques is currently in progress [67].

V. NUMERICAL DETAILS
Finally, after we have simplified the original nucleuselectron-photon problem by the explicit and then the small polariton approximation, we want to investigate the reliability of these approximations. We will do so by considering a real-space system that we then approximate by a different number of basis functions for the photonic and electronic subspaces. To keep the numerical calculation tractable we will focus on the polaritonic subsystem and freeze the nuclear coordinates. By assumption, it is only the electronic subsystem that is directly affected by the coupling to the photons and the nuclei merely feel the presence of the photonic modes as a change in the potential-energy surfaces. We therefore investigate the accuracy of the different levels of approximation in terms of the polaritonic observables and wavefunctions only. Before we do so, we want to give a brief recapitulation of the general explicit polariton approximation and its numerical implementation. This work-flow can be straightforwardly extended beyond the explicit polariton to the full Born-Huang expansion of Eq. (24) and (25).

A. Computational scheme for the explicit polariton
The derived equations incorporate the photonic degress of freedom for equilibrium QED calculations in a consis-tent way into the common Born-Oppenheimer quantum chemistry approximation. The additional computational effort remains limited as the photonic subspace can be represented as bare excitations of electronic ϕ n (r, {R n }) and nuclear χ µ ν (R n ) wavefunctions. Within, for example, configuration-interaction or coupled-cluster techniques [59], the excitations are part of the minimization procedure of the fermionic equation anyway and can be used directly within the minimization procedure as schematically presented in Fig. 1. The computational work flow of the explicit polariton approximation then becomes: In the first instance 1), we have to clarify how many modes with corresponding frequency ω α are essential to describe the field-matter interaction within the cavity. Commonly a single mode is assumed but there are situations where more are relevant [10]. This can be approached by comparison of cavity and matter spectral function. It determines which modes and frequencies are essential as well as the nature of matter excitations, that is, dominantly nuclear or electronic excitations. In instance 2), we use the insight from step 1) to calculate or select the necessary set of electronic states for a sufficient domain of parametric nuclear positions {R n }. The selected set of electronic states does not have to be a set of eigenstates of the electronic problem but it will render the following step simpler. Next, in step 3), those bare electronic states are mixed by Eqs (18) or (29) into coupled electron-photon states (polaritons). As discussed before, due to the structural similarity to Floquet theory, small adjustments in existing implementations can lead to efficient solutions of this step with marginal additional effort. Finally, we solve the nuclear component (19) in step 4) on the polaritonic energy-surfaces, potentially including non-adibatic coupling elements from electronic, photonic or mixed excitations. This procedure has to be iterated until self-consistency is reached. Higher excitations are typically stronger delocalized which is more likely to appear with photonic interaction. The ensemble of parametric nuclear values {R n } in the polaritonic subspace has to be potentially adjusted. Beyond the long-wavelength approximation, this would also take place in the photonic subspace. Alternatively for step 3), the Coupled-Cluster technique or any other modern quantum chemistry scheme to deal with electron-correlations could be directly applied to a polaritonic Slater-determinant in Born Oppenheimer approximationΨ 0 ({R n }, r, q) ≈ ψ 0 0 (r, {R n })Φ 0 (q, {R}) with a possible cluster excitation operator for electronic and photonic system  Here Hp is the photonic subspace, Hs is the electronic subspace and Hn is the nuclear subspace.
Step 1) is a purely conceptional selection of energetic domains, 2) corresponds to the very common process of solving the electronic many-body problem with parametric nuclear dependence while 3) represents the essential and computationally novel step. Solving the polaritonic step 3) can for example be efficiently done with slight adjustments of available Floquet implementations.
Step 4) corresponds to the familiar nuclear solution with potentially novel and adjusted non-adiabatic couplings as discussed.
Recently, together with others, we presented how an efficient implementation of the QED ground-state correction can be accomplished by means of quantumelectrodynamical density-functional theory within the optimized effective potential (OEP) approach [24]. There, the Kohn-Sham system is factorizing photonic and electronic system, which accounts for the classical shift of the Harmonic oscillator basis. The remaining quantum nature of the photonic orbitals is contributed by an effective potential which can be derived from perturbative corrections Φ (1) iσ,α to the Kohn-Sham orbitals φ iσ . Those corrections carry a single photon and are the solution to a Sternheimer response equation, that is, they incorporate the effect of the response of the system to photonic fluctuations. We can identify a structural similarity of polaritonic associated states ψ l ν (r, {R n }) and perturbative orbital corrections. This is especially prominent in the mode occupation of section III C.
Since the accuracy of the nuclear Born-Oppenheimer approximation is well known, we focus here on the polaritonic contribution 3) to gather additional physical insight for the presented approach, that is, we skip the final step in Fig. 1. Additionally, 3) is the step that connects firstprinciples to quantum optical model systems and we will see how a reduced set of electronic states affects the qualitative and quantitative behavior. An efficient implementation of the explicit polariton can start from Eq. (29). To allow for an exact reference, that is, distill the effect of photonic interaction, we solve the electronic system by exact diagonalization and use the unperturbed states in Eq. (29). The coefficientmatrices c ν nl relate a set of bare electronic states ϕ n (r) with n ∈ {0, 1, ..., n max } to mode excitations l ∈ {0, 1, ..., l max }. The full space resolved polaritonic eigenfunctions ψ ν (r) = {ψ 0 ν (r), ψ 1 ν (r), ψ 2 ν (r), ..., ψ lmax ν (r)} are represented in the limited set of bare electronic eigenfunctions ϕ n (r). The resulting eigenstates represent the electronic as well as photonic contribution for each polariton eigenstate ν. For example, the ground-state density of the correlated electron-photon system is given by n λ>0 ν=0 (r) = nmax n,n =0 dr 2 dr 3 ...dr Ne ϕ * n (r, {R n }) The density is an important observable in quantum chemistry that allows for intuitive interpretations and gives information, for instance, on the nature of chemical bonding [68]. The change of the real-space resolved density in the correlated matter-photon system can be helpful in interpreting the effect of the emergence of polaritonic states. Such changes do also directly affect the non-adiabatic coupling elements and can prove useful to get a more detailed understanding of the influence of cavity photons on specific chemical processes, as for example transition states, solvation energies, and involved processes within energy transfer. A further observable of quantum-chemical interest is the energydifference between ground and first excited state ∆E 0→1 as those energies present the surfaces on which the nuclear wavepackets move dominantly. On the other hand, of quantum optical interest is for example the energydifference between first and third excited correlated state ∆E 1→3 which represents the Rabi-splitting in our numerical example for weak coupling of Sec. V C. The cavity mode occupation â †â as presented in Sec. III C is an additional observable of interest. The selected designation of weak or strong coupling are based on the question to which extend the electronic structure is distorted by the transversal photon interaction. As we consider a perfect cavity without loss, one might assign both to the strong coupling regime as dissipative effects are by construction not existing. It is important to note that except of the density all the shown observables are integrated quantities. It is usually much more involved to accurately capture non-integrated (space-resolved) quantities than integrated quantities. This will become clear in the numerical example discussed in Sec. V C. So the accuracy of the approximation strongly depends on the quantity that is considered and simple models, while reliable for certain quantities, are less so for other observables.
For computational convenience, we map the coefficient matrix on a vector c ν nl → c ν m , m = n · l max + l to recast the explicit polariton into a diagonalizable matrix equation.

C. 2D GaAs quantum ring
Now we will investigate in detail the different possible levels of approximations. By this we mean that we consider the numerical convergence with respect to the electronic n ∈ {0, 1, ..., n max } and the photonic subspace l ∈ {0, 1, ..., l max }. The computational demand in a truncated subspace is significantly lower than an exact diagonalization in the full electron-photon space while reaching a high level of accuracy up to strong coupling. This is a consequence of the fact, that the electronic eigenstates are a well suited basis-set for equation (29). It changes if the electron-photon correlation extends into (ultra-)strong coupling where the electronic structure is drastically distorted by the correlation. The limited sub-set n λ 1 max is then no longer sufficient to describe those distortions and we have to extend our minimal set n λ 1 max < n λ∼1 max to describe the correlated system. This will clarify, that a few-level approximation is problematic by construction in the (ultra-)strong coupling domain. We investigate a model of a GaAs quantum ring introduced in previous works [24,37,40]. Hereby, one electron is trapped in a two-dimensional 'Mexican hat' potentialĤ en = ξ 1 2 (r ·r) + ξ 2 · exp − (r ·r) ξ 2 with ξ 1 = 0.7827, ξ 2 = 17.70 and ξ 3 = 0.997. We couple this system to a single mode in resonance to the first excitation ω α = E λ=0 with polarization in the diagonal direction λ = λ(e x + e y ). We then vary the effective coupling strength. The electronic structure is solved on a 201-by-201 grid with spacing ∆x = 0.1 and derivatives approximated by forth-order finite-differences. The photonic excitation space is incorporated by 40 Fock numberstates for weak and 80 for strong coupling. For additional information, we refer the reader to Ref. [24].

Weak coupling solution
The exact solution has two limiting cases with respect to the effective coupling. For weak coupling λ = 0.005, the system remains dominantly in its initial configuration and avoids the self-polarization potential, therefore the density is elongated perpendicular to the polarization. This can be illustratively seen in the exact solution Fig. 2 a). While the linear component of the photonic interactionqλ ·R favors density-accumulation along the polarization, the quadratic self-polarization potential (λ ·R) 2 diminishes it. For the ground-state, where the bilinear contribution becomes relevant in second order perturbation theory, the correlated wavefunction is calculated in a potential with two competing contributions. In this weak-coupling example, the quadratic part is dominating and is essential to describe the qualitative behavior of the correlated ground-state. Therefore ignoring the self-polarization would lead to a qualitatively wrong result. The photonic contribution is here a weak perturbation of the system and can be very accurately recovered by effective single excitation approximations, for example the exact-exchange optimized-effective potential approach in the context of quantum-electrodynamical density-functional theory presented in Ref. [24].
The polariton equation effectively mixes momenta of different electronic eigenstates. As a consequence of the radial symmetry, there are two degenerate first excited solutions, one with dominant momentum parallel and one with dominant momentum perpendicular to the field polarization. The minimal basis set for electronic excitations has to include both states, otherwise the dominant contribution is not properly captured. If we take this into account, therefore extend the electronic excitations to two or more, already the small polariton l max = 1 does recover the exact energetic structure for weak coupling quite accurately as presented in Tab. V C 1.
Considering the first three rows in Tab. V C 1 which represent the small polariton approximation, we observe that increasing the number of electronic states does not drastically improve the accuracy at first glance. While the chemically important ground-excited energy- difference does slightly improve, the polariton split and the mode occupation do not or even become worse. If we allow for at most double photonic occupation, this picture changes. Increasing the electronic set does then indeed improve the energies. We can intuitively understand this from the following argument. As long as we allow just for a single mode excitation, the probability of reaching higher excited states is vanishingly small. Including them nevertheless can then lead to inaccuracies since their effect is not properly captured in the approximation. The moment we allow our theoretical description to have two or more mode excitations, we are able to reach higher excited states and those states can now participate by contributing energy and momentum to the correlated solution.
Especially interesting is that the mode occupation introduced in Sec. III C accurately reproduces the one in length form (ex-dE) as long as the the set of eigenstates include just the first set (1S, 1P, ...). As we incorporate the next order (2S, 2P, ...) n max > 8, the mode occupation resembles the one in momentum form (ex-pA). While other integrated quantities change marginally, the mode occupation does not only drastically depend on the number of excited states, it even qualitatively changes its physical interpretation depending on the selected electronic set.
If we compare in Fig. 2 the electronic density influence of the smallest possible set, that is, an effective two-level system with n max = 2 shown in b), to the exact solution a), we find a qualitatively different behavior. The small polariton approximation with n max = 2 cannot even capture qualitatively the right physical behavior even in the weak-coupling limit. Although the small polariton does only couple by bilinear components ∇ lk ∼p, it partially incorporates the self-polarization by construction. As discussed before in Sec. III B, an effective two-level system is not able to properly capture the self-polarization and the same holds true in momentum form. However, if we stay within the small polariton l max = 1 but extend the electronic set to ground-state plus four excited levels n max = 4, the solution presented in c) does switch into the correct orientation. This puts in evidence once more that self-polarization and A 2 part cannot be interchanged or regarded as similar, especially in a truncated electronic set. By further increasing the number of mode and electronic excitations, we start to quantitatively recover the exact result as can be seen in Fig. 2 d). As a result, it becomes evident that a strongly simplified electronic structure is not sufficient to capture the physical behavior of the correlated ground-state. That the electronic set has to be extended quite a bit is also caused by the distorted elliptic shape of the exact solution which is demanding to recover from a superposition of angular momentum states. Additionally, the set of states should be chosen consistently with respect to the distribution of momenta. A selection which is slightly larger but not balanced is in our experience not beneficial.

Strong coupling transition
For ultra-strong coupling λ = 0.4 and therefore λ/ √ 2ω = 0.8, the photonic contribution is no longer a weak perturbation but does drastically reshape the electronic landscape. This is very illustrative in the effect on the electronic density as presented in the exact solution Fig. 3 a). Here, ∆n(r) does no longer avoid the direction of polarization but is distorted along the polarization axis. The bilinear contribution is therefore dominating over the self-polarization and we enter a regime which we could assign to a super/sub-radiant phase as briefly discussed in Sec. IV B. Although the explicit polariton in its perturbative character becomes problematic, including a sufficient number of excited states, which is still negligible in relation to the full electronic Hilbert-space of 40401 eigenstates, can capture the ultra-strong coupling twist. Nevertheless, in this limit it becomes extensively cumbersome to recover the correct energies and densities, therefore the small polariton is no longer a satisfying approximation even for the integrated quantities. The effective two-level approximation b) is essentially scaling the weak coupling solution, deviating strongly in the integrated observ-ables. The energies and mode-occupations in Tab. V C 1 deviate drastically from the exact solution. Recovering the correct energies is now demanding even for larger electronic sets and convergence is slow. Especially, the qualitative behavior of the density demands the next manifold of excited states, that is, 2S, 2P and so on, therefore n max 8. Furthermore, subset c) and d) in Fig. 3 elucidate that it can be even counter-beneficial for the density while beneficial for integrated quantities to incorporate a large amount of photonic excitations without increasing the electronic set. From Tab. V C 1 and Fig. 3, we conclude that within the ultra-strong coupling domain, a fully consistent method is demanded that is able to incorporate a large part of the electronic and photonic Hilbert space.

VI. SUMMARY AND CONCLUSION
We have presented a consistent and comprehensive derivation of approximation strategies for coupled nucleus-electron-photon systems based on the Born-Huang expansion.
We highlighted the connection between length and momentum form with the help of a shifted harmonic oscillator basis for the photonic subsystem, provided an alternative perspective on the relevance of the self-polarization contribution, discussed connections to Floquet theory and deduced an approach how to solve the fully correlated nucleus-electron-photon many-body system for finite and periodic systems. Especially structural similarities to Floquet theory imply the possibility of efficient implementations with marginal effort and furthermore ways to circumvent common problems of Floquet theory. The resulting generalization of the Born-Oppenheimer approximation to including the photons we have named the explicit polariton approximation. Since we could solve the photonic subsystem analytically the numerical costs of this approximation is comparable to common Born-Oppenheimer calculations. If we restrict the number of photonic excitations to one, the resulting simplified equations we then called the small polariton, which we were able to solve analytically when also the electronic excitations where restricted. This allowed us to connect to well-known analytical models of quantum optics and discuss implications, for example, super-radiant phase transitions. Finally we investigated the reliability of the explicit polariton approximation for weak and strong coupling and for different number of basis functions. We found that for integrated quantities such as excitation energies in the weak-coupling regime already the simple quantum-optical models were very accurate but mode occupation and spatially resolved quantities such as the density needed more basis functions to be at least qualitatively correct. On the other hand, increasing the basis functions for the electronic subsystem while keeping the photonic basis functions limited to one excitation led Figure 2. Weak coupling. Photonic influence on the ground-state density ∆n(r) = n λ>0 (r) − n λ=0 (r). Exact solution a) with a 40 photon Fock-space. The two-level approximation is not capable to qualitatively reproduce the correct physical behavior, even in weak coupling. in turn to worse results in certain integrated quantities. A balanced extension of basis functions in electronic and photonic subspace is necessary. In the ultra-strong coupling regime the small polariton approximation, also for the integrated quantities, was qualitatively wrong and many more basis states need to be included.
The presented results clearly highlight how first-principle approaches based on non-relativstic QED in the longwavelength limit and models of quantum optics and polaritonic chemistry are connected. They show how simple models can be extended to converge to the first-principle results by including in a consistent manner more and more basis functions going beyond few-level approximations. To go beyond simplified models or at least to be aware of the limitations of these models becomes increasingly important in the context of polaritonic chemistry and material sciences [8,24,28,49,69]. A theoretical description of strong matter-photon coupling to influence material properties clearly needs a description of the matter subsystem that can genuinely represent the physical and chemical processes and a simplified few-level description might not be enough to capture the complex processes that can happen in real systems. The presented explicit polariton approximation together with similar generalizations of quantum-chemical methods [28,37] as well as generalizations of many-body methods such as density-functional theory for coupled matter-photon systems [13,21,24] paves the way for such studies from first principles [70]. tems as discussed in Sec. IV B. Here the rotating-wave approximation allows to calculate the Rabi splitting (which is used as a measure of how strong the coupling is) by where, as discussed the cavity quantities (the effective volume V and the frequency ω), the dipole-transition element between the ground |g and the excited state |e of the individual two-level systems, the number of two-level systems N and the number of involved photons n ph come into play. If we would like to express the resulting hybridization in terms of an effective coupling of a single two-level system to the mode, we therefore see that we have two collective knobs to turn. The effective coupling can be increased, on the one hand, in this simplified picture by increasing the number of identical systems N inside the cavity, such that we have an increase of √ N . This perspective suggests that the individual matter systems (here simplified to a two-level system) are only affected by the fundamental coupling and consequently do almost not change, and the splitting is a macroscopic quantity of the resulting collective (bright) state of symmetrized super-positions. This behavior is recovered in Sec. IV B, where the unperturbed photonic system (n ph = 0) couples to many identical matter systems.
On the other hand, we also see that we can increase the effective coupling by the number of involved photons. In the simplified Tavis-Cummings description, where the ground-state of the combined light-matter system is merely the tensor product of the bare electronic and photonic ground state, this is associated with higher-excited states of the combined light-matter system. However, as also seen in Sec. IV B beyond the rotating-wave approximation, even within a simplified model there is an increase of photon occupation in the ground state, which in turn can lead to an increase of the effective coupling strength. Physically this increase can be interpreted as the back-reaction of the perturbed photonic vacuum on the matter subsystem, that is, the Lamb shift. The polarization of the electromagnetic vacuum is stronger the more charges are placed within it, for example, in the single-excitation subspace the corresponding energy shift of the ground state goes as N λ 2 (see Sec. IV B). The effective increase of the coupling strength due to the vacuum polarization therefore can change the matter subsystem locally and really affect the electronic structure and the applied local few-level approximation. It is this type of increase of the effective coupling strength we consider if we scale the fundamental coupling for an individual matter system.
Let us give a slightly more formal explanation of those two effects. If we consider the Green's function (propagator) of a single electron or hole G(r 1 t 1 , r 2 t 2 ) = −i TΨ(r 1 t 1 )Ψ † (r 2 t 2 ) (see Refs. [16,17] for more de-tails), then it is influenced by the electronic environment via the Coulomb interaction (longitudinal photons) and the the exchange of transversal photons. The transversal photons themselves are described by their meanfield, that is, classical, contribution and the photonpropagator D α (t 1 , t 2 ) = −i T ∆q α (t 1 )∆q α (t 2 ) , where ∆q α (t 2 ) =q α (t 2 )−q α (t 2 ) is the deviation from the meanfield contribution. In turn, also the photon propagator is influenced by the environment, which results in a set of coupled Dyson equations of the form where we use the abbreviation 1 ≡ (t 1 , r 1 ) and correspondingly for the integrals. Furthermore Σ[G, D] is the matter self-energy and Π[G, D] is the polarization of the mode and G 0 and D 0 are the bare, non-interacting Green's functions of the matter and the photons, respectively. In lowest order in the matter subsystem, we can take the bare photon propagator and get a behavior that resembles the first effect of the collective interaction that goes as √ N . On the other hand, the matter polarization due to the bare matter propagator will change the photonic subsystem and already in lowest order we find that the photon number scales as √ N as well, which will lead to the second collective effect. This second effect will be further influenced by the longitudinal Coulomb interaction and if we go higher in perturbation theory we will get possible non-linear enhancements as well. If we look at microscopic systems where the Coulomb interaction can be substantial, this might enhance this local effect on the matter system even further [66].
This procedure holds true for the nuclear coupling with the adjustment The same is true for the Laplacien acting on the photonic component leading with Ne i=1 = N e to derivatives with respect to the electronic dipole. The element can now be directly calculated with the explicit form of the HO eigenstates and q 0 α = − 1 ωα λ α · R. We could now either directly use this explicit form of the HO eigenstates or the fact, that we explicitly know the translation-operatorD(q 0 α ) (22) such that ∇ rj φ α,k (q α − q 0 α ) = ∇ rjDα (q 0 α ({R}))φ α,k (q α ) = ∇ R exp −iq 0 α ({R})p α φ α,k (q α ) = i λ α ω αp α φ α,k (q α − q 0 α ) and vice versa for Withp α = +i ωα 2 (â † α −â α ) and l|â † −â|k α = k α + 1δ α l,k+1 − k α δ α l,k−1 , we reach the element presented in equation (12). This involves momentum-transfer between electronic states. The same procedure for the second order term leads to (13) which in contrast does not transfer momentum and represents within the same mode fluctuations and squeezings ∆ lk ∼ −p 2 ∼ (â−â † ) 2 ∼ −(2â †â +1)+â †â † +ââ and couples different modes. The second oder coupling term, which corresponds to the diamagneticÂ 2 contribution, can be intuitively understood as renormalization of the electronic or photonic excitation. For instance the collective energetic shift L in Sec. IV increases with increasing photon-number and introduces an additional detuning between excitations. This non-linearity is intuitively elucidating the discussion of App. A.
Appendix C: Implications for periodic systems The initial Hamiltonian (1) is not periodic for λ > 0 in the electronic positions r since the photonic interaction breaks the translational invariance in this direction. While this problem is non-trivial to tackle directly for (1), it is trivial to extend the polaritonic equation (27) to periodic boundary conditions. The additional momentum ∇ r in the linear coupling is well suited for periodicity and the electronic surfaces E ν (R n ) just have to be reinterpreted as polariton-Bloch quasi-particles based on the pure electronic eigenfunctions ϕ k (r, {R n }). Here the momentum k labels a continuous excitation from which the electronic bands in a periodic material emerge via back-folding into the first Brillouin-zone. We present a simple example to clarify the above statement. Assume we want to solve the periodic Kohn-Sham equations for a solid within a cavity and saturated spins. We start with a Fourier-ansatz for a single-particle orbital in the excited state k where G is the inverse lattice vector parametrically depending on the nuclear configuration and n max xyz the corresponding cut-off. By construction of Kohn-Sham theory [71,72], the single particle orbitals are solved in a nonlinear self-consistent procedure using an effective Kohn-Sham potential that mimics the many-body electronic interactions. The set of ϕ k , ε k ({R n }), c k (G) is then given as solution to This defines a continuous basis with eigenvalues related to momenta k which are now used as input in equation (27) or (29). As the simplest example, the small polariton is given as where the momenta are connected according to Within the long-wavelength approximation, that is, there is no momentum-transfer between photon and electron, the electronic momentum is a conserved quantity up to Umklapp scattering. Note that ∇ kk = 0 for real functions. As a consequence, the quantized photonic nature splits the bare electronic excited bands into polaritonic bands, that is, similar to the interpretation of upper and lower polariton as for a two-level system in Sec. IV B.
Solving this equation will already imprint first quantum features in periodic materials from a natural and simple perspective. We want to emphasize here, that this does not describe the correct many-body excitations since each of the Kohn-Sham orbitals is coupled to the field instead of many-body states as the Kohn-Sham Slaterdeterminants. The above procedure is consequently the one of an effective single particle one. Calculations for realistic two-dimensional systems are under investigation [39]. The correct description would involve coupling of many-body eigenstates corresponding to poles of the spectral function [15].