where the c coefficients are taken from the lowest root of the NH-CQED-CIS Hamiltonian matrix. We see an approximate quadratic dependence of the photon occupation in the ground-state as a function of relative coupling strengths independent of dissipa- . Photon occupations of the NH-CQED-CIS/cc-pVDZ groundstate of MgH+ at a bondlength of 2.2 Angstroms as a function of the relative coupling energy scale, g/Re( ω), and the dissipation energy scale γ. In this case, g ≈ ω 2 λ z µ z for MgH+ where λ z = 0.0125 atomic units, Re(h ω) = 4.75 eV. We see approximate quadratic scaling of the photon occupation with the relative coupling strength in the ground-state at all loss rates, in agreement with the analytical findings discussed in Ref. 6. tion energy scale (see Figure ). The results of this analysis agrees qualitatively with the analytical analysis of the groundstate populations of the ground-state wavefunctions of dissipative systems in the ultra-strong coupling regime by De Liberato in Ref. 6. An important point raised by this work is that features associated with the ultra-strong coupling regime (e.g. virtual photon occupation in the ground-state) persist even in the presence of very strong dissipation that would obscure the effects normally associated with strong-coupling (e.g. Rabi splitting) . Similar to our analysis illustrated in Figure , the exact analysis also showed photon occupation of the ground-state wavefunction increases quadratically with the relative coupling strength, where the relative coupling strength is defined as g/Re( ω), where g relates to the coupling energy scale. To determine if our NH-CQED-CIS obeys this scaling relationship for the photon occupation of the groundstate, we compute the photon occupations for the ground-state of MgH+ computed at the NH-CQED-CIS/cc-pVDZ level of theory with bondlength fixed at r = 2.2 Angstroms for relative coupling strengths g/Re( ω) in a range of values between 0 and 0.5, where in our case g ≈ ω 2 λ z µ z , where Re(h ω) = 4.75 eV throughout (resonant with the |X → |A transition at this geometry, λ z = 0.0125 atomic units, λ z ≈ 2.2 atomic units, and Im(h ω) is chosen relative to g. We emphasize this analysis involves coupling strengths that are much larger than those realized in any experiments we are aware of with molecular systems, and serve purely to analyze the behavior of our approximate NH-CQED-CIS method compared to the exact analysis provided in Ref. 6.