Coupled Cluster Theory for Molecular Polaritons: Changing Ground and Excited States

We present an ab initio correlated approach to study molecules that interact strongly with quantum fields in an optical cavity. Quantum electrodynamics coupled cluster theory provides a non-perturbative description of cavity-induced effects in ground and excited states. Using this theory, we show how quantum fields can be used to manipulate charge transfer and photochemical properties of molecules. We propose a strategy to lift electronic degeneracies and induce modifications in the ground state potential energy surface close to a conical intersection.


INTRODUCTION
Nowadays, manipulation by strong electron-photon coupling is becoming a popular technique to design and explore new states of matter [1]. Recent advances in experimental and theoretical research include new ways to generate exciton-polariton condensates [2], induce phase transitions [3][4][5], tune exciton energies in monolayers of 2D materials and interfaces [2,6,7], and even enhance the electron-phonon coupling with possible effects on superconductivity [8,9]. Nevertheless, general techniques for manipulating molecules via strong coupling have not yet reached maturity.
Chemistry is one of the fields that has witnessed most progress in strong light-matter coupling applications. In particular, Ebbesen and coworkers have found that strong coupling to vibrational excited states in molecules can inhibit [10][11][12], catalyze [13,14], and induce selective change in the reactive path of a chemical reaction [15]. These experiments use an optical cavity, the simplest device where entanglement between matter and light can be observed. In an optical cavity (see Fig. 1), the quantized electromagnetic field interacts with the molecular system, producing new hybrid light-matter states called polaritons [16]. These states exhibit new and interesting properties, leading to unexpected phenomena. In the past few years, improvements of optical cavities [17][18][19][20] have resulted in devices that can reach the strong and ultra-strong coupling limits. This has made polaritonic states accessible at room temperature, also for a small number of molecules [18,21].
Theoretical modeling is an essential tool to provide fundamental understanding and outline new strategies * henrik.koch@sns.it for applications in polaritonic chemistry. The challenge is to develop an accurate theoretical description of entangled light-matter systems. In the quantum optics community, several groups have developed model Hamiltonians to reproduce the main features of polaritonic physics [12,[22][23][24][25][26]. The objective of the current work is to formulate and implement a quantitative ab initio method for polaritonic chemistry.
Presently, the only available ab initio theory is quantum electrodynamical density functional theory (QEDFT) [1,[27][28][29], which can describe interacting electrons and photons, on an equal footing. This method is a natural extension of density functional theory (DFT) [30] to quantum electrodynamics (QED). In QEDFT, the Kohn-Sham formalism treats electrons and photons as independent particles interacting through an exchange correlation potential. The QEDFT method is computationally cheap and reproduces the main polaritonic features for large systems, even at the mean-field level. However, its accuracy is limited due to the unknown form of the exchange correlation functional for the electron-photon interaction. In particular, in a mean-field treatment, no explicit electron-photon correlation is accounted for in the ground state. The problem can be overcome with a properly designed exchange-correlation functional, but current functionals are still not sufficiently accurate. Recently, Rubio and coworkers [31,32] proposed an extension of optimized effective potentials (OEPs). However, the accuracy of this functional still needs to be assessed.
Coupled cluster theory is one of the most successful methods for treating electron correlation in both ground and excited states of molecular systems [33,34]. It is nowadays routinely applied to compelling chemical problems due to major advances in computational resources and several decades of algorithmic developments. Coupled cluster methods are available in several programs; we use the electronic structure program e T [35] for our developments. In this paper, we extend coupled cluster theory to treat strongly interacting electron-photon systems in a non-relativistic QED framework. We refer to the resulting method as quantum electrodynamics coupled cluster (QED-CC) theory. To the best of our knowledge, this is the first coupled cluster formalism that incorporates many-body electron-photon operators for an ab initio Hamiltonian. Recently, a different coupled cluster formalism was proposed by Mordovina et al. [36]. Their study was limited to model Hamiltonians and they used state-transfer operators, instead of many-body operators, to describe the photonic part of the wave function. We should also mention that studies using other electronic structure methods and model or semi-empirical Hamiltonians have been presented by other authors [37][38][39].
In addition to presenting the complete formulation and implementation of QED-CC, we consider some interesting applications in photochemistry. In particular, we demonstrate how light-induced charge transfer in small dye molecules, commonly used as prototypes for photovoltaic applications, can be modified by the quantized electromagnetic field. Furthermore, we show that the presence of the cavity can break molecular symmetry and change relaxation mechanisms. Suitably defined fields can induce significant changes in both ground and excited state properties. These results pave the way for novel strategies to control photochemical reaction paths.

II. COUPLED CLUSTER THEORY FOR ELECTRONS
In this section we introduce notation and important concepts needed to develop the electron-photon interaction model. For a complete outline of coupled cluster theory, we refer to Ref. [34]. In standard coupled cluster theory for singlet states, the many-body wave function is expressed using the exponential parametrization, where |HF is a reference wave function that is usually chosen to be the closed-shell Hartree-Fock (HF) determinant. The cluster operator T generates electronic excitations when operating on the reference and, in this way, the exponential produces a superposition of Slater determinants. In the case of purely electronic many-body states, the cluster operator is defined as where N e is the number of electrons. Each term corresponds to excitations (single, double, triple, and so on), i.e.
where E ai = a † aα a iα + a † aβ a iβ (here, a and a † denote fermionic operators and (α, β) denote spin projections) , are singlet excitation operators and the parameters t ai and t aibj are called cluster amplitudes. Furthermore, we let indices (i, j, k, l) and (a, b, c, d) label occupied and virtual HF orbitals, respectively. General orbitals are labeled (p, q, r, s). The cluster operator can be expressed as where the excitation operators τ µ generate an orthonormal set of excited configurations: Together with |HF , these configurations define a subspace of the Hilbert space in which the Schrödinger equation is solved. Inserting the coupled cluster wave function in Eq. (1) into the time-independent Schrödinger equation, we obtain where H e is the electronic Born-Oppenheimer Hamiltonian The quantities h pq and g pqrs are one-and two-electron integrals respectively and, for convenience, we have introduced the operator In coupled cluster theory, Eq. (7) is projected onto the set {|HF , |µ }. Consequently, the coupled cluster energy is given by and the cluster amplitudes are determined from the equations Ω µ = µ|H e |HF = 0.
Here we have introduced the similarity transformed HamiltonianH e = exp(−T )H e exp(T ).
Coupled cluster theory is equivalent to exact diagonalization of the Hamiltonian in Eq. (8), also called full configuration interaction (FCI), when the cluster operator in Eq. (2) is untruncated and contains all possible excitations. When the excitation space is truncated, we obtain different levels of approximation and, needless to say, reduced computational cost. For example, T = T 1 defines the coupled cluster singles model (CCS) and T = T 1 + T 2 the coupled cluster singles and doubles model (CCSD). Coupled cluster theory is manifestly size-extensive, also in its truncated forms, a property that ensures that the total energy of non-interacting subsystems is the sum of the subsystem energies. This is unlike similar truncation in configuration interaction theory, where extensivity errors can become arbitrary large for an increasing number of subsystems.
Another important feature of coupled cluster theory is the size-intensivity of excitation energies: for non-interacting subsystems, excitation energies in each subsystem do not change with the total system size [40]. The excitation energies are the eigenvalues of the Jacobian matrix Since this matrix is non-Hermitian, special attention is required at electronic degeneracies; at such points, the matrix can become defective/non-diagonalizable. Defects are therefore expected close to conical intersections, as discussed in more detail in Section III C. In coupled cluster theory there are two prevailing approaches to electronic excited states. One is coupled cluster response theory (CCRT), which is based on a time-dependent formalism [40]. This theory provides both size-extensive and size-intensive molecular properties, such as excitation energies and transition moments. The other is based on a time-independent formalism and is referred to as equation of motion coupled cluster (EOM-CC) theory [41]. In EOM-CC theory, the excitation energies are the same as in CCRT. However, some molecular properties are not guaranteed to scale correctly with system size. For instance, transition moments are not necessarily size-intensive [42]. For the purpose of the present developments, which mainly relates to ground and excited state energies, the EOM-CC formalism is sufficient. In EOM-CC, the similarity transformed Hamiltonian is expressed in the basis {|HF , |µ }, The left and right eigenvectors ofH define the excited state vectors, L k | and |R k , and the eigenvalues ofH are the energies of the states. We have used Eqs. (10), (11), and (13) in the last equality. To extend coupled cluster theory to electron-photon systems, we need to introduce a new parametrization of the cluster operator.

III. COUPLED CLUSTER THEORY FOR ELECTRON-PHOTON SYSTEMS
To describe the interaction of the electromagnetic field with atoms, molecules, and condensed matter systems, the low-energy limit of QED is usually sufficient [43,44]. In particular, the non-relativistic Pauli-Fierz Hamiltonian in the dipole approximation [1,43,45,46], is usually an accurate starting point for describing electronic systems in optical cavities. The Hamiltonian in Eq. (15) is expressed in Coulomb and length gauge [43][44][45]. We use the Born-Oppenheimer approximation and keep the nuclear positions fixed. The second term in the Hamiltonian is the purely photonic part, represented by a sum of harmonic oscillators, one for each frequency. We have neglected the zero-point energies. The operators b † α and b α are bosonic creation and annihilation operators, respectively. The third term is the dipole self-energy term, which ensures that the Hamiltonian is bounded from below [46] and independent of origin. The last term, the bilinear coupling, couples the electronic and photonic degrees of freedom. In the length gauge, the light-matter coupling is via the dipole operator which consists of an electronic and a constant nuclear contribution. In the electronic term, d pq denotes oneelectron dipole integrals. The coupling is described through the transversal polarization vector e multiplied by the coupling strength λ α : Here, ε 0 and ε r are the permittivities of the vacuum and the dielectric materials separating the cavity mirrors, respectively. The α mode quantization volume is denoted V α . The dipole approximation is usually valid when the wavelength of the electromagnetic field is significantly larger than the size of the electronic system. There are, however, cases where the dipole approximation is not sufficiently accurate. For instance, this occurs when the size of the system is comparable to the cavity wavelength, or when the matter interacts with a circularly or elliptically polarized field. These aspects will be discussed in a forthcoming paper.

A. The QED-HF method
In order to formulate QED-CC, we need to define a suitable reference wave function. We formulate an extension of the Hartree-Fock method to QED, hereafter referred to as QED-HF. The non-correlated electrons and photons in QED-HF are described by where |0 is the photon vacuum and c n are expansion coefficients for the photon number states. The coefficient c n , where n = (n 1 , n 2 , . . . ), corresponds to the state with n α photons in mode α. Starting from Eq. (18), and assuming that |R is normalized, the energy, is minimized with respect to Hartree-Fock orbitals and photon coefficients c n . Note that QED-HF can be considered a special case of QEDFT with an appropriately chosen exchange-correlation functional. For a given Hartree-Fock state, the energy can be minimized with respect to the photon coefficients. This can be achieved by diagonalizing the photonic Hamiltonian where the mean value is with respect to |HF . This Hamiltonian can be diagonalized by the unitary coherent state transformation [47] U with a suitable choice of z. In the transformed basis, the photonic Hamiltonian is If we choose the Hamiltonian reduces to The eigenvectors of this operator are the photon number states, and the lowest eigenvalue corresponds to the vacuum state. In the untransformed basis, that is, for the Hamiltonian in Eq. (21), the eigenstates are the generalized coherent states where z α is given in Eq. (24) and are the normalized photon number states for mode α.
Applying the unitary transformation to the original Pauli-Fierz Hamiltonian, Eq. (15), gives us the final expression for the Hamiltonian in the coherent state basis: This Hamiltonian will be used in QED-CC. Note that the operator is manifestly origin invariant, differently from Eq. (15), where the invariance is obtained through a gauge transformation. This is also true for charged systems, where the dipole moment operator depends on the choice of origin. In the coherent state basis, the bilinear coupling and self-energy terms depend on fluctuations of the dipole moment away from the mean value. As we have shown, the QED-HF reference state is now given by The Hartree-Fock equations are solved using standard techniques [34,48]. In every iteration, the Hartree-Fock orbitals are updated and used to evaluate z 0 , see Eq. (24). The inactive Fock matrix [34] used in the optimization is given by where F e is the inactive electronic Fock matrix [34] and the stationary condition is equivalent to F ia = 0. The QED-HF ground state energy can now be written as where the correction to the electronic energy can be understood as the variance in the dipole interacting with the photon field.
We should point out that the eigenvalues of F, normally interpreted as orbital energies, are origin dependent. As a consequence, applying concepts that depend on the orbital energies, like the Koopmans' theorem [49], will require a different choice of the occupied-occupied and virtual-virtual blocks of F. For the same reason, F cannot be used as a zeroth order Hamiltonian in perturbation theories such as CC2 [50] and CC3 [51,52]. However, the origin dependence of the eigenvalues of F does not imply a loss of origin invariance in QED-HF (see Appendix A).

B. The QED-CC method
Extending the exponential parametrization in Eq. (1) to QED requires that the cluster operator generates excitations both in the purely photonic and electron-photon coupling spaces [36]. The cluster operator can therefore be partitioned as where T e is the standard cluster operator for the electrons and T p and T int consist of photon and electron-photon operators, respectively. The purely photonic operator is defined as In this equation, γ n are photon amplitudes. The form of the photonic operator was chosen to expand the photonic part of the Hilbert space and to give commuting cluster operators (as in the electronic case). Since exp(T p ) does not terminate, the parametrization is able to incorporate many-body effects within the limitation imposed by the projection space. In contrast, Mordovina et al. [36] used a nilpotent photonic operator that only enters linearly in the expansion of the coupled cluster state. The excitations in the electron-photon interaction operator T int are defined as direct products of electronic and photonic excitations. Thus, the operator can be expressed as where, for instance, The new cluster amplitudes, γ n , s n ai , s n aibj , etc., are parameters that will be determined from a set of projection equations.
Hierarchies of approximations are formulated by truncating the cluster operator and the associated projection space. Here we implement the special case of a single photon mode, where we only include one photon in the cluster operator. The coupled cluster wave function in QED-CC is given by where |R is the QED-HF reference given in Eq. (29), and the cluster operator is Electronic excitations are described at the singles and doubles level, and the photon mode is coupled to these excitations through S 1 1 and S 1 2 , respectively. The electronic operators T 1 and T 2 are given in Eqs. (3) and (4). The photon and electron-photon operators are defined as This model is referred to as QED-CCSD-1 with one photon mode, where "1" refers to the photonic excitation order. More involved terminology will be required to describe the full hierarchy. In the notation used by Mordovina et al. [36], this is a QED-CC-SD-S-DT model. However, due to the difference in photonic excitation operators and the coherent state basis, the model described here is not directly comparable to the one in Ref. [36]. Even if only one photon creation operator is included, the exponential will partially incorporate two photon contributions into the wave function. Thus we expect the convergence with respect to photons will be faster than CIlike diagonalization in photon number states. Furthermore, notice that the generalized coherent states basis incorporates higher photonic excitations as well.
The projection space used in Eq. (11) is defined from the excitations included in the cluster operator. With the notation |HF, n = |HF ⊗ |n (42) |µ, n = |µ ⊗ |n , the projection basis is where |HF, 0 = |R and µ is restricted to single and double excitations. The derivation of the amplitude equations follows the same procedure as in the electronic case, and the truncation of the equations is determined by the projection space and the commutator expansion of the similarity transformed Hamiltonian. Explicit formulas are presented in Appendix B.
The formation of polaritons usually appears in the optical spectrum as a Rabi splitting (proportional to the coupling strength λ) of the electronic states due to the coupling to the quantum field. Hence, we must also describe the excited states of the coupled system. In coupled cluster theory, electronic excitation energies may be determined using EOM-CC theory, as described in Section II. The projection space in QED-CCSD is extended, relative to the electronic case, giving rise to additional blocks in the Jacobian matrix: In addition to the electronic Jacobian A e,e , there are coupling blocks between electronic (e), electronic-photonic (ep), and photonic (p) configurations; see Eqs. (38) and (44). Explicit formulas for the sub-blocks of the Jacobian, are presented in Appendix C, along with the corresponding sub-blocks of the η vector. The ground and excited state QED-CCSD equations are solved using standard methods [34,53,54]. Additional properties of the electron-photon system can be calculated from the left and right eigenvectors of H. For instance, we can evaluate the ground and excited state electronic EOM density matrices as For a description of other molecular properties, we refer the reader to the literature [40,41].

C. Some technical aspects
The Pauli-Fierz Hamiltonian in Eq. (15) is defined on the direct product Hilbert space H = H e ⊗ H p . In the truncated description, where H e and H p are finitedimensional, exp(T )|R has an effectively finite expansion due to the finite projection basis. For instance, in QED-CCSD-1, the terms in exp(T )|R that give non-zero contributions are up to quadruple electronic excitations and double photonic excitations. However, in the limit of an infinite-dimensional H p , special care might be required to define the exponential operator [55].
Due to the non-Hermiticity of coupled cluster theory, it is known to give non-physical complex energies close to conical intersections between excited states of the same symmetry [56][57][58]. The same issues can arise in QED-CC and were mentioned by Mordovina et al. [36]. The problem can be traced to defects in the Jacobian matrix for a truncated projection basis. In the untruncated case, the states satisfy generalized orthogonality relations, ensuring a correct description of conical intersections [58]. To obtain a physically correct description with a truncated excitation space, one can enforce the orthogonality conditions where P is some projection operator [59,60]. This approach, unlike a posteriori corrections [36,57], can be extended to analytical energy gradients and nona- diabatic coupling elements using well-established Lagrangian techniques [61,62]. In passing, we note that no defects or complex eigenvalues were encountered in the results reported in this work.

IV. MOLECULAR POLARITONS
In this section, the QED-CCSD-1 model is used to investigate cavity-induced effects on the chemistry of molecules. All calculations are performed using a development version of e T [35]. Molecular geometries are provided in Supplemental Material [63]. A single cavity mode is used throughout; this approximation typically breaks down for small values of ω cav , when several replica states overlap energetically with electronic states. This usually happens in large cavities.

A. Diatomic molecules
Interesting QED effects can be observed for small diatomic molecules, such as H 2 and HF. We also use these molecules to benchmark the coupled cluster model against the more accurate QED-FCI approach. The comparison shows an excellent agreement, see Appendix E for a detailed discussion. For the calculations here we use a Gaussian basis set, in particular, cc-pVDZ [64].
The potential energy curves for the ground and excited states of these systems are shown in Fig. 2. The conical intersections and avoided crossings in the UV range define the photochemical properties of these molecules. An optical cavity set in resonance with one of the excited states can completely restructure the excitation landscape and redefine the photochemistry of the system. The color map in Fig. 2 indicates the electronic/photonic contributions to the states. The electronic states are highlighted in blue, while the photonic states are transparent white. For more details, see Appendix D.
Considering first H 2 set in resonance with the first singlet excited state (at the ground state equilibrium geometry), we are able to induce significant changes of the potential energy curves. Here we focus on the main Rabi splitting around −0.6 a.u., where in particular, the upper (UP) polariton is more bound than the bare electronic state. Hence, it should be possible to trap the molecule in the UP state. In contrast, the bare electronic state has, to a larger extent, a dissociative character.
Similar conclusions can be drawn for the HF molecule. Differently from H 2 , it has a permanent dipole moment. However, this seems to have minor effects on the general landscape. This is consistent with the fact that the total energy depends only on the fluctuations in the dipole moment (see Section III A).
From the above analysis, we see that the application of a quantum field inside an optical cavity can be used to fine-tune the excited state properties of molecules. This opens the way towards new decay paths, the possibility of trapping systems in excited dark polaritons, and many other effects re-designing the molecular photochemistry.

B. Charge transfer molecules
Recently, some research groups have suggested that quantum fields can have a significant impact on the charge transfer [65,66] and energy transfer [65,67] properties in molecules. These preliminary studies are based on model Hamiltonians; thus, only qualitative interpretations of the phenomena are provided. Here we present a quantitative analysis of cavity-induced effects on a charge transfer process.
We investigate p-nitroaniline (PNA), a simple amine often used as a prototype dye for solar energy applications, for instance in dye-sensitized solar cells [68,69]. This molecule has an intense low-lying charge transfer excitation (at about 3-4 eV) that potentially can be used to inject charge in a semiconductor and produce a current. Developing strategies to control the charge transfer process is of fundamental importance to increase photovoltaic efficiencies.
In our calculations on PNA, we have used the cc-pVDZ basis and oriented the polarization (λ = 0.05) along the principal axis of the molecule. The molecular structure was optimized with DFT/B3LYP using the 6-31+G* basis set [70].
Initially we investigate the effects of the cavity on the electronic ground state by analyzing the electron density. This can conveniently be carried out using the charge displacement analysis [71][72][73]. The charge displacement function is defined as where we integrate over the electron density difference  has been moved along the z coordinate. In particular, if ∆q(z) is positive, charge is transferred from right to left, and if negative, charge is transferred in the opposite direction.
In Fig. 3, we show the charge displacement function and isosurface for the ground state electron density difference with and without cavity, ∆ρ = ρ cav gs −ρ nocav gs . The cavity is set in resonance with the most intense low lying charge transfer excited state, ω cav = 0.178 a.u. Although the charge displacement is small, a clear cavity-induced charge reorganization is observed in the ground state. Specifically, we have a charge transfer of about 0.005 e − going from the acceptor (NO 2 ) to the nitrogen atom of the donor (NH 2 ) group. This counter-intuitive effect is in agreement with previous QEDFT/OEP studies from Flick et al. [32]. The cavity field accumulates more charge in the high-density regions. In this way, the variance of the dipole operator is reduced, see Eq. (31).
In the excited states of PNA, the cavity-induced effects are more evident. In Fig. 4 we show the dispersion of the low-lying excitation energies of PNA with respect to the cavity frequency ω cav . Due to the large transition dipole moment for the charge transfer excitation, a large Rabi splitting is observed when the cavity is resonant with this state. As discussed in the previous section, the cavity can induce significant changes in the excited states. Specifically, we have a state inversion between the lower polariton and the first two excited states. This effect in PNA and other dye molecules could be important for photovoltaic applications, where a proper alignment of charge transfer states with the states in a semiconductor is essential to optimize solar cell efficiency.
In Fig. 5, we show a charge displacement analysis of the charge transfer state, with a resonant quantum field. Here we use the density difference between the ground and excited state, ∆ρ = ρ es − ρ gs , with and without the cavity. Both the lower and upper polariton have charge transfer character and are shown separately.
Differently from the ground state, now a sizable charge transfer of almost 0.4 e − are moved from the donor to the acceptor group. The cavity divides the total charge transfer between the polaritons; thus, a compromise must be made between energetically aligning states and maintaining the charge transfer character. We note that the cavity field slightly reduces the total charge transfer (see the black dashed line in Fig. 5). This is because the charge transfer state also contributes to the other excited states, not just the polaritons.
In Fig. 6 we show how fine-tuning the cavity frequency can be used to change the degree of charge transfer in the lower and upper polaritons. This opens another possibility for charge transfer control, with potential for technological applications.

C. Photochemical processes
We now turn our attention to photochemical processes and the possibility of changing the ground state potential energy surface using an optical cavity. For this purpose, we choose the pyrrole molecule that exhibits conical intersections between the ground state and two low lying excited states. A detailed analysis of these conical intersections, and of the involved relaxation mechanism, can be found in Ref. [74,75]. In C 2v symmetry, the ground state is 1 A 1 while the first two excited states are 1 A 2 and 1 B 1 . The equilibrium geometry is calculated with CCSD and cc-pVDZ basis set.
We investigate the behavior of the potential energy curves when the NH bond distance R is varied, preserving the C 2v symmetry (Fig. 7). The polarization of the coupling strength is set to λ = 0.05 and the cavity frequency is set in resonance with the 1 B 1 state at R = 2.0 A(ω cav = 0.039 a.u).
In Fig. 8, we show the potential energy curves along the coordinate R with and without the cavity. Without the cavity (λ = 0), the CCSD model mostly reproduces the accurate potential energy curves calculated in Refs. [74,75]. The main difference is an inverted ordering of the 1 A 2 and 1 B 1 states close to the conical intersection. The correct ordering can be recovered by including triple excitations in the model, as we have confirmed with CC3 [52] calculations (see Appendix F). Since the ordering of the states does not change the conclusions, we performed the analysis at the QED-CCSD-1 level.
We first observe the lifting of the degeneracy as shown in Fig. 8d. In the coupled system, now in C 1 symmetry, all states can interact as they have the same symmetry. This unequivocally demonstrates that the cavity can also significantly impact the potential energy surface of the electronic ground state. This is the first time this phenomenon has been demonstrated in a molecule using an ab initio Hamiltonian. Previously, a similar observation was made with a model Hamiltonian for graphene [3].
Interestingly, a coupling strength (λ = 0.05) that produces a relatively small Rabi splitting of 0.001 a.u. around 2Å is able to open a gap twice this size, 0.002 a.u., at the conical intersection. We can rationalize this observation in the following way. The Rabi splitting is mainly due to the bilinear term in the Hamiltonian, whereas the lifting of the degeneracy is mainly due to the self-energy term. In Appendix F we show that the symmetry breaking is insensitive to basis sets and cavity frequencies.
A consequence of the above observation is the possibility to change relaxation pathways in chemical reactions. For instance, the relaxation through conical intersections to the ground state can be reduced, such that the relaxation can be dominated by radiative processes (which are generally slower). Note also that the majority of molecular orientations allow for this type of symmetry breaking. This should make the gap opening experimentally observable in standard temperature conditions where the molecules can rotate freely.

V. CONCLUDING REMARKS
We have developed a coupled cluster theory that explicitly incorporates quantized electromagnetic fields, denoted as QED-CC. This non-perturbative theory can describe molecular photochemistry inside an optical cavity. The QED-CC model is a natural extension of the well established coupled cluster model, used in electronic structure theory. The method provides a highly accurate description of electron-electron and electron-photon correlation, at least in regions where the electronic ground state is dominated by a single determinant. These correlations are not accounted for in commonly used model Hamiltonians and mean field methods. The accuracy is demonstrated by comparison with exact diagonalization (within an orbital basis) in a truncated photon space (QED-FCI). Unlike QED-FCI, the QED-CC hierarchy is computationally feasible also for larger molecules.
Initially, we investigated the restructuring of potential energy curves in diatomic molecules. In particular, we found that the interaction with the cavity creates polaritons which are more bound than the corresponding bare electronic excited states. Clearly, polaritons have crucial implications on the photochemistry. For instance, they can alter the relaxation pathways and trap molecules in dark excited states. A further study of these phenomena would be very interesting.
Cavity-induced effects on charge transfer processes were also investigated quantitatively for PNA. We explained how the cavity restructures the charge inside the molecule, and how these effects could be applied in photovoltaics.
Finally, we demonstrated how the cavity field can be used to manipulate conical intersections in molecules. We showed that the quantum field is able to lift degeneracies between ground and excited state. This analysis suggests that new experimental strategies can be developed to manipulate ultrafast molecular relaxation mechanisms through conical intersections.
When averaging over the photon vacuum state |0 and using Eq. (9) we obtain This operator has the same form as the electronic Hamiltonian in Eq. (8), with modified integrals and constants. These modified integrals can be inserted directly into the expression for the inactive electronic Fock matrix, and we obtain Eq. (30) for a single mode. This procedure is easily generalized to the multimode case. We now consider the origin invariance in QED-HF by shifting the dipole by a constant, d → d + ∆d. For neutral molecules the dipole moment operator is origin invariant, ∆d = 0, but for charged molecules this is not the case. Nevertheless, since the Hamiltonian (A1) is invariant, the energy is also invariant. On the other hand, the inactive Fock matrix in Eq. (30) is not invariant. As seen from the shifted Fock matrix, the occupied-virtual block of the Fock matrix, F ib and F aj , are unchanged, whereas the purely occupied and virtual blocks, F ij and F ab , are dependent on the origin. The electronic and total projection operators are here defined as Excited states weights, w e es , are calculated in an approximate way using the norm of the electronic part of the excitation vector R, In principle an equation equivalent to Eq. (D1) should be used also in this case, substituting |CC with |R . Considering that we are only interested in a qualitative estimate of the weights, the approximate form Eq. (D4) is adequate.

Appendix E: Comparison of QED-CCSD-1 and QED-FCI
We here compare QED-CCSD-1 and QED-FCI. For the latter method, we performed an exact diagonalization of the Hamiltonian in Eq. (15) with one photonic excitation, in order to obtain a consistent comparison.
We used a 3-21G basis [76] for H 2 and a STO-3G [77] basis for HF, with internuclear distances R H2 = 1.0Å and R HF = 0.917Å. In both cases, the coupling value λ = 0.05 is used.
In Fig. 9, we show the energy dispersion with respect to the cavity frequency ω cav . In this figure we use the QED-FCI linear response spectral function, where |Ψ n are the eigenfunctions of the QED Hamiltonian. In this case, we calculated the density-density spectral function instead of the more appropriate optical spectrum (with transition dipole moments) in order to compare the states independently from the selection rules. Coupled cluster results are displayed in Fig. 9 as red lines with electronic weights.
For both molecules we observe an excellent agreement. The only noticeable difference is the absence of a few electronically excited states in the FCI spectrum. These states have large double excitation character, and thus, small contributions to the spectral function. Notice here that QED-CCSD-1 is very accurate also for HF, where the electronic structure is not exact in CCSD. We also provide some additional results that are important to validate the accuracy of QED-CCSD-1 for pyrrole. Firstly, we address the basis set appropriateness, as electromagnetic fields may require more diffuse basis functions. In Fig. 10 potential energy curves calculated using an aug-cc-pVDZ [64] basis are shown. The inclusion of diffuse functions does not change the general qualitative picture described in Section IV C, confirming the accuracy of our predictions. In particular, the position of the intersections is nearly unaffected by the size of the basis and the qualitative shape of the potential energy curves is unchanged.
In Section IV C, we noted that CCSD gives and incorrect ordering of the 1 A 2 and 1 B 1 excited states. This can be rectified by including triple excitations in the electronic treatment, as shown with CC3 in Fig. 11.
The opening of the conical intersection is quite robust with respect to changes in the cavity frequency, ω cav , as seen in Fig. 12. This supports the claim that the dipole self-energy term is mainly responsible for lifting the degeneracy. The position of the Rabi splitting is, as expected, highly sensitive to the cavity frequency.