Note that under the vibrational eigenbasis {|ν i } (see Eq. 22), the Ĥren operator will be a constant matrix because it is only a function of R. As a consequence of the above system-bath partition, the system keeps its dimension the same as the bare matter part, which has greatly reduced computational cost compared to the conventional treatment. Moreover, no truncation approximation for photon Fock states is explicitly made, and the full Hilbert space of the cavity subsystem is taken into account.