Ab-initio Optimized Effective Potentials for Real Molecules in Optical Cavities: Photon Contributions to the Molecular Ground state

We introduce a simple scheme to efficiently compute photon exchange-correlation contributions due to the coupling to transversal photons as formulated in the newly developed quantum-electrodynamical density functional theory (QEDFT). Our construction employs the optimized-effective potential (OEP) approach by means of the Sternheimer equation to avoid the explicit calculation of unoccupied states. We demonstrate the efficiency of the scheme by applying it to an exactly solvable GaAs quantum ring model system, a single azulene molecule, and chains of sodium dimers, all located in optical cavities and described in full real space. While the first example is a two-dimensional system and allows to benchmark the employed approximations, the latter two examples demonstrate that the correlated electron-photon interaction appreciably distorts the ground-state electronic structure of a real molecule. By using this scheme, we not only construct typical electronic observables, such as the electronic ground-state density, but also illustrate how photon observables, such as the photon number, and mixed electron-photon observables, e.g. electron-photon correlation functions, become accessible in a DFT framework. This work constitutes the first three-dimensional ab-initio calculation within the new QEDFT formalism and thus opens up a new computational route for the ab-initio study of correlated electron-photon systems in quantum cavities.

O ver the past decades, methods in computational material science and quantum chemistry have been successfully applied to accurately model material properties. Such material properties usually depend on the electronic structure of the system of interest that is dictated by the laws of quantum mechanics. Recently it has been demonstrated that by hybridizing light strongly with the electronic structure of the system, novel effects appear providing a promising route for a new design of material properties. Such recent experiments include, matter−photon coupling for living systems, 6 vibrational strong coupling, 7 changes in chemical reactivity, 8 symmetry protected collisions of strongly interacting photons, 9 the Bose−Einstein condensation 10 or the room-temperature polariton lasing 11 of exciton-polaritons, and ultrastrong coupling in circuit-QED 12 to mention a few. Condensed matter systems driven out of equilibrium provide optional possibilities for novel properties, for example, the creation of Floquet-Bloch states 13 and Floquet-Weyl semimetals. 14 Additionally, the Floquet-scheme enables the development of new time-dependent DFT functionals with explicit memory dependence. Recently, we and our co-workers have introduced a novel density-functional approach (QEDFT) to describe such complex dynamics of strongly interacting electrons, photons, and phonon systems, [1][2][3][4][5]15 all on the same theoretical footing.
The framework of QEDFT is the first attempt to deal with the electron-photon interaction from first-principles and has been demonstrated for the first time in refs 1, 2, 5, and 16. Together with new experiments on chemical systems in optical cavities, 7,8,17,18 this work now opens up the field of quantumelectrodynamical (QED) chemistry and QED materials. 8,19,20 In this new field, so far model Hamiltonian schemes have also been used to successfully describe experiments, 21,22 but for an ab initio description a full real-space picture is necessary. QEDFT additionally allows to study multimode effects 19 that have been recently observed in experiment 23 and theory. 24 As in any density-functional theory, the practical applicability of QEDFT is build upon the underlying approximations for the exchange-correlation (xc) functional. In contrast to traditional density-functional theory, 25 where a whole range of different approximation schemes for the xc functional are available, 26 QEDFT still lacks a practical method to construct such approximations. Previous works 1,5,15 have opened the path to the development of such exchange-correlation functionals. Different routes are possible, for example, functionals based on, for example, the electron density, the electronic orbitals, or the electron current, 27 ultimately leading to the first quantitative accurate semilocal QEDFT functional. To close the gap, in this work, we introduce a simple, yet accurate, computational scheme to calculate the ab initio xc potential for electronic systems coupled to quantized photon modes based on the occupied electronic orbitals. This method is based on the optimized effective potential (OEP) approach introduced by some of us in ref 15. OEP has been previously used for purely electronic systems in DFT. 28−30 In ref 15, the construction of the OEP functional relies on the calculation of occupied and unoccupied orbitals. In particular, the calculation of unoccupied states is computationally very demanding due to the large configuration space in any realistic many-body simulation and therefore hampers the practicability of the scheme. Here, we introduce a scheme that overcomes this limitation and does not involve the calculation of any unoccupied orbital. As a consequence, we find our scheme to be numerically highly efficient and thus we are able, for the first time, to calculate realistic molecular systems interacting with quantized photon modes from first principles. To achieve this goal, we employ the Sternheimer scheme 31 that allows us to construct the electron− photon OEP equation in an efficient manner. In this way, we only require the calculation of occupied orbitals together with solving linear equations that makes this approach computationally superior to the previous formulation. As a consequence, our proposed scheme fits within the definition of a Kohn− Sham (KS) DFT as proposed by Axel Becke 32 that defines KS-DFT as occupied-orbital-only. Similar schemes have been used in the context of density-functional perturbation theory 33 and in many body-perturbation theory using the GW self-energy approach. 34 This paper is structured as follows: in section II, we introduce the formal framework to construct the ground-state xc potential using the OEP scheme. In section III, we apply the scheme to three different numerical examples and demonstrate the accuracy and the numerical feasibility for large-scale calculations. In the first example, we employ a model system for a GaAs quantum ring. 2,35 For this example, which is also exactly numerically solvable, we assess the accuracy of the scheme. We identify limiting cases when to expect reliable results from the approximation. In the second example, we apply our method to a three-dimensional system, the azulene molecule in full three-dimensional real space. We demonstrate the effects of the correlated electron-photon interaction on the ground-state density. Additionally, we construct a mixed electron-photon correlation function that illustrates necessary ingredients for novel correlated electron-photon spectroscopies. The last example of this paper treats chains of sodium dimers that allow us to systematically study the effects of many molecules in optical cavities. The latter two examples are the first QEDFT calculation for realistic molecules. All these calculations demonstrate the reliability and applicability of the proposed numerical scheme. With realistic systems now computationally accessible, a promising avenue in the design of QED materials is introduced.

■ THEORY
We start by introducing the general nonrelativistic setup of the correlated electron-photon systems considered in the present work following previous works. 1,5,15,16 Let us consider N e = ∑ σ=↑,↓ N σ = N ↑ + N ↓ interacting electrons of spin ↑ or ↓ located in an optical cavity. The electrons are coupled to N p quantized electromagnetic modes, that is, photon modes. Each photon mode is identified by its cavity frequency ω α and polarization direction e α . We describe the matter−photon coupling in the Coulomb gauge, dipole approximation (long-wavelength approximation) and the length gauge. 5,36 The hereby emerging electron−photon coupling strength parameter λ α is projected on the photon polarization direction λ α = e α λ α . While in Coulomb gauge, the matter−photon interaction is explicitly described by the transversal degrees of freedom, the longitudinal degree of freedom leads to the Coulomb potential that describes the two-particle electron−electron interaction 1/|r i − r j |. In this setup, the total electron−photon Hamiltonian reads (in atomic units) 5,15 ∑ ∑∑ where each photon mode is associated with a bosonic creation and annihilation operator (aα † , aα) that creates and destroys photons in the mode α. The transversal electron−photon interaction Ĥe p (α) consists of two terms that read explicitly In dipole approximation, the electron−photon coupling comprises the electron dipole operator R = R 0 + ∑ i=1 N e r i and the photon displacement coordinatê=̂+α The electronic coordinates r i are defined with respect to the center of charge of the system R 0 . As has been outlined in earlier work, using the creation and annihilation operators, we can setup the photon displacement and photon momentum operators qα and pα. 2 Physically these two operators are directly connected to the electric displacement field and the magnetic field, if summed over all modes. 2,19 The electron−photon coupling strength is given by where S α denotes the mode function, for example, a sinefunction for the case of a cubic cavity 1,15 and k α is the wave vector. The effect of the nuclei employing the frozen-nuclei approximation enters the electron-photon Hamiltonian of eq 1 via the external potential v ext (r). The effect of a static permanent dipole moment due to the nuclear charge can be neglected, since the two terms of eq 2 compensate each other in that case. For nuclei effects beyond the frozen-nuclei approximation, we refer the reader to ref 19, where electrons, nuclei, and photons are treated on the same quantized footing. Comparing QEDFT to standard DFT, we note that in QEDFT we have two sets of internal variables, that is, the electron density n(r) and the photon displacement variables q α . It can be shown 3 that these two internal variables are in an oneto-one correspondence to the external variables v ext (r) and j ext (α) . Here j ext (α) corresponds to the first-order time-derivative of an external charge current at time zero projected on the mode α, that is, ∫ d 3 rS α (k α ·r)e α ·∂ t j ext (r,t) at t = 0. 1,5 The reason for the appearance of the time-derivative is the length-gauge transformation that rewrites the coupling to the photons in terms of the displacement field instead of in terms of the vector potential. 1,5 Since the displacement field corresponds to the electric field minus the polarization, 19,37 and the electric field is

ACS Photonics
Article the time derivative of the vector potential, the conjugate external variable to q α needs to contain a time-derivative as well. Since we can unitarily transform (by polarizing the photon vacuum) the Hamiltonian with j ext α ≠ 0 into one with j ext α = 0, we can construct by the very same simple transformation the solutions for the inhomogeneous case from the case with j ext α = 0 (see ref 38 for details). By exploiting the one-to-one correspondence of QEDFT, we find that all observables (electronic, photonic, and mixed) become functionals of the internal variables. Formulated differently, any change in the internal variables will lead to changes in experimentally accessible observables.
In this work, we use the KS scheme 25 of density-functional theory introduced for electron-photon problems in refs. 1,2,5 and commonly used in all DFT calculations. The KS scheme allows us to describe interacting many-body problems by the following set of effective noninteracting equations 5 for N σ Kohn−Sham orbitals φ iσ (r) with spin σ. The effective Kohn−Sham potential v sσ (r) is given by and can be divided into the external potential v ext (r), the usual Hartree-exchange-correlation (Hxc) potential v Hxcσ (r) that accounts for the electron−electron interaction and the modedependent meanfield-exchange-correlation potential (Mxc) v Mxcσ (α) (r). (In general electron−photon systems, we find that contributions due to the kinetic energy can not be attributed solely to the electron−electron or electron−photon interaction.) Both Hxc and Mxc contain the unknown exchangecorrelation parts that have to be approximated. In exact KS-QEDFT, these parts are chosen such that the electron density n(r) that is the sum of the spin-resolved electron densities n σ (r) = ∑ i=1 N σ φ iσ *(r)φ iσ (r) is equivalent in the interacting and the noninteracting system. In the ground state, we have a simple connection between the exchange-correlation energy that includes contributions from the electron−electron interaction (ee) and the electron−photon interaction (α) and the corresponding xc potential that reads as follows 30 This connection will be now exploited to setup the OEP equation. Throughout this work, we use the exchange-only approximation, that is, While we use the standard definition for E x (ee) , 28,29 that is, the Fock energy, we focus in the following on the exchange energy due to the electron−photon interaction E x (α) . The interaction terms in eq 2 contain the electron−photon coupling strength λ α in first-order and second-order. For the ground state the first-order exchange energy vanishes, 15 if the photons are not exposed to an external current j ext (α) . Therefore, the leading order becomes the secondorder in λ α . While the second part of eq 2 (the dipole selfinteraction part) is time-local just as in the typical Coulombic exchange, the explicit electron−photon interaction part involves the absorption and emission of a single photon, and therefore, the propagation of a single photon state that generates a frequency dependency of the corresponding electronic self-energy. As a consequence, the exchange energy can be written as an orbital functional as 15 where c.c. refers to the complex conjugate of all former terms. Additionally, we define the projected dipole operator dα = λ α ·r.
As does the electron−photon interaction Hamiltonian in eq 2, also the electron−photon exchange energy E x (α) consists of two parts, both containing different electronic orbital shifts. The first orbital shift is the solution of the following Sternheimer equation 1 ( ) (9) with the matrix element d ijσ (α) = ⟨φ iσ |dα|φ jσ ⟩ and the orbital shift can be written explicitly as The second orbital shift is defined by While both orbital shifts can be formulated explicitly in terms of all KS orbitals (in eqs 10 and 11, respectively), only the second orbital shift Φ iσ,α (2) can be formulated explicitly in terms of solely occupied orbitals given by eq 12. However, the shift Φ iσ,α (1) can be defined implicitly by a Sternheimer equation that only invokes occupied orbitals as given in eq 9.
Since the exchange energy given in eq 8 scales with λ α 2 , the exchange energy is the Lamb shift of the ground state. 15 Thus, all corrections for the ground state are in their magnitude on the order of the Lamb shift. For electron-photon problems, we find that E x (α) , as defined by eq 8, has a functional dependency on all occupied orbitals and both orbital shifts. The standard route to obtain the OEP equation involves the calculation of functional derivatives of the orbitals and accordingly has to be generalized for the electron-photon case. In this case, we need consequently also the functional derivatives of the orbital shifts.
Nevertheless, as will be demonstrated in the following, the standard route to construct the OEP equation can be adapted to accommodate this subtle difference. Having defined the total exchange energy E x in eq 6, we now proceed to calculate the corresponding Kohn−Sham potential using functional derivatives. From eq 7, we can setup the following OEP equation by using the chain rule of functional derivatives

ACS Photonics
, , (2) p (13) The OEP equation of eq 13 contains several different terms that need an individual point-wise evaluation. First, we start discussing the functional derivatives of the exchange energy.
These terms can be calculated straightforwardly using eq 8 and are given as follows (please note that, for brevity, we do not explicitly evaluate the E x (ee) contributions, but state its As the next step, we need to calculate the functional derivatives of the KS orbitals and orbital shifts with respect to the Kohn−Sham potential v s . In the case of the KS orbitals, this derivative is given by 29 (17) where the sum runs over all orbitals, except i = j. All remaining terms in eq 13 are functional derivatives of the orbital shifts. We start by discussing Φ iσ,α (2) (r), since it is conceptually simpler to obtain, than Φ iσ,α (1) (r). From eq 12, for an infinitesimal change in (r), by keeping only first-order terms and combining with eq 17, we obtain The derivation of the remaining functional derivative of the first orbital shift, that is, δΦ iσ,α (1) (r)/δv s (r′) is given in full detail in the appendix and we only state the final result here Combining all these terms brings us to an alternative formulation of the OEP equation. By now plugging all ingredients into eq 13 an alternative OEP equation can be derived that is given by the simple equation 20) with the Kohn−Sham Green's function 29 Due to the energy dependence of E x (α) , we find that the nonvanishing additional inhomogeneity 30 Λ iσ (r) is given by , (22) and the orbital shift M iσ (r) by The orbital shift M iσ (r) contains in the first line the electron−electron interaction, we choose the exchange-only approximation, that is, and E x (ee) is the usual Fock exchange energy. The following lines are corrections due to the correlated electron-photon interaction that induce density changes in the electronic system. 15

Article
The main advantage of the present reformulation is that we can write the OEP equation for electron-photon problems in a simple form. This formulation is similar to refs. 28,30 that provide the formulation for electrons-only. (25) and the orbital shifts ψ iσ *(r) can be obtained using a Sternheimer equation This equation has to be solved self-consistently with eq 9. By this procedure, we have replaced the problem of calculating the OEP equation using all unoccupied states by a problem of solving N p + 1 Sternheimer equations that only involve occupied orbitals.

■ NOVEL TYPES OF OBSERVABLES
One of the advantages of QEDFT over DFT is the correct treatment of the quantum nature of the photon field and its interaction with a correlated many-body electron system. Thus, by exploiting the one-to-one correspondence of the internal variables to the external variables, 3 the photon observables (and the electronic observables) become functionals of the internal variables. Therefore, if we know the internal variables and their functional dependency, we can construct arbitrary observables. In the case of orbital functionals, we can use the KS orbitals to construct these observables. In this section, we now introduce new types of observables into the DFT framework, that is, photonic observables and observables of mixed matter−photon character. The first example for a photonic observable is the number of photons in each mode. This observable can be calculated in terms of KS orbital shifts as follows Physically, we can attribute the orbital shifts that are calculated by eq 12 with wave function corrections that carry each a single photon. Thus, we can use these shifts to calculate the photon number in each photon mode. While the first term in eq 27 is due to the quantum fluctuations of the photon mode, the latter term is a classical contribution due to a nonvanishing R 0 . By performing this connection, we assume that the photon number is dominated by contributions originating from single-photon processes. To this end, we can expect a good quality of this photon number observables if this is the case, while if many-photon processes contribute we expect poorer results.
Examples for mixed electron-photon observables 20 are electron-photon correlation functions. For instance, the charge-density-displacement-field correlation function A (α) (r) we define as where Ψ 0 denotes the many-body ground state of the system. The given expression is the leading term in an expansion in orders of λ α . Physically this correlation function corresponds to the local forces that the displacement field of the photons exerts on the electrons. 2,5 If we perturb the photon field, the change of these local forces will rearrange the charges in an intricate manner. While for a classical field A (α) (r) merely becomes a product of the (positive) electronic density and the value of the displacement field and is therefore only positive or negative, in the quantum case, A (α) (r) can locally change its sign. Consequently, probing this correlation function spectroscopically could allow to obtain novel insights into structural properties of complex systems.

■ KRIEGER-LI-IAFRATE APPROXIMATION
As will be demonstrated in the application section, the OEP equation leads to accurate results. However, since the xc potential is given only implicitly by eqs 25 and 26, it may be hard to converge. One way to circumvent this problem and to obtain an explicit formula for the xc potential is the Krieger-Li-Iafrate (KLI) approximation. 39−41 In contrast to the common Coulomb OEP equation, 28 in the case of correlated electron−photon coupling, an additional inhomogeneous contribution appears, that is, Λ iσ . The consequence of this structural deviation from the well-known OEP equation in the electronic case, where no inhomogeneity is present, complicates the common approximation schemes. A direct energy denominator approximation does not only leave an arbitrariness on the remaining energy denominator but the corresponding approximations leave divergent contributions uncanceled. The reformulation in terms of Sternheimer shifts avoids unbalanced approximations in divergent contributions. If we multiply eq 25 with the Kohn−Sham potential 28 and decompose eqs 25 and 26 starting from with eq 23 and we arrive at the exact reformulation The common approximation scheme now assumes (hŝ σ (r) − ε iσ )ψ iσ *(r) = 0, which is exact for a single electron if no inhomogeneity would be present in eq 25. A corresponding substitution involving ψ iσ *(r) ≈ Λ iσ (r)/φ iσ (r) leads in the general case to nodal points. The variety of possibilities result in different deficiencies and inconsistencies (see also Engel et al. 42 ). To remain as consistent as possible we decide to assume (hŝ σ (r) − ε iσ )ψ iσ *(r) = 0 and the KLI equation reads then

ACS Photonics
This reformulation avoids the solution of eq 26 and can be solved explicitly for the exchange potential as a linear equation. This improves the stability with respect to the initial guess and represents in many cases a valuable starting-point for the OEP calculation. The KLI effectively neglects off-diagonal contributions to the response function mediated by the exchange potential. Connecting to this, the accuracy reduces with increasing local dipole-moment which will especially manifest in the overestimation in local density perturbation under cavity influence.

■ NUMERICAL DETAILS
We have implemented the OEP equation of eq 25 and the corresponding KLI equation of eq 33 in the OCTOPUS package. 43−45 The OEP equation can be solved using standard algorithms as, for example, described in ref 28, that is, in a selfconsistent field (SCF) cycle. To obtain the numerical algorithm, we reformulate eq 25 as follows The quantity S σ (r) becomes a measure for convergence, since it is vanishing in the case of convergence (compare eq 34 and eq 25). To obtain the new potential in the SCF step, we use Different schemes to calculate c(r) are possible, 46 for example, dividing by the electron density, 47 using the Barzilai-Borwein minimizer 48 or connecting to conjugate-gradient algorithms. 46 However, for our purpose, we find that choosing a constant value 28 is already sufficient and already provides the most stable and reliable algorithm. Thus, we choose c(r) = 0.1 au for the azulene molecule and c(r) = 20 au for the sodium chains.
As in the case of electronic OEP, 28,41 we also find for the photon OEP that we can add an arbitrary (spatial-independent) constant to the exchange potential that does not alter the physical results. If we follow the lines of the electronic OEP 28,41 and enforce the condition N , we find that in the single electron case, the single particle Kohn−Sham energy deviates from the total energy. From a physical point-of-view it is desirable that both coincide to connect to ionization energies. We find by fixing to the contribution of the highest occupied orbital to the exchange energy defined in eqs 6 and 8, that is, , we can restore this connection. However, we note that for the electronic OEP 28,41 both routes coincide. Furthermore, since in the present study we focus on energy differences, the arbitrary constant only modifies the offset in the presented xc potentials.

■ NUMERICAL APPLICATION
As numerical applications, we analyze different examples in 2D and 3D. The first example is used to demonstrate the accuracy of the employed method. To this end, we benchmark the OEP scheme with an exactly solvable model system, that is, a GaAs quantum ring located in an optical cavity, 2,49 where we have published exact results previously. 2,49 In this way we are able to validate the presented scheme before in the next examples, we apply it to real systems. Thus, in the second example, we solve the electron−photon OEP equation for the first time in full three-dimensional real space. We study the azulene molecule and report the changes in the ground-state density if the molecule is located inside an optical cavity. The last example deals with realistic ensembles of molecules in optical cavities.
Here we study the ground-state density of chains of sodium dimers with different length. The different examples studied in this work are schematically depicted in Figure 1.

■ GAAS QUANTUM RING IN AN OPTICAL CAVITY
We start by discussing the model for a GaAs semiconductor quantum ring coupled to a single photon mode in an optical cavity. The model consists of a single electron restricted to two spatial dimensions in real-space (r = r x e x + r y e y ) interacting with the single photon mode with frequency ℏω α = 1.41 meV and polarization direction e α = (1,1). The polarization direction enters via the electron−photon coupling strength, that is, λ α = λ α e α and depends on the specific experimental setup. We choose the photon mode frequency in resonance with the first electronic transition. The external potential of the single electron is given by the following formula with parameters ℏω 0 = 10 meV, V 0 = 200 meV, d = 10 nm, and m 0 = 0.067m e . 35,49 For the electron−photon coupling strength, we choose two values, that is, in weak coupling with λ α = 0.0034 meV 1/2 /nm and in strong coupling λ α = 0.1342 meV 1/2 / nm. The effective three-dimensionality of this problem (twodimensional electron and one-dimensional photon mode) is low enough that an exact solution is still accessible via exact

ACS Photonics
Article diagonalization. 50 To obtain the exact ground state, we employ a two-dimensional grid of N = 127 grid points in each direction with a grid spacing of Δx = 0.7052 nm to describe the single electron. The photon field is represented in the photon number eigenbasis and we include up to 41 Fock number states. Using the exact wave function, we can numerically construct the exact exchange-correlation potential. 2, 51 We start by discussing the weak-coupling limit, where λ α = 0.0034 meV 1/2 /nm. In Figure  2a, we show the exact ground-state density obtained by exact diagonalization. Compared to the bare electronic ground-state (for λ α = 0) that also has a ring structure in the weak-coupling limit, we find small distortions. 49 In Figure 2b, we show the OEP ground-state density and in (c) the difference between the exact and the OEP ground-state density. The deviation of the OEP ground-state density to the exact ground-state density is very small and in the order of magnitude of 10 −10 , that is, close to our numerical precision. This high precision of the approximate electron density has its origin in the high quality of the OEP approximation for the xc potentials. The exact and the OEP xc potential are plotted in (d) and (e), respectively. In (f) we plot the difference of the exact to the OEP potential and find significant differences ( ∼ − 10 5 ) only in low-density regions, that is, at the border of our grid. In contrast, the inner high-density regions are well approximated leading to the accurate description of the electron density. This larger error can also be attributed to the algorithm, since low density regions are harder to converge in the OEP scheme. However, since low density regions do not contribute much to observables such as the total energy, this error will effectively not influence the overall result. In Figure 3 we show how the KLI approximation (Sec.) performs in the weak-coupling limit for the single-electron case. In (b) we plot the KLI approximated electron density and in (c), we show the difference to the exact reference. We find errors in the electron density in the order of ∼ − 10 7 that are due to the KLI xc potential. The KLI xc potential is shown in (e). We find that in comparison to the exact reference shown in (d) the overall shape of the potential is approximated correctly, while the peak in the middle of the potential is missing. The deviations can be also seen in (f), where we plot the difference between the exact and the KLI xc potential. To quantify the differences for this example, we print the results of our calculations in Table 1. The first three rows show the exact, OEP and KLI results for the total energy E tot and the photon number n pt in the weakcoupling limit using the external potential of eq 37. Overall, we find a very good performance, of the OEP and KLI approximations. The OEP performs slightly better, but also the KLI gives accurate energies and photon numbers. Let us now analyze the strong-coupling limit. In Figure 4, we show the results obtained in the strong-coupling regime (λ α = 0.1342 meV 1/2 /nm), where we find a deviation in the exact groundstate electron density from the ring structure in the weakcoupling regime to a double-well structure 19 as shown in Figure   Figure 2. Exact (a) and OEP approximated (b) electron density in the weak-coupling limit (λ α = 0.0034 meV 1/2 /nm). The difference is shown in (c). The corresponding xc potentials are shown in (d) and (e), respectively, and (f) shows the difference in the xc potentials.

ACS Photonics
Article 4a. This splitting is accompanied by a higher peak in the xc potential in the center of the grid as shown in Figure 4d. Although in the weak-coupling limit, we find a very high accuracy of the OEP approximation, in the strong coupling limit, we observe the breakdown of the OEP approximation. In Figure 4b, we find that the OEP predicts an electron density that is located in only one of the two subwells with a wrong xc potential shown in Figure 4e. Consequently, the error of the OEP density and the potential shown in Figure 4c,f are very high. The origins of this failure of the OEP can be understood by calculating the photon number ⟨â †a⟩ and the double occupancy ⟨â †â †aâ⟩ in the photon mode shown in Figure 5.
Scaling the electron-photon coupling strength λ α from the weak to the strong coupling limit, we find that two-photon processes become the dominant contributions to the ground state, when the electron density splits. 49 Since the OEP approximation by construction only considers single photon processes, its failure in this region is a natural consequence of the higher weight of two (and more) photon processes in the setting of the xc potential. In ref 49 we have calculated the exact eigenvalues and find a close degeneracy of the ground state and the first-excited state in the strong-coupling limit (reminiscent of static correlation in quantum chemistry 52 ). From a numerical point of view, this degeneracy introduces an instability in the selfconsistency procedure. Similarly as in the electron-only case, where static correlation can be described by including correlation effects beyond exact exchange, in correlated electron−photon problems, we conclude that in the strong coupling limit going beyond exact exchange, that is, single photon processes, to higher order processes, that is, twophoton, three-photon, and so on is required to accurately describe this limit. In the last example, we study an asymmetric example in the strong-coupling limit (λ α = 0.1342 meV 1/2 /nm), where the external potential is given bỹ with V̅ 0 = 0.1123 meV/nm. The cavity frequency is again chosen to be in resonance to the first electronic excitation, that is, ℏω α = 2.72 meV. The results are shown in Figure 6. We find while the density is approximated accurately, with errors in the order of 10 −6 , also observables such as the photon number listed in Table 1 are approximated quite accurately, since in this regime the mean-field contribution in eq 27 becomes dominant in comparison to the fluctuations.
In conclusion, we have demonstrated in this section that the photon OEP is capable of describing a wide range of parameters correctly. In the weak-coupling regime, we have found highly accurate results. Additionally, we find in the strong coupling regime a failure of the OEP in the symmetric setup, while in the asymmetric setup, we have again an accurate description of the electron density. Having at hand a scheme to derive approximations for general functionals, we can also investigate novel types of observables that are not accessible with traditional DFT but need a QEDFT calculation. In the case at hand we find, for instance, good agreement for the photon number, where both the OEP and KLI approximation yield reliable results. Next, after we have assessed the quality of our approximations, we turn to real systems and show that QEDFT is an efficient ab initio scheme to determine properties of complex systems coupled to photons.
■ AZULENE (C 10 H 8 ) MOLECULE IN AN OPTICAL CAVITY Our next example is the first real application of the QEDFT framework to a three-dimensional molecule, that is, the azulene (C 10 H 8 ) molecule (Figure 7). To find a reliable equilibrium structure and determine the cavity frequency, we follow the following route. First, we obtain the 3D conformer structure for azulene from the PubChem database 53 (CID: 9231). Second, we use the geometry optimization of the OCTOPUS package

ACS Photonics
Article employing the LDA functional 25,54 to calculate a relaxed ground-state structure. Third, using this relaxed structure, we use the electron OEP to calculate a HOMO−LUMO gap of 2.41 eV that serves as the cavity frequency, that is, ℏω α = 2.41 eV. (The ground state results are not sensitive to a resonance. 55 ) The electron−photon coupling includes the polarization direction of the photon field that is polarized along the x-direction with a strength of λ α = 37.47 eV 1/2 /nm (0.08 au), that is, λ α = 37.47 eV 1/2 /nm e x . [For standard experimental parameters, e.g., for a single trapped-atom cavity, as described in ref 56 (Figure 3, V = 18.148 μm 3 ), we deduce an experimental value of λ 0 = 3.9 × 10 −7 eV 1/2 /nm (8.32 × 10 −10 a.u.).] In this example, we want to investigate the question how the correlated electron−photon interaction alters the electronic ground-state density n 0 (r). To numerically calculate the ground-state density of the azulene molecule, we use a grid of dimensions 32 × 36 × 16 Bohr in xyz directions. The grid spacing is chosen to be Δx = 0.11 Å and to describe the core electrons of the carbon and hydrogen atoms we use LDA Troullier-Martins pseudopotentials. 57 Thus, we explicitly describe the 48 valence electrons in our calculations amounting to 24 doubly occupied Kohn−Sham orbitals. As described in the previous section, to describe the electron−electron interaction, we use the Fock exchange energy 28 also in the OEP setting. In Figure 8, we show in the top panel the groundstate density of the molecule in an optical cavity as a cut in the x−y plane. The electrons are highly localized in-between the nuclei. The aromatic ring structure induced by the arrangement of the carbon atoms is inherited in the ground-state electron density that has naturally the same symmetry. The middle plot of Figure 8 shows the difference in the electronic ground-state density exposed to electron−photon coupling to the bare electronic ground-state density, that is, the change in the density by going from gas phase to the case inside the cavity. The figure shows a rich fine structure in the center of the molecule, but also a pronounced accumulation region of electronic density at the top and bottom rim of the molecule. The plot on the bottom of Figure 8 show the results of the KLI approximation. While the KLI seems to fail to correctly describe the rich inner structure of the density differences Δn(r), it correctly predicts the density accumulation regions at the top and bottom of the molecule. However, these regions are overestimated by a factor ∼4. To quantify the effects of the quantized electron−photon interaction on many-electron systems, we have provided numerical results in Table 2. For different levels of theory, we print the energy difference between lowest and highest occupied orbitals (24−1), the HOMO−LUMO gap (25−24), the total energy E tot , and the electronic and the electron−photonic part of the exchange energy E x (ee) and E x (α) , respectively. For the given parameters, the electron-photon exchange energy is in the order of ∼3.8 eV and 2 orders of magnitude smaller than the electronic exchange energy E x ee that is roughly ∼500 eV. As could be expected, the changes due to the coupling to the vacuum of the cavity are small in the ground state, that is, we have determined the Lamb shift. However, due to the electron−photon coupling we now have access to novel types of observables. To connect to the novel mixed electron−photon observables within the framework of QEDFT, we calculate the correlator A (α) (r), as defined in eq 28, without the mean-field contributions in Figure 9. We find that the resulting local-force map due to the coupling to the photons indeed shows a complex structure with local sign changes. It indicates the forces that the electrons experience due the displacement field. Indeed, the local forces nicely agree with the rearrangement of the charge density upon coupling to the vacuum field. If we would perturb the photon mode, the electrons would experience forces in different directions depending on their position in the molecule. In contrast, a classical field in dipole approximation would only induce forces in one direction. In conclusion, in this section we have presented the distorting effects of the quantized electron- Figure 8. From top to bottom as a cut in x−y plane: OEP ground-state density of azulene, difference of electron−photon OEP and electron OEP ground-state density, and difference of electron−photon KLI and electron KLI ground-state density.

ACS Photonics
Article photon interaction on molecules in cavities. We find that in QEDFT new observables become numerically accessible that could allow for novel experimental spectroscopic setups. 20

■ CHAIN OF SODIUM DIMERS
The last example that is studied in this paper is a chain of sodium dimers of variable length, that is, up to 10 sodium dimers. We use this setup to highlight that QEDFT allows to investigate problems of quantum optics from first principles. For instance, we can consider the reliability of the ubiquitous Dicke model, 58 where many two-level systems are coupled to a cavity mode. This model predicts that due to the collective behavior of the two-level systems the usually weak coupling of the matter to the photon mode is effectively increased. This collective effect is one way of reaching the strong coupling limit in experiment. Still, due to the many simplifying assumptions employed in deriving this model some implications are debated, for example, the super-radiant phase transition. 59,60 With a firstprinciples approach such as QEDFT many of these assumption can be avoided which could shed new light on these issues. Here we will not target these more intricate problems but rather show that we can recover from first-principles the increase in the effective coupling strength. We do so by analyzing the behavior of the correlated electron-photon ground-state, when more and more emitters are coupled to the cavity field For this setup, we use the parameters for a sodium dimer given in ref 61. For the optical cavity frequency, we choose the energy of the 3s−3p transition, that is, ω α = 2.19 eV. We assume that the photon field is polarized along the direction of the sodium chains with a strength of λ α = 2.95 eV 1/2 /nm (0.006 au), that is, λ α = 2.95 eV 1/2 /nm e y . To calculate the chain of sodium dimers (Na 2 ), we use the sodium pseudopotentials and equilibrium distances for a single sodium dimer of ref 61. For the real-space grid, we use dimensions 60 × min(60, 2N c × 10) × 60 Bohr with grid spacing 0.5 Bohr, where N c is the chain length. The distance between the sodium dimers is chosen as d = 10 Bohr. The case of three sodium dimers is illustrated in Figure 10. As in the previous example, in the top panel we show a cut of the ground-state electron density in the x−y plane. Each sodium dimer contains two electrons and the electrons are localized between the sodium nuclei. In the middle plot, we show the difference in the electron density of the system with and without the cavity. The lower plot in Figure 10 shows a cut along the y-axis in blue against the ground-state electron density in the cavity in gray. We find three maxima for density accumulation and four minima from where the density has been rearranged. Further, we find that the electron-photon interaction pushes the electron density onto high-density regions. This density accumulation stems from regions inbetween the dimers, where the amount of density is decreasing in the cavity. The next figure, Figure 11 shows the case of 10 sodium dimers. The first plot shows the electron density of the ground state. In the second plot we see the difference of the electron density of the system inside the optical cavity to the bare electron density. Between the maxima, we find local minima where electron density is rearranged, as shown in the plot in the bottom. We compare to the KLI approximation in Figure 12.
Here we find the KLI strongly overestimates the effects of the electron-photon interaction. In the last figure of this section, Figure 13. We plot the number of photons in the correlated electron-photon ground state using the functional presented in eq 27. In total, we find for the KLI and the OEP a linear behavior, where the KLI overestimates the number of photons slightly. From eq 27 we also see that ⟨aα † aα⟩ ∼ λ α 2 . Thus, we can capture this behavior alternatively by defining a new effective Figure 9. Correlation function as a cut in x−z plane A (α) (r) as defined in eq 28, calculated for the azulene molecule. Figure 10. From top to bottom as a cut in x−y plane: OEP groundstate density of three sodium dimers, difference of electron-photon OEP and electron OEP ground-state density, and sum along the x-axis of the difference between the electron-photon OEP density and the electron OEP density in blue against the electron-photon OEP density in gray. Please note that the latter has been reduced by a factor of 1/ 2000.

ACS Photonics
Article coupling constant λ̃∼ α N c that scales with the square-root of the chain length. This example nicely illustrates the collective coupling of matter to the cavity field in the weak-coupling regime. This result agrees with predications based on the Dicke model, where the coupling strength scales with the square root of the number of two-level systems. However, still more work needs to be done to properly characterize the emergence of collective phenomena due to the strong light-matter coupling in a set of N emitters.

■ SUMMARY AND CONCLUSION
In conclusion, in this work, we have illustrated how the cavitymediated electron-photon interaction is capable of rearranging the electron density in two-and three-dimensional systems. We find that our OEP approach accurately describes situations in the weak coupling limit. In the strong coupling limit, we find broken symmetry solutions which can be attributed to the accuracy of the employed approximate transversal energy orbital functional. The overall effect of the transversal photons on the ground state density is minor as expected from the magnitude of the Lamb-shift-type-energy correction. However, it allows to investigate problems of quantum optics from firstprinciples, such as the increase of the effective matter-photon coupling strength upon increasing the number of molecules inside a cavity. Furthermore, the present work lays the foundation for the ab initio construction of excited states and Figure 11. From top to bottom as a cut in x−y plane: OEP groundstate density of ten sodium dimers, difference of electron-photon OEP and electron OEP ground-state density, and sum along the x-axis of the difference between the electron-photon OEP density and the electron OEP density in blue against the electron-photon OEP density in gray. Please note that the latter has been reduced by a factor of 1/ 2000. Figure 12. From top to bottom as a cut in x−y plane: KLI groundstate density of ten sodium dimers, difference of electron-photon KLI and electron KLI ground-state density, and sum along the x-axis of the difference between the electron-photon OEP density and the electron OEP density in blue against the electron-photon OEP density in gray. Please note that the latter has been reduced by a factor of 1/100.

ACS Photonics
Article new functionals for QEDFT. While the small contribution of transversal photons on the ground state is reaffirming standard DFT calculations that neglect coupling to transversal photons, we have found that the effect of transversal photons on excited states, such as, for example, Rabi splittings, can be substantial. The present work constitutes the first mandatory step toward such studies of excited states of strong light-matter coupled quantum systems. Work beyond the exchange approximation, i.e. to include multiphoton processes is currently under investigation. Additionally, this approach could also benefit standard electronic DFT, since similar ideas, that is, expressing the exchange-correlation energy in terms of orbital shifts could also be applied to the correlation part in the xc approximation.
We have introduced our QEDFT approach as a viable tool to predict and describe the emerging field of QED chemistry, and QED materials, where chemical systems are placed in optical cavities.
■ APPENDIX Derivation of the Functional Derivative of Second Orbital Shift