Non-Hermitian Cavity Quantum Electrodynamics-Configuration Interaction Singles Approach for Polaritonic Structure with ab initio Molecular Hamiltonians

We combine ab initio molecular electronic Hamiltonians with a cavity quantum electrodynamics model for dissipative photonic modes and apply mean-field theories to the groundand excited-states of resulting polaritonic systems. In particular, we develop a restricted Hartree-Fock theory for the mean-field ground-state and a non-Hermitian configuration interaction singles theory for mean-field excited-states of the molecular system strongly interacting with a photonic mode, and apply these methods to several paradigmatic polaritonic systems. We leverage the Psi4Numpy framework to yield open-source and accessible reference implementations of these methods.


I. INTRODUCTION
The interaction between molecular excitations and nanoconfined photons can produce the requisite strong interactions for polaritonic chemistry 1-29 . Motivated by a desire to provide a realistic picture of the molecular structure under the influence of strong photonic interaction, there has been a recent surge in activity focused on merging ab initio molecular electronic structure theory with cavity quantum electrodynamics (ab initio CQED) to provide an accurate and predictive model of polaritonic chemistry 9,24,27,[30][31][32][33][34][35] . Such approaches can provide access to potential energy surfaces, couplings, and other properties of interest for simulating the structure and reactivity of polaritonoic chemical systems. Here we present a simple ab initio QED method for treating ground-and excited-polaritonic states with explicit inclusion of photonic lifetimes via a non-Hermitian cavity quantum electrodynamics -configuration interaction singles approach (NH-CQED-CIS). We implement this approach in the coherent state basis which results from solution of the CQED-Hartree-Fock (CQED-HF) equations. We endeavor to provide a detailed picture of the key equations and algorithmic considerations for both the CQED-HF and NH-CQED-CIS approach, and also provide reference implementations through the Psi4Numpy project 36 . We apply both methods to the analysis of polaritonic structure of several paradigmatic systems.

II. THEORY
We start with the Pauli-Fierz Hamiltonian in the dipole approximation in the length gauge, written in atomic units, following 24,27,31 :Ĥ wherê (2) withT e (x i ) denoting the electronic kinetic energy operator for electron i,V eN (x i ; X A ) the (attractive) coulomb operator for electron i and nucleus A,V ee (x i , x j ) the (repulsive) coulomb operator for electrons i and j, and V N,N is the total (repulsive) coulomb potential between all of the nuclei. Within the Born-Oppenheimer approximation, V N,N is a constant, the nuclear kinetic energy is neglected, and the electron-nuclear attraction depends parametrically on the fixed nuclear coordinates. The photonic contribution is captured by the complex energŷ and the photon-molecule intereaction contains a bilinear coupling term,Ĥ and a quadratic dipole self energy term, In the above,b † andb are the bosonic raising/lowering operators for the photonic degrees of freedom, andω = ω − i γ 2 is a complex frequency of the photon with the real part ω being related to the energy of the photon, and the imaginary part γ being related to the dissipation rate of the photonic degree of freedom 18,25,37 . The term µ represents the ground state molecular dipole expectation value which has cartesian components ξ ∈ {x, y, z}. A given ξ component of the dipole operator has the formμ ξ = ∑ N e iμ is an operator that depends on electronic coordinates and within the Born-Oppenheimer approximation, we treat µ ξ nuc (x A ) as a function of the nuclear coordinates rather than a quantum mechanical operator. Note that the shift of the Hamiltonian by µ results from the transformation to the coherent state basis 24 . To compute µ ξ , we can use a single Slater determinant approximation to the ground electronic state, where µ ξ ii = i|μ ξ |i denotes the ξ -component of the occupied molecular dipole integrals and ∑ is the ξcomponent of the nuclear dipole moment defined by the nuclear coordinates and atomic charges in the molecule. As we will see in the following section, solving the CQED-RHF equations will entail iterative updates to µ with the CQED-RHF orbitals.

A. CQED-RHF
As our first step in solving for the energy eigenstates of Eq. 1, we follow Ref. 24 and 27 and introduce a product wavefunction between an electronic Slater determinant (which in practice may be initialized using a canonical RHF wavefunction) and a zero-photon number state, To develop CQED-RHF theory, we examine the expectation value of Eq.1 with respect to Eq. 7, where we see that the terms involvingĤ p andĤ ep vanish, and the expectation value ofĤ e is analogous to the ordinary RHF energy. To evaluate the expectation value ofĤ dse , we can first expandĤ dse in terms of the dipole operator (with electronic and nuclear contributions) and dipole expectation values as follows: In the above expansion ofĤ dse we have specifically indicated that the product of electronic dipole operators contains 2electron contributions when i = j, and 1-electron quadrupole contributions when i = j. The quadrupole contributions arise from the fact thatμ ξ ( Furthermore, a one-electron term arises that contains the electronic dipole operator scaled by λ · µ nuc − λ · µ , where again µ will be iteratively updated during the QED-RHF procedure.
To solve the CQED-RHF equations, the additional oneelectron terms above will be added to H core and the additional two-electron terms above will be included in the densitymatrix dependent terms in the Fock operator: where and leading to the total QED-RHF energy being where For clarity, we briefly outline the CQED-RHF algorithm below: B. Non-Hermitian QED-CIS in the CQED-RHF basis A mean-field description of the excited states of the molecular system strongly interacting with photonic degrees of freedom may be obtained through a configuration interaction singles (CIS) ansatz for the excited-states. Here we formulated a non-Hermitian version of such an ansatz, NH-CQED-CIS, that incorporates the dissipative features of the photonic degrees of freedom. In our presentation, we formulate NH-CQED-CIS in the coherent state basis using the orbitals that result from the CQED-RHF approach outlined above.
The polaritonic energy eigenfunctions for state I in the NH-CQED-CIS ansatz can be written as a linear combination of the CQED-RHF ground-state and products of all possible single excitations out of the CQED-RHF ground state. The CQED-RHF ground state involves the product of an electronic Slater determinant with the photon vacuum state, single excitations can occur as electronic excitations from an occupied orbital i to a virtual orbital a, the raising of the photon number state from |0 → |1 , or both. We therefore write the NH-QED-CIS wavefunction for state I as where the coefficients c denote the contribution of a given term to the wavefunction, where we have denoted the electronic excitations in the subscript and the photonic excitations in the superscript of these coefficients. These coefficients, and the corresponding energy eigenvalues for a given NH-CQED-CIS state I may be obtained by diagonalizing the Hamiltonian matrix built in the basis of the CQED-RHF and singly-excited states from Eq. 15. We use a spin-adapted basis of singly excited electronic states, such that |Φ a i = 1 √ 2 |Φ aα iα + |Φ aβ iβ . There are three classes of matrix elements that contribute to the Hamiltonian matrix; we write these matrix elements after shifting the total Hamiltonian matrix in Eq. 1 by E CQED−RHF . Matrix elements involving the CQED-RHF electronic Slater determinant |Φ 0 and photonic states |s and |t , where s,t ∈ {0, 1} involve only the (complex) photonic energy, Matrix elements coupling |Φ 0 |s to |Φ a i |t involve only thê H ep contributions: Matrix elements coupling different singly excited electronic and/or photonic states involve all terms of the Hamiltonian, including the usual CIS terms:

III. REFERENCE IMPLEMENTATIONS
We provide reference implementations using Psi4Numpy 36 , which provides a simple NumPy interface to the Psi4 38 quantum chemistry engine. The code for these reference implementations can be freely accessed in the hilbert package 39,40 . Furthermore, to provide a no-installation option for interested users to experiment with these implementations, we utilize the ChemCompute project 41 to host the illustrative calculations discussed in the Results section below. Interested users can navigate to https://chemcompute.org/register to register for a free ChemCompute account. Following registration, interested users can run calculations described in Figure 1 and Table I in the results section using the link in Ref. 42, and can run calculations described in Figure 2 and 3 using the link within Ref. 43.

IV. RESULTS
We apply the CQED-RHF and NH-CQED-CIS approaches to a few simple paradigmatic polaritonic chemical systems.
First we examine the ground-state of formaldehyde strongly coupled to a single photon mode, which has been explored in by several groups that have been developing density functional theory-based ab initio-QED methods 9,35 . We optimize the geometry of lone formaldehyde at the RHF/cc-pVDZ level and perform all calculations at that geometry. At this level, the RHF ground-state has a dipole moment oriented purely along the z−axis with µ z = −1.009 atomic units. We examine the impacts of photon polarization along both the y−, z− (or both) axes in Figure 1. In particular, we scan values of the magnitude of |λ | in the range [0, 0.2] atomic units in increments of 0.02 atomic units for three different polarizations: λ y = (0, |λ |, 0), λ z = (0, 0, |λ |), and λ yz = 0, |λ | √ 2 , |λ | √ 2 . We note that the electric field vector enters explicitly into the CQED-RHF equations while the photon frequency does not, so we do not report a specific value ofω for these calculations.
We see, not surprisingly, that the CQED-RHF energy increases monotonically with the magnitude λ z , and interestingly, shows the same qualitative behavior (albeit slightly less dramatically) with increasing magnitude of λ y (see top panel of Figure 1). We first decompose the total CQED-RHF energy into the canonical RHF contribution, and the Pauli-Fierz contribution that includes the remaining terms in the CQED-RHF Fock operator traced against the converged CQED-RHF density matrix (with elements D µν denoted above).
Not surprisingly, we see that the Pauli-Fierz contribution shows the most dramatic scaling with the electric field magnitude, comprising the majority of the energy change, but we also see that the canonical RHF contribution is modifed by the electric field, as well. We present a more granular decomposition of these energetic contributions only for the largest field magnitude |λ | = 0.2 atomic units in Table I. From Eq. 7, we see that only the photon vacuum contributes to the CQED-RHF wavefunction. However, in Eq. (17), we see that the CQED-RHF wavefunction can couple to states which involve singly-excited electronic configurations and singly-occupied photon states. This coupling can lower the energy of the lowest energy eigenstate of the CQED-CIS Hamiltonian relative to the ground-state determined by the CQED-RHF method.
We illustrate this point through the ground-state potential  energy surface of the MgH + diatomic cation coupled to a photon that is on resonance with the ground-state singlet (|X ) to first excited-state singlet (|A ) when the MgH + bondlength is approximately 2.2 Angstroms 23 . We consider the molecule to be oriented along with z axis, which is parallel to the polarization vector of the photon that has magnitude λ z = 0.075  2. Ground-state potential energy surface of MgH+ coupled to a photon resonant with X → A transition. We see that the photonic coupling raises the energy relative to the lone molecular system (as computed at the RHF/cc-pVDZ level). Coupling between the CQED-RHF wavefunction and singly-excited electronic and photonic configurations stabilizes the CQED-CIS ground-state relative to the CQED-RHF ground-state. atomic units. We compute the potential energy surface of this system at the RHF/cc-pVDZ, CQED-RHF/cc-pVDZ, and CQED-CIS/cc-pVDZ level and plot the results in Figure 2.
As a final illustrative example, we consider the upperpolariton (|UP ) and lower-polariton (|LP ) states that emerge from coupling MgH + to a photon resonant with the |X → |A transition. We again consider a z-polarized photon with magnitude λ z = 0.0125 atomic units. This time, we allow the photon to have a complex energy where the imaginary part accounting for the finite lifetime can also be related to the energy uncertainty of the photon. We consider the real part of the photon energy to be 4.75 eV, and we consider the imaginary part to be either 0 eV or 0.22 eV as shown in the top and bottom panels of Figure 3, respectively. In addition to computing these polariton surfaces at the CQED-CIS/cc-pVDZ level, we also fit a 3-level model Pauli-Fierz Hamiltonian from ordinary CIS/cc-pVDZ potential energy surfaces: The polaritonic surfaces obtained from diagonalizing Eq. (20) are referred to as the 'Model LP' and 'Model UP' surfaces in Figure 3. We see with a pure real photon energy, the |LP and |UP surfaces experience a strong splitting in the region where the |X, 1 state (the ground-state plus a photon) crosses the |A, 0 state (the first excited-state without a photon). However,  for the strongly dissipative photon, we see that both the model and CQED-CIS curves closely approximate the CIS curves for the lone molecules, which signals that this system is not in the strong-coupling regime because of the lossiness associated with the photon 25 .

V. CONCLUSIONS
We combined ab initio molecular electronic Hamiltonians with a cavity quantum electrodynamics (CQED) model for dissipative photonic modes and apply mean-field theories to the ground-and excited-states of resulting polaritonic systems. In particular, we developed a restricted Hartree-Fock theory for the mean-field ground-state and a non-Hermitian configuration interaction singles theory for meanfield excited-states of the molecular system strongly interacting with a photonic mode, and applied these methods to several paradigmatic polaritonic systems. We leveraged the Psi4Numpy 36 framework to yield open-source and accessible reference implementations of these methods.

ACKNOWLEDGMENTS
JM and JJF Acknowledges support from the Research Corporporation for Scientific Advancement Cottrell Scholar Award and the NSF CAREER Award CHE-2043215. Computational resources were provided in part by the MERCURY consortium (http://mercuryconsortium.org/) under NSF grants CHE-1229354, CHE-1662030, and CHE-2018427. We gratefully acknowledge A. E. DePrince III for numerous helpful discussions and access to benchmark results for the CQED-RHF method.