In Figure of the text, we see deviations between the polaritonic surfaces computed by the NH-CQED-CIS/cc-pVDZ method and those computed by a simpler 3x3 model Hamiltonian that is parameterized by CIS/cc-pVDZ calculations of MgH+. In Figure , we explore the eigenvalues of these two approaches for the lower-and upper-polariton surfaces of MgH+ with bondlength fixed at 2.2 Angstroms coupled to a photon with energy 4.75 eV with variable values of λ z between 0 and 0.025 atomic units, which is approximately twice the value considered in Figure . We see that deviations between the model Hamiltonian and the NH-CQED-CIS results become evident when the fundamental coupling strength λ has a magnitude between 0.01 and 0.015 atomic units, which is in the range of values we considered for Figure in the main text, and that is also in line with experimentally-realized coupling strengths reported in Ref.