Molecular orbital theory in cavity QED environments

Coupling between molecules and vacuum photon fields inside an optical cavity has proven to be an effective way to engineer molecular properties, in particular reactivity. To ease the rationalization of cavity induced effects we introduce an ab initio method leading to the first fully consistent molecular orbital theory for quantum electrodynamics environments. Our framework is non-perturbative and explains modifications of the electronic structure due to the interaction with the photon field. In this work, we show that the newly developed orbital theory can be used to predict cavity induced modifications of molecular reactivity and pinpoint classes of systems with significant cavity effects. We also investigate electronic cavity-induced modifications of reaction mechanisms in vibrational strong coupling regimes.

T he use of strong light-matter coupling to modify molecular properties and reactivity is nowadays a very popular topic in physics and chemistry [1][2][3][4] . Many groundbreaking works have shown that interaction with confined fields can impact many matter features ranging from the modification of absorption and emission spectra [5][6][7][8] to the alteration of photochemical processes 9,10 . In particular, the chemistry community has focused on how and to which extent reactions can be engineered by coupling with light [11][12][13][14] . Interaction with quantum fields can indeed significantly affect molecular processes both in the ground state and excited states 15,16 with reported examples of reactions slowing down 14 , speeding-up 17 or becoming selective towards one product 18 . The easiest way to achieve strong coupling between light and matter is through optical cavities (See Fig. 1a), where the frequency of the electromagnetic field is determined by the geometrical features of the device 19 . Inside the cavity, the photonic vacuum couples to the molecular system creating polaritonic states 20 with distinct features 21 . Most importantly, the properties of the mixed matter-photon states can be engineered tuning the photonic part of the system, which means that polaritons represent a very effective way to modulate matter properties in a non-invasive way 22 .
A detailed theoretical description of the strong coupling regime is urgently needed to develop an intuitive picture of cavity chemistry, that would also greatly ease the experimental design. However, this is a challenging task because in the strong coupling regime photons become a critical component of the quantum system and they must therefore be treated as quantum particles following quantum electrodynamics (QED). Only recently, a wider interest from the chemistry community in electron-photon systems has led to the introduction of semi empirical methods 23 , variational theories 24 as well as many others 25 . Only four ab initio methods have been proposed so far: QED Hartree Fock (QED-HF) 26 , quantum electrodynamical density functional theory (QEDFT) [27][28][29] , QED coupled cluster (QED-CC) 26,[30][31][32] and QED full configuration interactions (QED-FCI) 26,33 . Molecular orbital theory for the strong coupling regime has so far not been proposed, although all the methods mentioned above use an orbital basis to parametrize the wave function 26,34,35 . The molecular orbital (MO) concept is a powerful theoretical tool used to  Fig. 1 Orbital properties in optical cavities. a Pictorial representation of an optical cavity with an injection input, the coated glasses, and a spacer to regulate the distance between the mirrors. The frequency of the cavity fields is proportional to the inverse of the distance between the two mirrors. The lambda parameter quantifies the strength of the light-matter coupling. The cavity picture is based on a video realized by the Ebbesen group. b Origin invariance of the SC-QED-HF orbitals compared with the origin dependence of QED-HF orbitals. If the methoxy ion is displaced in space, both the orbital energies and the orbital shapes changes for QED-HF orbitals, while the SC-QED-HF orbitals remain unaltered. c Orbital modifications due to changes in the cavity parameters. Changing the cavity parameters (frequency and coupling) the electronic ground state and therefore the molecular orbitals are changed. Particularly, we show that an avoided crossing like situation can be observed among orbitals.
develop correlated theories and to provide a qualitative and simple interpretation of molecular properties 36 . In particular, since MOs also display local properties of the system like the electrophilicity or nucleophylicity of atomic fragments, reactivity can in many cases be more easily rationalized in terms of orbitals, rather than with electron densities. The introduction of a MO theory for systems in QED environments may therefore significantly enhance our chemical intuition about experiments in strong coupling regimes.
In standard quantum chemistry, molecular orbitals are obtained as solutions to the Hartree Fock (HF) equations, however a proper extension of this concept to QED environments has proven to be non-trivial. Indeed, the orbitals provided by QED-HF display unphysical behaviour for charged systems with both the orbital energies and the shapes changing when the molecule is translated (See Fig. 1b). The orbital construction scheme is therefore not reliable with critical quantities like the HOMO-LUMO gap 37 being ill-defined. Moreover, standard perturbation theory to describe electron-electron and electron-photon correlation cannot be used without a proper orbital theory 26 . In this work, we introduce a new ab initio electron-photon framework called strong coupling quantum electrodynamics Hartree Fock (SC-QED-HF) which naturally leads to the first fully consistent MO theory in QED environments. Our analysis reveals that in order to obtain a well behaved molecular orbital description the photonic field contribution must be accounted for in a nonperturbative way, ensuring that in the limit of infinite coupling the exact wave function is obtained. The newly developed method captures electron-photon correlation allowing us to display field induced effects on the electronic ground state while keeping the computational cost lower than post HF approaches.

Results
The photon-matter interaction is described through the well known dipole Hamiltonian [38][39][40] where H e is the electronic Born-Oppenheimer Hamiltonian, d is the molecular dipole operator, λ ¼ ffiffiffiffiffi ffi is the coupling constant between the molecule and the field, where V is the quantization volume. The ω and ϵ are the frequency and polarization of the electric field, and the bosonic operators b and b † respectively annihilate and create photons. Besides H e , the dipole Hamiltonian is composed of three main terms: the purely photonic term, an interaction term and the dipole self energy λ 2 2 d Á ϵ ð Þ 2 . The self energy is often neglected, however this term is critical to ensure the Hamiltonian is bounded from below and displays the correct scaling with the size of the system 41 . To describe a molecule in a cavity, we must solve the eigenvalue problem for the dipole Hamiltonian where ψ is the electron-photon wave function. Within the SC-QED-HF approach the wave function is approximated as Here a y pσ and a pσ respectively create and annihilate an electron in orbital p with spin σ, HF j i denotes an electronic Slater determinant and j0 ph i is the photonic vacuum. The electron-photon correlation basis defining the orbitals p is obtained by diagonalizing the operator (d ⋅ ϵ). The η p parameters in Equation (3) are orbital specific coherent state coefficients as explained in more details in the Supplementary Methods. Using SC-QED-HF we obtain origin invariant and size-extensive energies, and most importantly, we can construct fully origin invariant molecular orbitals, see Fig. 1b. Non size-extensive effects in the cavity can be modeled by optimizing the correlation basis and including more photons in the reference wave function in Equation (3). This modification could be crucial if long range collective effects need to be modeled. Related approaches have been proposed by Mandal et al. 25 , Ashida et al. 39 and Schäfer et al. 42 . These methods also use a unitary electron-photon transformation to partially decouple the electrons and photons. In particular, the approaches in Refs. 39 and 42 also start from the solution of the infinite coupling case and we therefore expect the same eigenvalues in this limit. The method in Ref. 25 , on the other hand, is based on a CI like scheme where the dipole operator is diagonalized in a limited configuration space. However, developing an MO theory based on this method is not trivial.
Importantly, since the unitary SC-QED-HF transformation in Equation (3) is variationally optimized, the description of the system for finite coupling is dramatically improved (See Supplementary Fig. 1). We emphasize that SC-QED-HF can also be used to study any kind of strongly coupled systems such as plasmonic nanocavities 43,44 and superconducting circuits 45 .
As displayed in Fig. 1c, both frequency and coupling variations can induce significant modifications of the MOs. Both λ and ω can be varied by changing the geometrical features of the cavityopening the possibility of engineering ground state molecular properties by tuning the field parameters. We point out that while the coupling is intrinsically a light-matter property, the frequency is a field property only. The frequency dependent orbitals of SC-QED-HF are dressed with the photon field according tõ where U is the unitary transformation from the HF basis to the correlated basis. Both occupied and unoccupied orbitals are affected by the cavity parameters, however, larger changes are typically observed for the latter. Molecular orbitals are used to qualitatively predict and interpret molecular excitations as well as bond formations, therefore the cavity-induced modification of the orbital character reported in Fig. 1c suggests significant changes in the molecular properties. Orbital avoided crossings and orbital mixing are other effects that can be induced by coupling to the photon field. In particular, when two orbitals have the same symmetry, an "avoided crossing" like situation can be observed (the orbitals mix but their energies never become degenerate, see Fig. 1c). Conversely, when two almost degenerate orbitals have different symmetry, a crossing is allowed and the orbital ordering is simply exchanged, see Supplementary Fig. 2. Whenever the HOMO-LUMO gap is small we expect to see particularly large cavity effects and to achieve significant control over molecular properties. In these cases the near degeneracy between HOMO and LUMO can be enhanced or lifted by the cavity fields 26 with considerable impact on the ground state wave function. A reliable scheme to model these kinds of systems might require combining the SC-QED-HF Hamiltonian with an active space method 46 .
One of the main applications of molecular orbital theory is the rationalization of molecular reactivity. Specifically, MOs are used to predict the critical steps of a reaction, as well as the effects that a modification of the conditions can induce on the final outcome.
Here we investigate whether the molecular orbitals provided by SC-QED-HF can be used to predict cavity effects on the formation of a complex from its initial constituents. As a case study, we chose a prototypical example of a nucleophilic substitution (SN2): a methyl chloride reacting with an ammonia molecule to produce methylamine and hydrogen chloride. A particularly interesting configuration along the reactive path is the one reported in Fig. 2 where both the C-Cl and the C-N bonds are not fully formed nor broken. Since this arrangement is similar to the transition state of the SN2 reaction, its cavity induced stabilizations or destabilizations play a critical role in determining a reaction speedup or slowdown. Analyzing the MOs of the complex we notice that there is only one occupied orbital with a significant bonding component over both the chloride and the nitrogen, HOMO-5. The HOMO-5 orbital of the complex is composed of the nitrogen lone pair, the C-Cl sigma bond and the three C-H sigma bonds. Projection on the orbitals of the fragments allows us to monitor field induced changes in the character of HOMO-5. Specifically, in Fig. 2 we display how HOMO-5 and its projections change for three different values of λ: 0.0, 0.1 and 0.2. Coupling values of 0.1 or 0.2 are extremely high, far superior to what is nowadays experimentally achieved for molecules interacting with a single mode of the quantum field. However, larger values are used here to approximate multimode and collective effects 33 (See methods). As the coupling increases, HOMO-5 becomes more localized on the two different fragments. Specifically, the sigma bond between the nitrogen and the carbon disappears as observed from the orbitals of the fragments. Indeed their overlap slightly decreases with the coupling. We notice in particular that for λ = 0.2, the HOMO-5 orbital has the character of a sigma bond between the carbon and the chloride plus a non-bonding lone pair from the nitrogen. The orbital analysis therefore suggests that cavity induced effects make the configuration in Fig. 2 less favourable inside the cavity than outside. This statement is confirmed using an energy argument with the formation energy of the complex (ΔE(λ) = E complex (λ) − E fragments (λ)) increasing by 5.30 kJ/mol for λ = 0.1 and 17.36 kJ/mol for λ = 0.2, compared to outside the cavity.
In panels b and c of Fig. 2, we display the difference between the ground state electron densities computed with and without the cavity for different couplings (Δρ = ρ λ − ρ nocav , λ = 0.1, 0.2).
In particular, the isosurface shows changes of the electronic density due to the cavity. To further analyze the charge readjustment we report an integration of the density difference along the field polarization in Fig. 2 26,33 . We observe that as the coupling increases there is a significant charge depletion among both the C-Cl and the C-N bonds, as expected from the orbital analysis. This behaviour confirms the conjecture presented in Refs. 26,29,41 that upon interaction with the field, electrons reorganize to minimize the electronic dipole component along ϵ. We note that the electronic density redistribution observed in Fig. 2   In dark blue we plot the integral of the density differences to further display cavity induced charge reorganization 26,33 . The polarization of the field is placed along the C-Cl bond and the cavity frequency is 13.6 eV.
not induced by resonance effects between a molecular excited state and the cavity field. It is due to the interaction with the quantum field only. This shows that even in the non-resonant case, modifications of molecular properties can be expected. This charge reorganization might actually play a role in determining cavity induced modifications of molecular reactivity. This idea was also reported by Galego et al. 16 and tested using a transition state method for the strong light-matter coupling case. Strong coupling between molecular vibrations and cavity fields, vibrational strong coupling (VSC), has proven to be a particularly effective way to affect the outcome of molecular reactions 14,17,18,47 . Although the number of experiments displaying to what extent reactions can be engineered through VSC is constantly increasing, several aspects still remain unclear. Indeed, while the formation of vibrational polaritons certainly contributes to determine the observed effects, questions still remain on the role played by the electrons and the solvent 14,48 . Moreover, it has been suggested that an important part might be played by collective effects 49 and by dark states 50 . Just recently, multiple works have suggested that the reaction slowdown in the VSC regimes can be described through nuclear effects only. Explanations use either the dynamical caging of the photon fields 51 or cavity induced energy redistribution among normal modes 52 . However, in both cases electron-photon interaction was approximated preventing an accurate quantification of the cavity induced electronic effects on the chemical reactivity. Since SC-QED-HF captures some electron-photon correlation we will use it to assess the relevance of the electronic changes in VSC. We will focus on electronic effects and disregard the coupling between the fields and the nuclear vibrations. In particular, this means that the formation of the vibrational polaritons is not considered in the following analysis. We compute the potential energy surfaces (PES) at the HF and SC-QED-HF level of theory for two different reactions that are known for respectively being slowed down and catalyzed under vibrational strong coupling conditions: the deprotection of 1-phenyl-2-trimethylsilylacetylene (PTA) with tetra-n-butylammonium fluoride (TBAF) 14 (Fig. 3a) and the solvolysis of para-nitrophenyl acetate (PNPA) with TBAF 17 (Fig. 3b). The field polarization was aligned along the Si-C and the C-O bonds, respectively, while the frequencies of the cavity field are chosen to match the ones reported in Refs. 14 and 17 . A first inspection of the potential energy surfaces (See Supplementary Figs. 3, 4) reveals that for both reactions the overall shape is not changed by the cavity. Specifically, in the case of PTA-deprotection an intermediate is formed before the Si-C bond is broken (second step in Fig. 3a) and a transition state barrier must be overcome. Conversely, in the case of the PNPAsolvolysis the reaction is concerted. In Fig. 3 we show the difference between the SC-QED-HF and HF potential energy surfaces. The colour maps in Fig. 3c, d visualize the QED effects highlighting which configurations are stabilized by the cavity and therefore whether we should expect a reaction slowdown or a speedup. In particular, the dark zones are stabilized and the bright zones being destabilized upon coupling to the photons. In both the analyzed cases the reagents are more stabilized than the products pointing towards an overall slowdown of both the reactions. In the PTA-deprotection case our prediction is in agreement with what is experimentally observed, while in the case of PNPA, electronic effects would predict a slowdown instead of a speedup. Most importantly, however, for the two reactions in Fig. 3 the QED modifications of the electronic structure are so small they do not seem to affect in a significant way the shape of the PES. Stabilizations reach a maximum of 0.01 eV in the case of the PTA-deprotection and to about 0.07 eV in the PNANsolvolysis case. Moreover, we point out that rate changes in the experiments are only observed if resonance between the cavity mode and the nuclear vibrations is achieved. Conversely, similar cavity induced modifications of the PESs are observed in a wide range of cavity frequencies (top-right corner of Fig. 3c, d). In conclusion the observed cavity induced electronic effects are not large enough to justify the experimental outcome. In particular we point out that no resonance condition is found in the electronic structure properties. We therefore conjecture that electronic effects alone are of secondary importance in the vibrational strong coupling regime and that no significant variations of the reaction mechanism or of the transition state should be expected for the reported examples. This is in disagreement with what was reported by Galego et al. Although different molecular systems were studied in Ref. 12 , the discrepancy with our findings is due to the inclusion of the dipole self energy contribution in the Hamiltonian in Equation (1) 41 . Summarizing, our studies show that electronic effects can be relevant under the appropriate experimental conditions (see Fig. 2, near degeneracy systems) but do not seem to play an important role in resonant VSC experiments. This is reasonable considering that nuclear properties are only indirectly related to electronic effects.

Discussion
In this work we have introduced a non-perturbative electronic structure method for the strong coupling regime, strong coupling QED-HF. The approach effectively captures part of the electronphoton correlation and scales as N 3 using an optimal implementation. Moreover, it leads to the first fully consistent molecular orbital theory for QED environments. The method is completely general and can therefore be applied to any kind of quantized field environment 44 . The obtained orbitals have the correct fermionic/bosonic statistics and are dressed with the photonic field as shown in Equation (4). Furthermore, in the infinite coupling limit the SC-QED-HF wave function is an exact eigenfunction of the light-matter Hamiltonian. Similarly to standard molecular orbitals, the dressed orbitals can be used to rationalize molecular reactivity. From the considered examples we confirm that electrons interacting with the quantized field get more localized and, in this way reduce the electronic dipole in the electric field direction. The importance of cavity induced electronic effects in vibrational strong coupling can also be investigated using the new method. Our study shows that QED effects on the electronic ground state can be significant under the proper experimental settings (see Fig. 2), however they only play a secondary role for the analyzed VSC cases. Therefore a cavity induced change of the transition states in Fig. 3 is unlikely; in accordance with the conclusions reported in other works 51,53 . Furthermore, SC-QED-HF opens the way for developing perturbative methods to capture additional electron-electron and electron-photon correlation with a lower computational cost than, for instance, QED-CCSD-1. This should in particular allow us to describe non size-extensive effects in cavities that are critical to reproduce strong coupling modifications of intermolecular interactions 33 . Finally, the introduction of a molecular orbital theory opens to the possibility of developing multi-level methods to simulate strong light-matter coupling for large molecular systems and include collective effects in an ab initio framework.

Methods
The SC-QED-HF approach is a fully ab initio electron-photon correlated framework. The method has been implemented in a development version of the eT program 54 . The optimization of the wave function is performed using the gradients of the energy with respect to the η p parameters and the density matrix. The gradients are computed analytically and their explicit expressions are reported in the Supplementary Methods. The DIIS extrapolation procedure is adopted to speed up convergence 46 . The electronic density is initialized as a superposition of atomic densities while the starting guess for the η p parameters is chosen to equal the eigenvalues of d ⋅ ϵ. We have tested other starting values of the η p parameters, however the final result is not affected by the initial guess. SC-QED-HF allows us to incorporate some electron-photon correlation in the wave function without invoking the complexity of methods like coupled cluster or full configuration interaction. In particular, the current implementation of SC-QED-HF scales as N 5 due to the calculation of the two electron integrals in the dipole basis. Future optimization of the method, including screening of the two electron integrals, will reduce the scaling to N 3 . Additional theoretical details regarding SC-QED-HF are reported in the Supplementary Methods. All geometries shown in this paper have been optimized using the ORCA package 55 . We used DFT-B3LYP/def2-SVP basis set to optimize the geometries for the small molecules (methoxyde, methanol) as well as for the structure in Fig. 2. For the optimization of the complex in Fig. 2 the distance between the nitrogen and the carbon was fixed to 2.1 Å while the positions of all the other atoms were reoptimized. The fragments in Fig. 2 were obtained from the geometry of the complex without any additional optimization. In this way we isolated cavity induced effects neglecting modifications of the orbitals induced by the geometry relaxation. Finally, the complex formation energy ΔE was computed from the energies of the fragments at infinite separation. The geometries for the PNPA potential energy surfaces were optimized using DFT-PBE0/def2-SVP basis set while the PTA geometries were computed using DFT-B3LYP/6-31G* basis set similarly to what has already been reported in literature 52 . All the calculations reported in this paper have been performed using a cc-pVDZ basis set.
New experimental techniques keep pushing the boundaries of the ultrastrong light-matter coupling regime by extreme confinement of the field quantization volume, using for instance nanoplasmonic picocavities 43 , where quantization volumes can be reduced down to the nm 3 limit. However, some of the calculations presented in this work have been performed for values of the light-matter coupling that can seem unrealistically high, with effective quantization volumes of around 0.2 nm 3 . This choice is justified because inside optical cavities molecules interact with more than one mode at a time and field induced effects from different modes sum with each other. The net effect can be approximated using a larger value of a single mode coupling to the molecule, λ eff . In particular, in the QED-HF case, the multimode calculation is identical to using a single mode with an effective coupling, λ eff ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi N modes p λ 12,33 . This relation does not hold exactly for SC-QED-HF because of the electron-photon correlation, however the net effect of a multimode approach still remains an increase of the light-matter coupling.

Data availability
The data that support the findings of this study are available from the corresponding author upon request.

Code availability
The custom code used for this study is available from the corresponding author upon request. A publicly available version of the code will be released following the time schedule of the eT program "https://etprogram.org/".  In the top corners we show that the energy for the PES configurations is monotone with respect to the cavity frequency and therefore no resonance effects are observed. The coupling parameters are λ = 0.025 and λ = 0.044 respectively, while the frequencies are in resonance with the experimental Si-C and the C=O vibrational frequencies. In the insets of panel c and d the resonance features of SC-QED-HF are investigated. The yellow star indicates the resonant frequency. The QED effects have been shifted so that at resonance they are equal to zero.