  Confined electromagnetic modes strongly couple to collective excitations in
ensembles of quantum emitters, producing light-matter hybrid states known as
polaritons. Under such conditions, the discrete multilevel spectrum of
molecular systems offers an appealing playground for exploring multiphoton
processes. This work contrasts predictions from the Tavis-Cummings (TC) model,
in which the material is a collection of two-level systems, with the
implications of considering additional energy levels with harmonic and
anharmonic structures. We discuss the exact eigenspectrum, up to the second
excitation manifold, of an arbitrary number $N$ of oscillators collectively
coupled to a single cavity mode in the rotating-wave approximation. Elaborating
on our group-theoretic approach [New J. Phys. 23, 063081 (2021)], we simplify
the brute-force diagonalization of a gigantic $N^2\times N^2$ Hamiltonian
(where $N=10^6-10^{10}$, as experiments suggest) to the diagonalization of, at
most, $4\times4$ matrices. We thoroughly discuss the eigenstates and the
consequences of weak and strong anharmonicities. Furthermore, we find resonant
conditions between bipolaritons and anharmonic transitions where two-photon
absorption can be enhanced. Finally, we conclude that energy shifts in the
polaritonic states induced by anharmonicities become negligible for large $N$.
Thus, calculations with a single or few emitters qualitatively fail to
represent the nonlinear optical response of the collective strong coupling
regime. Our work highlights the rich physics of multilevel anharmonic systems
coupled to cavities absent in standard models of quantum optics. We also
provide concise tabulated expressions for eigenfrequencies and transition
amplitudes, which should serve as a reference for future spectroscopic studies
of molecular polaritons.
