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).1-5 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, for example, electron-photon correlation functions, become accessible in a density-functional theory (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.

Over 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−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 quantum-electrodynamical (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 Ne = ∑σ=↑,↓Nσ = N↑ + N↓ interacting electrons of spin ↑ or ↓ located in an optical cavity. The electrons are coupled to Np 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 – r|. 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 (âα†, âα) that creates and destroys photons in the mode α. The transversal electron–photon interaction Ĥep(α) consists of two terms that read explicitly

In dipole approximation, the electron–photon coupling comprises the electron dipole operator R = R0 + ∑r and the photon displacement coordinate . The electronic coordinates r are defined with respect to the center of charge of the system R0. 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 bywhere Sα denotes the mode function, for example, a sine-function 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 via the external potential vext(r). The effect of a static permanent dipole moment due to the nuclear charge can be neglected, since the two terms of eq 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 one-to-one correspondence to the external variables vext(r) and jext(α). Here jext(α) corresponds to the first-order time-derivative of an external charge current at time zero projected on the mode α, that is, ∫d3rSα(kα·r)eα·∂jext(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 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 jextα ≠ 0 into one with jextα = 0, we can construct by the very same simple transformation the solutions for the inhomogeneous case from the case with jextα = 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 φ(r) with spin σ. The effective Kohn–Sham potential v(r) is given byand can be divided into the external potential vext(r), the usual Hartree-exchange-correlation (Hxc) potential vHxcσ(r) that accounts for the electron–electron interaction and the mode-dependent meanfield-exchange-correlation potential (Mxc) vMxcσ(α)(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 exchange-correlation 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) = ∑φ*(r)φ(r) is equivalent in the interacting and the noninteracting system. In the ground state, we have a simple connection between the exchange-correlation energythat 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, Exc ≈ E(ee) + ∑α=1E(α). While we use the standard definition for E(ee),[28,29] that is, the Fock energy, we focus in the following on the exchange energy due to the electron–photon interaction E(α). The interaction terms in eq 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 jext(α). Therefore, the leading order becomes the second-order in λα. While the second part of eq (the dipole self-interaction 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 , also the electron–photon exchange energy E(α) consists of two parts, both containing different electronic orbital shifts. The first orbital shift is the solution of the following Sternheimer equationwith the matrix element d(α) = ⟨φ|d̂α|φ⟩ 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 and 11, respectively), only the second orbital shift Φ(2) can be formulated explicitly in terms of solely occupied orbitals given by eq . However, the shift Φ(1) can be defined implicitly by a Sternheimer equation that only invokes occupied orbitals as given in eq .

Since the exchange energy given in eq 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(α), as defined by eq , 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 in eq , we now proceed to calculate the corresponding Kohn–Sham potential using functional derivatives. From eq , we can setup the following OEP equation by using the chain rule of functional derivatives

The OEP equation of eq 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 and are given as follows (please note that, for brevity, we do not explicitly evaluate the E(ee) contributions, but state its implications if necessary)

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 vs. In the case of the KS orbitals, this derivative is given by[29,30]where the sum runs over all orbitals, except i = j. All remaining terms in eq are functional derivatives of the orbital shifts. We start by discussing Φ(2)(r), since it is conceptually simpler to obtain, than Φ(1)(r). From eq , for an infinitesimal change in Φ(2)(r), that is, δΦ(2)(r), by keeping only first-order terms and combining with eq , we obtain

The derivation of the remaining functional derivative of the first orbital shift, that is, δΦ(1)(r)/δvs(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 an alternative OEP equation can be derived that is given by the simple equationwith the Kohn–Sham Green’s function[29]

Due to the energy dependence of E(α), we find that the nonvanishing additional inhomogeneity[30] Λ(r) is given byand the orbital shift M(r) by

The orbital shift M(r) contains in the first line the electron–electron interaction, we choose the exchange-only approximation, that is,and E( 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]

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.and the orbital shifts ψ*(r) can be obtained using a Sternheimer equation

This equation has to be solved self-consistently with eq . By this procedure, we have replaced the problem of calculating the OEP equation using all unoccupied states by a problem of solving Np + 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 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 is due to the quantum fluctuations of the photon mode, the latter term is a classical contribution due to a nonvanishing R0. 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 aswhere Ψ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 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, Λ. 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 with the Kohn–Sham potential[28] and decompose eqs and 26 starting fromwith eq andwe arrive at the exact reformulation

The common approximation scheme now assumes (ĥ(r) – ε)ψ*(r) = 0, which is exact for a single electron if no inhomogeneity would be present in eq . A corresponding substitution involving ψ*(r) ≈ Λ(r)/φ(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 (ĥ(r) – ε)ψ*(r) = 0 and the KLI equation reads then

This reformulation avoids the solution of eq 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 and the corresponding KLI equation of eq 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 self-consistent field (SCF) cycle. To obtain the numerical algorithm, we reformulate eq as follows

The quantity Sσ(r) becomes a measure for convergence, since it is vanishing in the case of convergence (compare eq and eq ). 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 , 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 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 .

Figure 1

Overview of the three molecular systems in an optical cavity studied in the present work: (I) GaAs quantum ring, (II) single azulene molecule, (III) chain of ten Na2 dimers all of which are coupled to a single photon mode with frequency ωα and electron–photon coupling strength λα.

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 = re + re) 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 formulawith parameters ℏω0 = 10 meV, V0 = 200 meV, d = 10 nm, and m0 = 0.067m.[35,49] For the electron–photon coupling strength, we choose two values, that is, in weak coupling with λα = 0.0034 meV1/2/nm and in strong coupling λα = 0.1342 meV1/2/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.[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 meV1/2/nm. In Figure a, 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 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 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 () 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 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 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 . The first three rows show the exact, OEP and KLI results for the total energy Etot and the photon number npt in the weak-coupling limit using the external potential of eq . 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 , we show the results obtained in the strong-coupling regime (λα = 0.1342 meV1/2/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[19] as shown in Figure a. This splitting is accompanied by a higher peak in the xc potential in the center of the grid as shown in Figure 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 breakdown of the OEP approximation. In Figure 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 Figure e. Consequently, the error of the OEP density and the potential shown in Figure c,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 Figure . 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 self-consistency 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, two-photon, 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 meV1/2/nm), where the external potential is given bywith 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 . 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 are approximated quite accurately, since in this regime the mean-field contribution in eq becomes dominant in comparison to the fluctuations.

Figure 2

Exact (a) and OEP approximated (b) electron density in the weak-coupling limit (λα = 0.0034 meV1/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.

Figure 3

Exact (a) and KLI approximated (b) electron density in the weak-coupling limit (λα = 0.0034 meV1/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.

Table 1

Results of the Self-Consistent KS Calculation for the GaAs Quantum Ring in an Optical Cavity: The Total Energy Etot and Photon Number npt for Different Levels of Theory, Coupling Strength, and Symmetric (s) and Asymmetric (a) Potentials (pot)

theory	pot	λα (meV1/2/nm)	Etot (meV)	npt
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

Figure 4

Exact (a) and OEP approximated (b) electron density in the strong-coupling limit (λα = 0.1342 meV1/2/nm). The difference is shown in (c). The corresponding xc potentials are shown in (d) and (e), respectively, (f) shows the difference in the xc potentials.

Figure 5

Plot of the photon number occupation ⟨â†â⟩ and double photon number ⟨â†â†ââ⟩ for the GaAs quantum ring as a function of the electron-photon coupling strength λα. The inset shows the region λα ∈ [0.06, 0.09] meV1/2/nm. When the double occupancy becomes significant, the OEP approximation starts to fail (see text for more detail).

Figure 6

Exact (a) and OEP approximated (b) electron density in the strong-coupling limit (λα = 0.1342 meV1/2/nm). The difference is shown in (c). The corresponding xc potentials are shown in (d) and (e), respectively, (f) shows the difference in the xc potentials.

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 (C10H8) 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 (C10H8) molecule (Figure ). 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 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 eV1/2/nm (0.08 au), that is, λα = 37.47 eV1/2/nm e. [For standard experimental parameters, e.g., for a single trapped-atom cavity, as described in ref (56) (Figure , V = 18.148 μm3), we deduce an experimental value of λ0 = 3.9 × 10–7 eV1/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 n0(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 , 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 Figure 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 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 . 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 Etot, and the electronic and the electron–photonic part of the exchange energy E(ee) and E(α), 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 Eee 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 , without the mean-field contributions in Figure . 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.[20]

Figure 7

Azulene (C10H8) molecule in an optical cavity, λα denotes the polarization direction of the photon field.

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.

Table 2

Results of Self-Consistent KS Calculation for Real 3D Azulene in an Optical Cavity: Energy Difference between the Highest Occupied Orbital (HOMO; 24th Orbital) and the Lowest Occupied Orbital (1st orbital), Energy Difference between the Lowest Unoccupied Orbital (LUMO; 25th orbital) and the Highest Occupied Orbital (HOMO), the Total Energy Etot, the Exchange Energy E(ee), and the Photon Exchange Energy E(α) for Different Levels of Theory

theory	24–1 (eV)	25–24 (eV)	Etot (eV)	Ex(ee) (eV)	Ex(α) (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

Figure 9

Correlation function as a cut in x–z plane A(α)(r) as defined in eq , calculated for the azulene molecule.

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 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 (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 eV1/2/nm (0.006 au), that is, λα = 2.95 eV1/2/nm e. To calculate the chain of sodium dimers (Na2), 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, 2Nc × 10) × 60 Bohr with grid spacing 0.5 Bohr, where Nc 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 . 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 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 in-between the dimers, where the amount of density is decreasing in the cavity.

Figure 10

From top to bottom as a cut in x–y plane: OEP ground-state 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.

The next figure, Figure 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 . Here we find the KLI strongly overestimates the effects of the electron-photon interaction. In the last figure of this section, Figure . We plot the number of photons in the correlated electron-photon ground state using the functional presented in eq . In total, we find for the KLI and the OEP a linear behavior, where the KLI overestimates the number of photons slightly. From eq we also see that ⟨âα†âα⟩ ∼ λα2. 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 N emitters.

Figure 11

From top to bottom as a cut in x–y plane: OEP ground-state 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 ground-state 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.

Figure 13

Photon occupation ⟨âα†âα⟩ for the case of variable chain length of sodium dimers.

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