We consider the non-relativistic photon-matter Hamiltonian, ĤPEN (r, q, R), hereafter referred to as the photonelectron-nuclear (PEN) Hamiltonian. As described in Ref. [ 34], this Hamiltonian is treated within the dipole approximation in the length gauge, i.e., it is derived by applying the Power-Zienau-Woolley 87 gauge transformation to the minimal coupling Hamiltonian in the Coulomb gauge , resulting in ĤPEN (r, q, R) = ĤM (r, R) + ĤP (q) + ĤPM (r, q, R) (1) with symbols r indicating the set of N e electronic spatial coordinates, R indicating the set of N N nuclear coordinates and q representing the displacement coordinates of the N P cavity modes. Atomic units ( = e 2 = m e = 1) will be used throughout. The matter (M) Hamiltonian is the sum of the nuclear kinetic energy operator TN (R) and of the Born-Oppenheimer (BO) Hamiltonian ĤBO (r, R), containing the electronic ki-netic energy and all the interactions Ve,n (r, R), namely