In the above, b † and b 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 γ/2 being related to the dissipation rate of the photonic degree of freedom . 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 μξ (x i ) + µ ξ nuc where μξ (x i ) is an operator that depends on electronic coordinates and within the Born-Oppenheimer approximation, we treat the Cartesian components of the nuclear dipole moment µ ξ nuc as functions 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 . In the NH-CQED-CIS theory we present, we will utilize an orbital basis that arises from solving the CQED-RHF equations arising from a Hermitian total Hamiltonian where only the real part of ω is retained. The CQED-RHF method has been described elsewhere , though we provide a brief outline in the Appendix. Also provided in the Appendix is an expansion of the Ĥdse term that explicitly shows where the 1-electron dipole, 1-electron quadrupole, and 2-electron dipole operator terms that appear in the NH-CQED-CIS matrix elements come from. We also note that the formulation presented here considers only a single photonic mode; generalizations to multiple photonic modes can be formulated by introducing additional frequencies and coupling parameters for each additional modes in the Ĥp , Ĥep , and Ĥdse terms. In particular, the multimode version would read