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

###### Abstract

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) Ruggenthaler et al. (2014); Flick et al. (2015); Ruggenthaler (2015); Ruggenthaler et al. (2011); Tokatly (2013). 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.

## I Introduction

Over the past decades, methods in computational material science and quantum chemistry have been successfully applied to accurately model materials 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 Coles et al. (2017), vibrational strong coupling George et al. (2016), changes in chemical reactivity Thomas et al. (2016), symmetry protected collisions of strongly interacting photons Thompson et al. (2017), the Bose-Einstein condensation Kasprzak et al. (2006) or the room-temperature polariton lasing Kéna-Cohen and Forrest (2010) of exciton-polaritons, and ultra-strong coupling in circuit-QED Murch (2016) to mention a few. Condensed matter systems driven out of equilibrium provide optional possibilities for novel properties, e.g. the creation of Floquet-Bloch states Wang et al. (2013) and Floquet-Weyl semimetals Hübener et al. (2017). 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 Ruggenthaler et al. (2011); Tokatly (2013); Ruggenthaler et al. (2014); Pellegrini et al. (2015); Flick et al. (2015); Ruggenthaler (2015), 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. Tokatly (2013); Ruggenthaler et al. (2014); Flick et al. (2015); Flick (2016). Together with new experiments on chemical systems in optical cavities Hutchison et al. (2012); Thomas et al. (2016); George et al. (2016); Ebbesen (2016), this work now opens up the field of Quantum-electrodynamical (QED) Chemistry and QED Materials Thomas et al. (2016); Flick et al. (2017a); Ruggenthaler et al. (2017). In this new field, so far model Hamiltonian schemes have also been used to successfully describe experiments Galego et al. (2015, 2016), but for an ab-initio description a full real-space picture is necessary.

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 Kohn and Sham (1965), where a whole range of different approximation schemes for the xc functional are available Marques et al. (2012), QEDFT still lacks a practical method to construct such approximations. Previous works Tokatly (2013); Ruggenthaler et al. (2014); Pellegrini et al. (2015) has opened the path to the development of such exchange-correlation functionals. Different routes are possible, e.g. functionals based on e.g. the electron density, the electronic orbitals or the electron current Schäfer et al. (2017) 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. Pellegrini et al. (2015). OEP has been previously used for purely electronic systems in DFT Kümmel and Perdew (2003a); Kümmel and Kronik (2008); Engel and Dreizler (2011). In Ref. Pellegrini et al. (2015), 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 Sternheimer (1951) 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 which 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 Becke (2014) that defines KS-DFT as occupied-orbital-only. Similar schemes have been used in the context of density-functional perturbation theory Baroni et al. (2001) and in many body-perturbation theory using the self-energy approach Giustino et al. (2010).

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 Räsänen et al. (2007); Flick et al. (2015). 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 calculatable, a promising avenue in the design of QED materials is introduced.

## Ii Theory

We start by introducing the general nonrelativistic setup of the correlated electron-photon systems considered in the present work following previous works Tokatly (2013); Ruggenthaler et al. (2014); Pellegrini et al. (2015); Flick (2016). Let us consider interacting electrons of spin or located in an optical cavity. The electrons are coupled to quantized electromagnetic modes, i.e. photon modes. Each photon mode is identified by its cavity frequency and polarization direction . We describe the matter-photon coupling in the Coulomb gauge, dipole approximation (long-wavelength approximation) and the length gauge Craig and Thirunamachandran (1998); Tokatly (2013). The hereby emerging electron-photon coupling strength parameter is projected on the photon polarization direction . 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 . In this setup, the total electron-photon Hamiltonian reads Tokatly (2013); Pellegrini et al. (2015)

(1) | ||||

where each photon mode is associated with a bosonic creation and annihilation operator (, ) that creates and destroys photons in the mode . The transversal electron-photon interaction consists of two terms that read explicitly

(2) |

In dipole approximation, the electron-photon coupling comprises the electron dipole operator and the photon displacement coordinate . The electronic coordinates are defined with respect to the center of charge of the system . As has been outlined in earlier work, using the creation and annihilation operators, we can setup the photon displacement and photon momentum operators and Flick et al. (2015). Physically these two operators are directly connected to the electric displacement field and the magnetic field, if summed over all modes Flick et al. (2015, 2017a). The electron-photon coupling strength is given by

(3) |

where denotes the mode function, e.g. a sine-function for the case of a cubic cavity Ruggenthaler et al. (2014); Pellegrini et al. (2015) and 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 . 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. Flick et al. (2017a), 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, i.e. the electron density and the photon displacement variables . It can be shown Ruggenthaler (2015) that these two internal variables are in an one-to-one correspondence to the the external variables and . Here corresponds to the first-order time-derivative of an external charge current at time zero projected on the mode , i.e., at Tokatly (2013); Ruggenthaler et al. (2014). 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 Tokatly (2013); Ruggenthaler et al. (2014). Since the displacement field corresponds to the electric field minus the polarization Flick et al. (2017a); Rokaj et al. (2017), and the electric field is the time derivative of the vector potential, the conjugate external variable to needs to contain a time-derivative as well. In this work, we choose , without loss of generality. For , or we can find a unitary transformation Dimitrov et al. (2017) that eliminates corresponding terms in Eq. 1. This transformation introduces instead a static external displacement field. 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 Kohn and Sham (1965) of density-functional theory introduced for electron-photon problems in Refs. Ruggenthaler et al. (2014); Tokatly (2013); Flick et al. (2015) 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 Tokatly (2013)

(4) |

for Kohn-Sham orbitals with spin . The effective Kohn-Sham potential is given by

(5) |

and can be divided into the external potential , the usual Hartree-exchange-correlation (Hxc) potential that accounts for the electron-electron interaction and the mode-dependent meanfield-exchange-correlation potential (Mxc) ^{1}^{1}1In 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 exchange-correlation parts that have to be approximated. In exact KS-QEDFT, these parts are chosen such that the electron density that is the sum of the spin-resolved electron densities is equivalent in the interacting and the noninteracting system. In the ground state, we have a simple connection between the exchange-correlation energy

(6) |

that includes contributions from the electron-electron interaction and the electron photon interaction and the corresponding xc potential that reads as follows Engel and Dreizler (2011)

(7) |

This connection will be now exploited to setup the OEP equation. Throughout this work, we use the exchange-only approximation, i.e. . While we use the standard definition for Kümmel and Perdew (2003a); Kümmel and Kronik (2008), i.e. the Fock energy, we focus in the following on the exchange energy due to the electron-photon interaction . 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 Pellegrini et al. (2015), if the photons are not exposed to an external current . Therefore, the leading order becomes the second-order in and the exchange energy can be written as an orbital functional as Pellegrini et al. (2015)

(8) | |||

where refers to the complex conjugate of all former terms. Additionally, we define the projected dipole operator . As does the electron-photon interaction Hamiltonian in Eq. 2, also the electron-photon exchange energy consists of two parts, both containing different electronic orbital shifts. The first orbital shift is the solution of the following Sternheimer equation

(9) | ||||

with the matrix element and the orbital shift can be written explicitly as

(10) |

The second orbital shift is defined by

(11) | ||||

(12) |

While both orbital shifts can be formulated explicitly in terms of all KS orbitals (in Eq. 10 and Eq. 11, respectively), only the second orbital shift can be formulated explicitly in terms of solely occupied orbitals given by Eq. 12. However, the shift 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 , the exchange energy is the Lamb shift of the ground state Pellegrini et al. (2015). Thus all corrections for the ground state are in their magnitude on the order of the Lamb shift. For electron-photon problems, we we find that 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 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

(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 ^{2}^{2}2Please note, that for brevity, we do not explicitly evaluate the contributions, but state its implications if necessary.

(14) | ||||

(15) | ||||

(16) |

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 . In the case of the KS orbitals, this derivative is given by Kümmel and Kronik (2008); Engel and Dreizler (2011)

(17) |

where the sum runs over all orbitals, except . All remaining terms in Eq. 13 are functional derivatives of the orbital shifts. We start by discussing , since it is conceptually simpler to obtain, than . From Eq. 12, for an infinitesimal change in , i.e. , by keeping only first-order terms and combining with Eq. 17, we obtain

(18) | |||

The derivation of the remaining functional derivative of the first orbital shift, i.e. is given in full detail in appendix VI.1 and we only state the final result here

(19) | ||||

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 Kümmel and Kronik (2008)

(21) |

Due to the energy dependence of , we find that the nonvanishing additional inhomogeneity Engel and Dreizler (2011) is given by

(22) | ||||

and the orbital shift by

(23) | ||||

The orbital shift contains in the first line the electron-electron interaction, we choose the exchange-only approximation, i.e.

(24) |

and 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 Pellegrini et al. (2015).

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. Kümmel and Perdew (2003a); Engel and Dreizler (2011) that provide the formulation for electrons-only.

(25) |

and the orbital shifts can be obtained using a Sternheimer equation

(26) |

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 +1 Sternheimer equations that only invoke occupied orbitals.

### ii.1 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 Ruggenthaler (2015), 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, i.e. 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

(27) | ||||

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 . By performing this connection, we assume that the photon number is dominated by contributions stemming 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 Ruggenthaler et al. (2017) are electron-photon correlation functions. For instance, the charge-density-displacement-field correlation function we define as

(28) | ||||

(29) | ||||

where 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 Tokatly (2013); Flick et al. (2015). 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 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 can locally change its sign. Consequently probing this correlation function spectroscopically could allow to obtain novel insights into structural properties of complex systems.

### ii.2 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 Eq. 34 and Eq. 35, it may be hard to converge. One way to circumvents this problem and to obtain an explicit formula for the xc potential is the Krieger-Li-Iafrate (KLI) approximation Krieger et al. (1990, 1992a, 1992b).

In contrast to the common Coulomb OEP equation Kümmel and Perdew (2003a), in the case of correlated electron-photon coupling an additional inhomogeneous contribution appears, i.e. .
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 Kümmel and Perdew (2003a) and decompose Eq. 25 and Eq. 26 starting from

(30) |

with Eq. 23 and

(31) |

we arrive at the exact reformulation

(32) | ||||

The common approximation scheme now assumes which is exact for a single electron if no inhomogeneity was present in Eq. 25. A corresponding substitution involving leads in the general case to nodal points. The variety of possibilities result in different deficiencies and inconsistencies (see also Engel et al. Engel et al. (1998)). To remain as consistent as possible we decide to assume and the KLI equation reads then

(33) | ||||

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.

### ii.3 Numerical details

We have implemented the OEP equation of Eq. 25 and the corresponding KLI equation of Eq. II.2 in the OCTOPUS package Marques et al. (2003); Castro et al. (2006); Andrade et al. (2015). The OEP equation can be solved using standard algorithms as e.g. described in Ref. Kümmel and Perdew (2003a), i.e. in a self-consistent field (SCF) cycle. To obtain the numerical algorithm, we reformulate Eq. 25 as follows

(34) |

The quantity 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

(35) |

Different schemes to calculate are possible Hollins et al. (2012), e.g. dividing by the electron density Kümmel and Perdew (2003b), using the Barzilai-Borwein minimizer Barzilai and Borwein (1988), or connecting to conjugate-gradient algorithms Hollins et al. (2012). However for our purpose, we find that choosing a constant value Kümmel and Perdew (2003a) is already sufficient and already provides the most stable and reliable algorithm. Thus, we choose a.u. for the azulene molecule and a.u. for the sodium chains.

As in the case of electronic OEP Kümmel and Perdew (2003a); Krieger et al. (1992b), 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 Kümmel and Perdew (2003a); Krieger et al. (1992b) and enforce the condition , we find that in the single electron case, the single Kohn-Sham energy deviates from the total energy. From a physical point-of-view it is desirable that both coincides to connect to ionization energies. We find by fixing to the contribution of the highest occupied orbital to the exchange energy defined in Eq. 8, i.e.

(36) | |||

we can restore this connection. However, we note that for the electronic OEP Kümmel and Perdew (2003a); Krieger et al. (1992b) 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.

## Iii 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, i.e. a GaAs quantum ring located in an optical cavity Flick et al. (2015, 2017b), where we have published exact results previously Flick et al. (2015, 2017b). 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 Fig. 1.

### iii.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 () interacting with the single photon mode with frequency meV and polarization direction . The polarization direction enters via the electron-photon coupling strength, i.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

(37) |

with parameters meV, meV, nm, Räsänen et al. (2007); Flick et al. (2017b). For the electron-photon coupling strength, we choose two values, i.e. in weak coupling with meV/nm and in strong coupling meV/nm. The effective three-dimensionality of this problem (two-dimensional electron and one-dimensional photon mode) is low enough that an exact solution is still accessible via exact diagonalization Flick et al. (2014). To obtain the exact ground state, we employ a two-dimensional grid of grid points in each direction with a grid spacing of nm to describe the single electron. The photon field is represented in the photon number eigenbasis and we include up to Fock number states. Using the exact wave function, we can numerically construct the exact exchange-correlation potential Flick et al. (2015); Helbig et al. (2009).

We start by discussing the weak-coupling limit, where meV/nm. In Fig. 2 (a), we show the exact ground-state density obtained by exact diagonalization. Compared to the bare electronic ground-state (for ) that also has a ring structure in the weak-coupling limit, we find small distortions Flick et al. (2017b). In Fig. 2 (b), 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 , i.e. 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 () only in low-density regions, i.e. 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 Fig. 3 we show how the KLI approximation (Sec. II.2) 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 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.

theory | pot | [(meV/nm)] | [meV] | |
---|---|---|---|---|

exact | s | 0.0034 | 33.8782 | 0.0004738 |

OEP | s | 0.0034 | 33.8782 | 0.0004730 |

KLI | s | 0.0034 | 33.8782 | 0.0004727 |

exact | s | 0.1342 | 35.3072 | 3.1926 |

OEP | s | 0.1342 | 35.3349 | 3.4011 |

exact | a | 0.1342 | 32.4816 | 2.2053 |

OEP | a | 0.1342 | 32.4875 | 2.2087 |

To quantify the differences for this example, we print the results of our calculations in Tab. 1. The first three rows show the exact, OEP and KLI results for the total energy and the photon number in the weak-coupling 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 Fig. 4, we show the results obtained in the strong-coupling regime ( meV/nm), where we find a deviation in the exact ground-state electron density from the ring structure in the weak-coupling regime to a double-well structure Flick et al. (2017a) as shown in Fig. 4 (a). This splitting is accompanied by a higher peak in the xc potential in the center of the grid as shown in Fig. 4 (d). Although in the weak-coupling limit, we find a very high accuracy of the OEP approximation, in the strong coupling limit, we observe the break-down of the OEP approximation. In Fig. 4 (b), 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 Fig. 4 (e). Consequently the error of the OEP density and the potential shown in Fig. 4 (c) and (f) are very high.

The origins of this failure of the OEP can be understood by calculating the photon number and the double occupancy in the photon mode shown in Fig. 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 Flick et al. (2017b). 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. Flick et al. (2017b), 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 Dimitrov et al. (2016)). 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, i.e. single photon processes, to higher order processes, i.e. two-photon, three-photon, etc. is required to accurately describe this limit.

In the last example, we study an asymmetric example in the strong-coupling limit ( meV/nm), where the external potential is given by

(38) |

with meV/nm. The cavity frequency is again chosen to be in resonance to the first electronic excitation, i.e. meV. The results are shown in Fig. 6. We find while the density is approximated accurately with errors in the order of , also observables such as the photon number listed in Tab. 1 are approximated quite accurately due to the dominant mean-field contribution in Eq. 27.

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

### iii.2 Azulene (CH) molecule in an optical cavity

theory | 24-1 [eV] | 25-24 [eV] | [eV] | [eV] | [eV] |
---|---|---|---|---|---|

KLI | 16.57 | 2.24 | -1648.39 | -501.79 | 0.00 |

OEP | 16.68 | 2.42 | -1648.53 | -503.04 | 0.00 |

KLI-PT | 16.48 | 2.25 | -1644.38 | -502.11 | 3.89 |

OEP-PT | 16.66 | 2.54 | -1644.71 | -503.67 | 3.79 |

Our next example is the first real application of the QEDFT framework to a three-dimensional molecule, i.e the azulene (CH) molecule. 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 Kim et al. (2015) (CID: 9231). Second, we use the geometry optimization of the OCTOPUS package employing the LDA functional Kohn and Sham (1965); Perdew and Zunger (1981) to calculate a relaxed ground-state structure. Third, using this relaxed structure, we use the electron OEP to calculate a HOMO-LUMO gap of eV that serves as the cavity frequency, i.e. eV. The electron-photon coupling includes the polarization direction of the photon field that is polarized along the x-direction with a strength of , i.e. ^{3}^{3}3For standard experimental parameters, e.g. for a single trapped-atom cavity as described in Ref. Khitrova et al. (2006) (Fig. 3), (V = m), we deduce an experimental value of .. In this example we want to investigate the question how the correlated electron-photon interaction alters the electronic ground-state density . To numerically calculate the ground-state density of the azulene molecule, we use a grid of dimensions Bohr in directions. The grid spacing is chosen to be Å and to describe the core electrons of the carbon and hydrogen atoms we use LDA Troullier-Martins pseudopotentials Troullier and Martins (1991). 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 Kümmel and Perdew (2003a) also in the OEP setting.

In Fig. 8, we show in the top panel the ground-state 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 Fig. 8 shows the difference of the electronic ground-state density exposed to electron-photon coupling to the bare electronic ground-state density, i.e. 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 Fig. 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 , it correctly predicts the density accumulation regions at the top and bottom of the molecule. However, these regions are overestimated by a factor . To quantify the effects of the quantized electron-photon interaction on many-electron systems, we have provided numerical results in Tab. 2. For different levels of theory, we print the energy difference between lowest and highest occupied orbitals (), the HOMO-LUMO gap (), the total energy , and the electronic and the electron-photonic part of the exchange energy , and , respectively. For the given parameters, the electron-photon exchange energy is in the order of eV and two orders of magnitude smaller than the electronic exchange energy that is roughly eV. As could be expected, the changes due to the coupling to the vacuum of the cavity are small in the ground state, i.e., 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 as defined in Eq. 28 without the mean-field contributions in Fig. 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-photon interaction on molecules in cavities. We find that in QEDFT new observables become numerically accessible that could allow for novel experimental spectroscopic setups Ruggenthaler et al. (2017).

### iii.3 Chain of Sodium Dimers

The last example that is studied in this paper is a chain of sodium dimers of variable length, i.e. up to ten sodium dimers. We use this set up 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 Garraway (2011), 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, e.g., the superradiant phase transition Knight et al. (1978); Vukics et al. (2014). With a first-principles 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. Kümmel et al. (2000). For the optical cavity frequency, we choose the energy of the 3s-3p transition, i.e. eV. We assume that the photon field is polarized along the direction of the sodium chains with a strength of , i.e. . To calculate the chain of sodium dimers (Na), we use the sodium pseudopotentials and equilibrium distances for a single sodium dimer of Ref. Kümmel et al. (2000). For the real-space grid, we use dimensions Bohr with grid spacing Bohr, where is the chain length. The distance between the sodium dimers is chosen as Bohr.

The case of three sodium dimers is illustrated in Fig. 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 Fig. 10 shows a cut along the y-axis in blue against the ground-state electron density in the cavity in grey. 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 in-between the dimers, where the amount of density is decreasing in the cavity.

The next figure, Fig. 11 shows the case of ten 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 Fig. 12. Here we find the KLI strongly overestimates the effects of the electron-photon interaction.

In the last figure of this section, Fig. 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 . Thus, we can capture this behavior alternatively by defining a new effective coupling constant 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 emitters.

## Iv Summary and Conclusion

In conclusion, in this work, we have illustrated how the cavity-mediated 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 first principles, 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 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 e.g. Rabi splittings, etc. can be substantial. The present work constitutes the first mandatory step towards such studies of excited states of strong light-matter coupled quantum systems. Additionally, this approach could also benefit standard electronic DFT, since similar ideas, i.e. 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, where chemical systems are placed in optical cavities.

## V Acknowledgements

We would like to thank Claudiu Genes, Camilla Pellegrini, and Ilya V. Tokatly for insightful discussions and acknowledge financial support from the European Research Council (ERC-2015-AdG-694097), by the European Union’s H2020 program under GA no.676580 (NOMAD), and the Austrian Science Fund (FWF P25739-N27).

## Vi Appendix

### vi.1 Derivation of the functional derivative of second orbital shift

The derivation of the functional derivative of the second orbital shift with respect to the Kohn-Sham potential can be derived analogously to the derivation of Eq. 17 discussed in Ref. Engel and Dreizler (2011). By keeping the first order terms, we find the following Sternheimer equation that defines the infinitesimal change in the orbital shift

(39) |

For , we can employ the following relation Engel and Dreizler (2011)

(40) |

The Sternheimer equation can be solved explicitly and has the solution