Polaritonic Coupled-Cluster Theory

We develop coupled-cluster theory for systems of electrons strongly coupled to photons, providing a promising theoretical tool in polaritonic chemistry with a perspective of application to all types of fermion-boson coupled systems. We show benchmark results for model molecular Hamiltonians coupled to cavity photons. By comparing to full configuration interaction results for various ground-state properties and optical spectra, we demonstrate that our method captures all key features present in the exact reference, including Rabi splittings and multi-photon processes. Further, a path on how to incorporate our bosonic extension of coupled-cluster theory into existing quantum chemistry programs is given.


I. INTRODUCTION
In recent years seminal experiments at the interface between quantum optics, quantum chemistry, and material sciences have shown that when photons and matter couple strongly the emergence of light-matter hybrid states, called polaritons, can substantially change chemical and physical properties of molecular systems [1][2][3][4][5][6][7][8][9]. This has led to the observation of changes in chemical reactions [10,11], suppression of photon degradation and photobleaching [12,13], tunable third-harmonic generation from polaritons [14], roomtemperature superfluidity in a polariton condensate [15], or modifications of intersystem crossings [16]. While experimentally the influence of strong coupling on matter, e.g., due to placing molecules inside a cavity or plasmonic nanostructure, has been firmly established, theoretical approaches to describe situations where photonic, electronic, and nuclear degrees of freedom become strongly mixed so far do not provide a detailed and general explanation of the observed effects [17][18][19][20].
When reliability and accuracy is concerned, coupledcluster (CC) theory [46][47][48][49] has become the method of choice in quantum chemistry. This wave-function method provides a hierarchy of approximations with truncation based on excitation order from a mean-field reference state. Size-consistent and size-extensive molecular electronic energies, as well as other ground and excited-state properties, can be calculated with chemical accuracy [50]. This makes an extension of CC theory to the case of strong matter-photon coupling highly desirable. Fortunately, there is nothing intrinsic in CC theory that imposes a restriction to purely electronic problems. The requirements for a computationally tractable CC theory can be satisfied for bosonic degrees of freedom as well, as demonstrated by pioneering applications to molecular vibrations [51].
In this article we develop the theoretical framework for CC calculations on systems of interacting electrons coupled to photons, and illustrate the potential of such an approach to describe the behavior of matter strongly coupled to cavity modes with a computational cost that scales polynomially with problem size.

II. METHOD
We consider fixed-nucleus molecular Hamiltonians coupled to a cavity in the dipole limit and in the length gauge [52][53][54] H =Ĥ e + α ω c,αâ † αâ α + γ α ω c,αd (â † α +â α ) + γ 2 α ω c,αd 2 , whereĤ e is the electronic Hamiltonian,â ( †) α are the creation and annihilation operators for a cavity mode with frequency ω c, α , andd is the electronic dipole operator. The coupling parameter γ α tunes the strength of the light-matter interaction; here we focus mostly on cases where γ α > 0.05, where the system is typically considered to be in the strong-coupling regime [19,31]. The theory can be extended to more general light-matter Hamiltonians [17,29,52,53] in a straight-forward manner.
CC theory is based on an exponential ansatz for the groundstate wave function where | 0 is an uncorrelated reference state (usually Hartree-Fock) andT is the cluster operator. The cluster operator is a weighted sum of excitation operatorŝ where μ labels a general excitation in the system, and the t μ amplitudes are to be determined. The exponential form makes the CC state multiplicatively separable, introducing the fundamental feature of size extensivity. It also has the property that even when the cluster operator is truncated to include only low-order excitations, higher-order effects are incorporated through the expansion of the exponential ofT . Standard CC theories are classified by the number of electrons excited in the list of operators τ μ . For example, the method that includes all single and double excitations is called CC singles and doubles (CCSD). The excitation operators are written in terms of fermionic creation and annihilation operatorsĉ ( †) , e.g., the operatorτ a i =ĉ † aĉ i excites an electron from occupied orbital i to unoccupied orbital a.
The amplitudes and the ground-state energy are obtained through projected equations where | μ =τ μ | 0 andĤ is the similarity-transformed Hamiltonian e −TĤ eT . The polynomial scaling of the theory is attributed to the fact that usually a polynomial number of excitation operators entersT and that the Baker-Campbell-Hausdorff expansion ofĤ truncates [48]. The latter condition is met when excitation operators are commutative and nilpotent. Hence, there is nothing special about electrons in the formulation and CC theory can be applied to general many-particle systems. The natural excitation operator for a bosonic mode is simply the creation operatorâ † , but the lack of nilpotency [(â † ) 2 = 0] means that a formalism based onâ † would lack the simple structure of electronic CC theory. However, if the number of photons in the system is limited to n max , which is necessary for numerical treatment of bosons in the basis of Fock number states (see Appendix B for convergence tests), it is possible to map the bosonic mode to a lattice of n max + 1 sites, |0 , |1 , . . . , |n max , each corresponding to a number state; then the excitation operatorŝ FIG. 1. Two types of excitation operators for a bosonic degree of freedom. Both can be used to access any number state, and both are sets of commutative operators. But the operatorsτ n are additionally nilpotent, simplifying CC formulations that use this form.
clearly fulfill both the commutativity and nilpotency condition, while allowing any number state to be addressed (Fig. 1). Now, to build CC theory for electron-photon systems, we introduce a more general cluster operator (6) in which the electronic termsτ μ are supplemented by purely photonic excitationsτ n and connected light-matter excitationŝ τμτñ.
To describe the truncation of each term in the cluster operator, we extend the terminology common in electronic CC theory: in the acronym CC-X -Y -Z, X , Y , and Z will specify the level of electronic, photonic, and mixed excitations respectively; CC-SD-S-0 refers to conventional electronic CCSD with additional single-photon excitations, CC-SD-S-D includes coupled excitations of one electron together with one photonic excitation additionally to CC-SD-S-0, and CC-SD-S-DT adds coupled double electronic together with one photonic excitations to CC-SD-S-D. In the one-mode calculations presented here the photonic excitation level is at most singles, but multiphoton excitations within that single mode are included [see Eq. (5)].
As the reference state we take a product of an electronic Slater determinant and the vacuum state for radiation modes In this form, polaritonic CC theory displays some similarities with the polaron ansatz in quantum optics [55,56]. However, as opposed to CC theory, the polaron ansatz is developed for two-level systems and does not target the electronic structure explicitly.
Additionally, since our method is in principle formulated in the full Hilbert space of the problem, it does not suffer per se from gauge-ambiguity issues exhibited by some effective quantum optical models [57,58].
Because the structure of polaritonic CC theory is closely analogous to that of conventional electronic CC theory, it is possible to anticipate many of the advantages of the formalism. The properties of the operators in Eq. (5) ensure truncation of the energy and amplitude equations to produce a polynomially scaling theory. Because the photon mode is modeled by a lattice, the situation is identical to one in which there is simply a single additional fermion beyond the ↑and ↓-spin electrons, making an even more explicit connection to electronic CC theory.

III. APPLICATION
As a proof-of-principle, we consider a half-filled foursite Hubbard chain with an additional dipole coupled to a single photon cavity mode with frequency ω c (see Fig. 2 and Ref. [59]). Here we consider three values for the light-matter coupling parameter γ , representing weak (γ = 0.01), strong (γ = 0.07), and ultrastrong (γ = 0.2) coupling regimes [18,19]. The Hamiltonian of the model system is as in Eq. (1), but without the sum over α and witĥ wheren iσ =ĉ † iσĉ iσ denotes the density of a spin-σ electron on site i, and t 0 and U are the usual hopping and on-site repulsion constants. The dipole operator of the system is given byd = i d i (n i↑ +n i↓ ). The results for the ground-state energy E 0 and photonic mode occupation â †â in the system are summarized in Table I for different levels of CC theory and compared to full configuration interaction (FCI) results for the coupled electron photon system and to FCI and CCSD results in the limit of vanishing electron-photon coupling, FCI(0) and CCSD(0), respectively. We observe an increasing impact of the cavity on the ground-state properties of the system for growing coupling strengths. For instance, the gap between CC-SD-S-0 and CCSD(0) energies widens due to increasing importance of the dipole self-interaction term. This behavior is captured very well with polaritonic CC theory for all coupling strength as soon as coupled excitation are included, namely with CC-SD-S-D and CC-SD-S-DT approximations.
The photon-mode occupation â †â , which is now accessible with CC theory, is zero for CC-SD-S-0 for all coupling strength, which shows the intrinsic mean-field character of the CC-SD-S-0 approximation. The photon-mode occupation is also captured well as soon as coupled excitation are included. The best agreement with FCI is achieved for both observables with CC-SD-S-DT.
However, a considerable effect of the cavity on the molecular ground state is observed only in the ultrastrong coupling regime ( E 0 ≈ 0.02), which is captured well with CC theories that include coupled excitations. For the other two regimes this impact is, as was to be expected, rather small [30,45,60]. The power of CC theory therefore also lies in the treatment of excited (polaritonic) states that we show in the following.
A key experimentally accessible feature to study is the ground-state absorption spectrum. The matter absorption cross section is given by [17,61] where | k are many-body eigenstates ofĤ with energȳ hω k , ω is the frequency of incident light, and η is a (small) broadening parameter accounting for the finite lifetime of the state.
Here equation of motion CC (EOM-CC) theory will be used to access excited-state information [62][63][64][65]. In EOM-CC each excited state of the system is produced by applying a linear excitation operator to the correlated ground state: The coefficients entering the excitation operator R k are found by diagonalizingĤ in the subspace of excited states | μ  SinceĤ is non-Hermitian, we obtain a biorthogonal set of left (L j ) and right (R k ) operators satisfyinĝ and we define the left eigenstates as With these states, the absorption cross section can be approximated as A final detail is that the non-Hermitian form of the CC theory can lead to problems for close-lying or degenerate states (such as at conical intersections), that either have complex eigenvalues or exhibit significant overlaps among the right (or left) eigenvectors which can cause numerical instability. Both issues can be resolved following a correction method based on Ref. [66] and details can be found in Appendix C.
Leaving the weak-coupling case to Appendix A, we start our discussion of excited-state properties with the case of strong light-matter coupling of γ = 0.07. Absorption spectra for different approximations are shown in Fig. 3 as functions of cavity frequency ω c .
For reference, the panels at either end show the FCI spectrum for the purely electronic problem (left) and for the full light-matter problem (right). In between we show the results of EOM-CC calculations with various levels of truncation of the cluster operator. The bare electronic spectra [FCI(0) and CCSD(0)] feature horizontal nondispersive absorption lines; as expected, approximate EOM-CC accurately reproduces low-energy features in the spectrum.
When coupled to a cavity, the spectrum includes matter absorption lines and additional linear dispersive branches for one-photon (lines where ω = ω c ) and many-photon processes (ω = 2ω c , etc.).
In the strong-coupling case of Fig. 3, we see at most three-photon processes with significant amplitude. At points where there are matter excitations resonant with cavity modes we observe avoided crossings in the absorption lines. In these regions hybrid light-matter states (polaritons) are formed, accompanied by the signature Rabi splitting of an absorption peak. They can be clearly seen in various regions of the FCI spectrum, and are well captured by the EOM-CC approximation that includes coupled excitations. The systematic improvement in the CC treatment as the cluster operator is extended is clear, and remaining deviations from FCI are largely caused by the truncation of the electronic part of the cluster operator.
The exact (FCI) spectrum in the ultrastrong coupling case is much more complicated (Fig. 4). Additional features are seen beyond the simple combination of nondispersive matter lines and linear dispersive photon lines with avoided crossings.
These include induced transparencies (dark states) in regions of crossings of various absorption lines with the lowest one-photon line, and complicated structures in the high energy part of the spectrum where many multiphoton processes overlap in the spectrum.
The CC-SD-S-DT calculation captures much of this complex structure, with remaining deviations mainly caused by the truncation of the electronic cluster operator. The lowenergy part of the spectrum is reproduced extremely accurately, including the induced transparencies. The qualitative features of the high-energy part of the spectrum are also captured.

IV. CONCLUSION
Overall, we have shown that subtle light-matter correlations that appear in strong and ultrastrong cavity experiments can be captured in the framework of CC theory, paving the way for high-accuracy modeling and interpretation of experiments in these regimes.
Although this work has focused on model systems, extension to real ab initio Hamiltonians is straightforward. In this work we simulated the CC equations in the framework of exact diagonalization, but the formulation of the theory in terms of photon excitation operatorsτ n means that affordable polynomial scaling implementations will very closely mirror standard electronic CC codes, but with extra channels to describe the additional quantum degrees of freedom.
The individual photonic creation and annihilation operators are quadratic expressions in the excitation operatorsτ n [for exampleâ = n √ n (τ n−1τ † n )] so that the light-matter interaction in Eq. (1) becomes a four-point interaction, exactly analogous to the four-point two-particle interaction in the electronic ab initio Hamiltonian. Thus the equations for the polynomial-scaling implementation of polaritonic CC-SD-S-DT theory are effectively identical to a subset of the conventional CCSDT equations, but with removal of exchange diagrams and inclusion of alternative values in place of two-electron integrals. The theory scales as O(N 6 n max ) that is roughly the same as the O(N 6 ) scaling of conventional CCSD theory, since n max is usually much smaller than the number of electronic orbitals N.
Following this insight, we can conclude that the polaritonic CC method be applied to study linear and nonlinear optical processes occurring in molecules in the gas phase. We further believe that the applicability of our approach can be extended to liquid or condensed-matter systems by employing polarizable continuum solvation models [67,68] or QM/MM techniques [69].
Furthermore, the formalism developed here can be extended to CC theories for coupling of electrons to polarization modes, phonons, or thermal reservoirs, including coupling to multiple boson modes and boson-boson interactions. Work along these lines is underway in our groups.

APPENDIX B: CONVERGENCE STUDY, PHOTON NUMBER CUTOFF
The Fock space for the bosonic modes is infinite, so a reasonable truncation has to be performed in practice by imposing a maximum number of quanta per mode n max ph . Here we demonstrate convergence with respect to n max ph for the FCI spectra of Figs. 3-5. The absorption cross sections σ (ω) are plotted in Fig. 6 for the cavity frequency ω c = 1.028, where the cavity is resonant to the first absorption peak of the bare electronic system, and for different photonic cutoffs n max ph for weak, strong, and ultrastrong coupling and compared to the bare spectrum. We observe a convergence of the spectra at n max ph = 1 (weak coupling), n max ph = 4 (strong coupling), and n max ph = 7 (ultrastrong coupling), respectively. These values of n max ph were used for all results presented in this paper.

APPENDIX C: CORRECTION METHOD FOR CLOSE-LYING EIGENSTATES
Coupled-cluster EOM provides a powerful framework for treating excited states that is systematically improvable towards exact solutions to the Schrödinger equation. One downside is that at finite truncation of the cluster operator, the theory is non-Hermitian, and this can lead to difficulties for near-degenerate excitations. For example, the non-Hermitian model Hamiltonian (Ĥ) has different right and left eigenvectors, and pairs of right (or left) eigenvectors associated with a near degeneracy can become almost parallel. Second, the eigenvalues ofĤ can develop nonzero imaginary components.
The issues can be resolved based on the analysis of Köhn and Tajti [66], and we refer the interested reader to their work for a detailed discussion. Here we briefly outline their method to fully specify the calculations we present.
The eigenvalue E i of the similarity-transformed Hamilto-nianĤ is associated with left and right eigenvectors i | and | i that fulfill the biorthogonality condition The right eigenvectors are conventionally normalized, but generally nonorthogonal. The elements of the metric matrix S are S i j = i | j , so we have S ii = 1 and S i j = 0 with i = j. The left and right eigenvectors are related via This relation between the two causes problems, when right vectors become (almost) parallel S i j → 1. The norm of the corresponding left eigenvectors diverges and for S i j = 1 the inverse of S simply does not exist, as the right vectors | i do not span the full space anymore.
In order to compensate for this behavior, we employ here a correction scheme based on the method of Köhn et al. [66]. The idea is basically to rotate the close-lying eigenstates within their subspace such that their overlap becomes smaller and the implications described above become less pronounced.
For simplicity, consider only two close-lying eigenstates of H:Ĥ We build a subspace matrix and shift it by = ε 2 −ε 1 2 , 023262-6 such that A is traceless and has either purely real or purely imaginary eigenvalues. We further write the metric matrix as where S is the overlap of the two right vectors S = 1 | 2 . If S is complex, the state | 2 has to be rotated such that S becomes real: We now use the overlap matrix S to first orthogonalize the right-hand basis within the subspace spanned by | 1 and | 2 by multiplying it with S −1/2 , We then choose a new smaller overlap (and the corresponding angle ϑ = arccos ). We define it as function of the old overlap S and eigenvalue difference and further impose the condition that it is monotonous and is bounded from below by 0 and from above by a maximum value max . The specific form of the function is not that important and we take the one used in Ref. [66]: = max tanh (|S|/ max ), if real, tanh (1/(|S| · max )), if imaginary.
The orthogonalized eigenstates are then rotated again, this time to have the new overlap = cos ϑ. The corresponding rotation matrix reads To summarize the derivation above, the new subspace basis is obtained via Then the vectors | 1 , | 2 have to be normalized and the left eigenvectors have to be transformed accordingly. So far we have just rotated the eigenstates within their subspace and the eigenvalues remained unaffected by these operations. Now also the eigenvalues have to be adapted with respect to the new overlap. This is done with following formula: corr = ± | sin ϑ| | sin ϕ| 1, if real, −iS, if imaginary. (C11) We refer the reader to Ref. [66] for a detailed derivation. Throughout the paper we use max = 0.2 and apply the correction scheme described above as soon as < 0.05 or imaginary. These parameters can be adapted if needed.