# Generalization of the Tavis-Cummings model for multi-level anharmonic
systems: insights on the second excitation manifold

J.Campos-Gonzalez-Angulo  Department of Chemistry and Biochemistry, University
of California San Diego, La Jolla, California 92093, USA  J. Yuen-Zhou
<EMAIL_ADDRESS>http://yuenzhougroup.ucsd.edu

###### Abstract

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\times 4$ 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.

††preprint: AIP/123-QED

##  I Introduction

The strong coupling (SC) of confined photonic modes and excitations in
semiconductors and molecular materials leads to a wealth of interesting
chemical and condensed-matter physics phenomena [1, 2, 3, 4, 5] such as room-
temperature Bose-Einstein condensation,[6, 7, 8, 9] long-range energy
transfer,[10, 11] modification of chemical reactivity,[12, 13, 14] and quantum
information processing.[15, 16]. Strong coupling leads to hybrid light-matter
quantum modes known as polaritons, originally characterized for electronic
excitations in bulk materials,[17, 18] and later on studied in cavity
systems.[19, 20] Within linear response, the formalism to understand these
systems is naturally based on a harmonic approximation for the material
degrees of freedom.[18] However, to understand nonlinear optical properties,
the material degrees of freedom are typically modeled as two-level quantum
emitters, leading to the very well-studied Jaynes-Cummings and Rabi models in
the case of a single emitter,[21, 22] or the Tavis-Cummings (TC) and Dicke
models in the collective regime.[23, 24] Recently, there has been an
increasing interest in the properties of vibrational polariton systems
resulting from the strong coupling of infrared (IR) cavity modes and ensembles
of localized high-frequency vibrational modes in molecules in condensed
phases.[25, 26, 27, 28, 29, 30, 31, 32] The multilevel anharmonic spectrum of
these vibrational modes implies that their accurate description should invoke
more than two levels per emitter. Indeed, consideration of more realistic
vibrational SC systems involving anharmonicities has been the subject of
theoretical explorations, such as the use of single-molecule models to explain
cavity-induced modifications to chemical reactivity,[33, 34, 35, 36, 37] as
well as the use of many-molecule models to explain non-linear response
experiments.[38, 39, 40, 41, 42, 43, 28] In particular, recent work [44, 45,
46, 47, 48, 49] has elucidated novel multiphoton absorption phenomena where
the TC description is insufficient, and a multilevel anharmonic spectrum of
the material is essential. These effects go beyond vibrational SC and should
have analogues in other electromagnetic ranges, such as those systems under
electronic SC where more than two electronic states per quantum emitter must
be considered.

While theoretical complexities emerging from the multilevel anharmonic
spectrum of the oscillators are expected, they can be overcome by taking
advantage of the permutational symmetries derived from the assumption that the
emitters behave identically.[50, 51, 52, 53] For regimes in which the
counterrotating (non-energy conserving) terms of the coupling can be
disregarded, the Hamiltonian can be separated according to the total number of
excitations allocated in the many-body states that conform the basis of the
Hilbert space.[54, 55] The general strategy to simplify and solve the
Schrödinger equation for this instance can be found in ref. 56, work on which
the present manuscript elaborates by describing in detail the wavefunctions
and the energy spectrum resulting from diagonalizing the Hamiltonian of a
collection of anharmonic multilevel quantum emitters dipolarly coupled to a
single cavity mode. In particular, we focus on the structure of the subspace
built from states bearing two excitations and compare it with the well-known
solutions of the first excitation manifold. For systems under collective SC,
the number of emitters is very large at $N\approx 10^{6}-10^{10}$, and a
brute-force numerical diagonalization to solve for the eigenspectrum of these
hybrid light-matter systems is unattainable. Our method exploiting the
permutational symmetry arising from the consideration of $N$ identical
emitters dramatically reduces the problem to the diagonalization of very small
matrices.

The organization of this manuscript is as follows: in section II we revisit
the first-excitation manifold and present the tools and vocabulary of
polaritonic states that permeate the remaining of the paper. Section III opens
by introducing generalities of the states with two excitations and proceeds
with the presentation of the eigenstates for the TC model. These solutions are
contrasted with those obtained after including an additional energy level per
emitter, in the harmonic and anharmonic regimes, respectively. Tabulated
expressions are also provided so that this work serves as a future reference
for work on molecular polariton spectroscopy. Finally, the conclusions are
presented in section IV.

##  II The first excitation manifold

For a collection of $N$ identical multi-level quantum emitters interacting
with a confined electromagnetic mode of frequency $\omega_{0}$ in the regime
where the rotating-wave approximation is valid, the total number of
excitations in the system, $n_{\textrm{exc}}$, is a conserved quantity that
defines the so-called excitation manifolds.

When $n_{\textrm{exc}}=0$, all the components of the ensemble are in their
respective ground states, and the system is characterized by the ($N+1$)-body
state $\ket{0}$.

A basis can be defined for states with $n_{\textrm{exc}}=1$, such that the
excitation is localized on each of the emitters. To be specific, we have

$\displaystyle\ket{1_{0}}=$ $\displaystyle\hat{a}_{0}^{\dagger}\ket{0},$ (1a)
and $\displaystyle\ket{1_{i>0}}=$
$\displaystyle\quantity(\hat{\sigma}_{i}^{(0)})^{\dagger}\ket{0},$ (1b)

where $\hat{a}_{0}^{\dagger}$ is the creation operator acting on the EM mode,
and $\hat{\sigma}_{i}^{(v)}$ is the local operator acting on the $i$th emitter
projecting the state $\ket{v+1}$ onto $\ket{v}$. The notation in eq. 1 implies
that all particles not explicitly indicated are in their respective ground
states.

The Hamiltonian $\hat{H}_{1}$ describes the system in the singly excited
manifold. In the basis defined in eq. 1, the matrix elements of $\hat{H}_{1}$,
in units of $\hbar$, are given by

$\matrixelement{1_{i}}{\hat{H}_{1}}{1_{j}}=\begin{cases}\omega_{0}&\text{if
}i=j=0\\\ \omega_{10}&\text{if }i=j>0\\\ g_{01}&\text{for }i=0,\enskip j>0\\\
g_{10}&\text{for }i>0,\enskip j=0\\\ 0&\text{otherwise}\end{cases},$ (2)

where $\omega_{uv}$ is the excitation frequency between molecular energy
levels $u$ and $v$, and
$g_{uv}=\sqrt{\frac{\omega_{uv}}{2\hbar\epsilon_{0}\mathcal{V}}}\bra{u}\hat{\mu}\ket{v}$
is a coupling constant with $\mathcal{V}$ the mode volume, and $\hat{\mu}$ the
transition dipole moment operator ($\hbar$ is the reduced Planck’s constant,
and $\epsilon_{0}$ is the permittivity of the vacuum).

The system remains unchanged upon permutations of the emitters; therefore, its
description can be simplified using the SU(2) collective operators[57, 58]

$\displaystyle\hat{J}_{-}^{(v)}=$
$\displaystyle\sum_{i=1}^{N}\hat{\sigma}_{i}^{(v)},$ (3a)
$\displaystyle\hat{J}_{+}^{(v)}=$
$\displaystyle\quantity(\hat{J}_{-}^{(v)})^{\dagger},$ (3b) and
$\displaystyle\hat{J}_{0}^{(v)}=$
$\displaystyle\frac{1}{2}\quantity[\hat{J}_{+}^{(v)},\hat{J}_{-}^{(v)}],$ (3c)

which adhere to the angular momentum algebra, i.e.,

$\displaystyle\quantity[\hat{J}_{0}^{(v)},\hat{J}_{\pm}^{(v)}]=$
$\displaystyle\hat{J}_{\pm}^{(v)},$ (4a)
$\displaystyle\quantity[\hat{J}_{+}^{(v)},\hat{J}_{-}^{(v)}]=$ $\displaystyle
2\hat{J}_{0}^{(v)}.$ (4b)

Since these operators are permutationally invariant, they carry the totally-
symmetric or trivial irreducible representation (irrep) of the symmetric group
$S_{N}$,[50, 59] which will be denoted by A throughout this manuscript.

The global ground state, $\ket{0}$, is also permutationally invariant;
therefore, the states $\ket{1_{0}}$ and

$\begin{split}\ket{1_{\textrm{A}}}=&\frac{\hat{J}_{+}^{(0)}}{\sqrt{N}}\ket{0}\\\
=&\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{1_{i}},\end{split}$ (5)

carry the totally-symmetric irrep as well. On the other hand, the states

$\ket{1_{\textrm{B}(k)}}=\sum_{n=1}^{N}c_{n}^{(k)}\ket{1_{n}},$ (6)

where the coefficients $c_{n}^{(k)}$ fulfill $\sum_{n=1}^{N}c_{n}^{(k)}=0$ and
$\sum_{n=1}^{N}c_{n}^{(k)*}c_{n}^{(k^{\prime})}=\delta_{kk^{\prime}}$, span
the $N-1$-dimensional Hilbert subspace corresponding to the standard irrep of
$S_{N}$, denoted by B, which is orthogonal to the totally-symmetric one.

Since the subspace that the wavefunctions $\ket{1_{\textrm{B}(k)}}$ span is
highly degenerate, the choice of basis is not unique.[60] The most prominent
examples of these bases are the Fourier basis:

$c_{n}^{(k)}=\frac{1}{\sqrt{N}}\exp\left(2\pi i\frac{(k-1)n}{N}\right),$ (7)

in which $2\pi(k-1)/N$ is a wave vector, and the Schur-Weyl basis, where

$c_{n}^{(k)}=\frac{\alpha_{n}^{(k)}}{\sqrt{k\quantity(k-1)}},$ (8)

with

$\alpha_{n}^{(k)}=\begin{cases}-1&\text{for }1\leq n<k\\\ k-1&\text{if }n=k\\\
0&\text{for }k<n\leq N\end{cases}.$ (9)

The states in the symmetry-adapted basis in eqs. 5 and 6 are also known as
Dicke states,[23, 61] $\ket{J_{(0)},M_{J}^{(0)}}$, with quantum numbers given
by the equations

$\displaystyle\hat{J}_{(0)}^{2}\ket{J_{(0)},M_{J}^{(0)}}$
$\displaystyle=J_{(0)}\quantity(J_{(0)}+1)\ket{J_{(0)},M_{J}^{(0)}},$ (10a)
and $\displaystyle\hat{J}_{0}^{(0)}\ket{J_{(0)},M_{J}^{(0)}}$
$\displaystyle=M_{J}^{(0)}\ket{J_{(0)},M_{J}^{(0)}},$ (10b)

where

$\hat{J}_{(0)}^{2}=\quantity(\hat{J}_{0}^{(0)})^{2}+\frac{1}{2}\quantity(\hat{J}_{+}^{(0)}\hat{J}_{-}^{(0)}+\hat{J}_{-}^{(0)}\hat{J}_{+}^{(0)}).$
(11)

For reasons that will become evident later, we found our labeling of the
irreps more convenient; however, the identification of our notation with the
corresponding quantum numbers from angular momentum operators can be found in
table 1.

Table 1: SU(2) and SU(3) identifiers of the molecular symmetry-adapted states. State | $J_{0}$ | $M_{J}^{(0)}$ | $y$ | SU(3) multiplet
---|---|---|---|---
$\ket{0}$ | $N/2$ | $-N/2$ | $N/3$ | $\quantity(N,0)$
$\ket{1_{\textrm{A}}}$ | $N/2$ | $1-N/2$ | $N/3$ | $\quantity(N,0)$
$\ket{1_{\textrm{B}(k)}}$ | $N/2-1$ | $1-N/2$ | $N/3$ | $\quantity(N-2,1)$
$\ket{1_{\textrm{A}}^{2}}$ | $N/2$ | $2-N/2$ | $N/3$ | $\quantity(N,0)$
$\ket{1_{\textrm{B}(k)}^{2}}$ | $N/2-1$ | $2-N/2$ | $N/3$ | $\quantity(N-2,1)$
$\ket{1_{\textrm{C}(k\ell)}^{2}}$ | $N/2-2$ | $2-N/2$ | $N/3$ | $\quantity(N-4,2)$
$\ket{2_{\textrm{A}}}$ | $\quantity(N-1)/2$ | $\quantity(1-N)/2$ | $N/3-1$ | $\quantity(N,0)$
$\ket{2_{\textrm{B}(k)}}$ | $\quantity(N-1)/2$ | $\quantity(1-N)/2$ | $N/3-1$ | $\quantity(N-2,1)$

In the basis
$\quantity{\ket{1_{0}},\ket{1_{\textrm{A}}},\ket{1_{\textrm{B}(2)}},\ket{1_{\textrm{B}(3)}},\ldots,\ket{1_{\textrm{B}(N)}}}$,
the Hamiltonian becomes

$\begin{split}H_{1}=&\begin{pmatrix}\omega_{0}&\sqrt{N}g_{01}&0&0&\cdots&0\\\
\sqrt{N}g_{10}&\omega_{10}&0&0&\cdots&0\\\ 0&0&\omega_{10}&0&\cdots&0\\\
0&0&0&\omega_{10}&\cdots&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&0&\cdots&\omega_{10}\end{pmatrix}\\\
=&\mathbf{H}_{1}^{(\textrm{A})}\oplus\quantity(\mathbf{1}_{N-1}\otimes\mathbf{H}_{1}^{(\textrm{B})}),\end{split}$
(12)

where $\mathbf{H}_{1}^{(\textrm{A})}$ is a $2\times 2$ matrix,
$\mathbf{H}_{1}^{(\textrm{B})}=\omega_{10}$, and $\mathbf{1}_{d}$ is the
$d$-dimensional identity matrix.

From eq. 12, it becomes clear that the dark states $\ket{1_{\textrm{B}(k)}}$
are eigenfunctions of $\hat{H}_{1}$. The other eigenstates are in the totally-
symmetric sub-space:

$\ket{1_{\pm}}=\pm h_{\pm}\ket{1_{0}}+h_{\mp}\ket{1_{\textrm{A}}},$ (13)

where $+$ and $-$ label the upper and lower polariton, respectively, and
$h_{\pm}=\sqrt{\frac{1}{2}\quantity(1\pm\frac{\Delta}{\Omega_{10}})}$ are
known as Hopfield coefficients. The corresponding eigenfrequencies are

$\begin{split}\expectationvalue{\hat{H}_{1}}{1_{\pm}}=&\omega_{\pm}\\\
=&\frac{\omega_{0}+\omega_{10}\pm\Omega_{10}}{2},\end{split}$ (14)

where $\Omega_{10}=\sqrt{\Delta^{2}+4Ng_{01}^{2}}$ is the Rabi frequency, and
$\Delta=\omega_{0}-\omega_{10}$ is the detuning. Notice that the use of
permutational symmetry arguments to derive the analytical solution of the
eigenspectrum of $\hat{H}_{1}$ is a drastic simplification of a brute-force
numerical diagonalization of a $(N+1)$-dimensional matrix, which in the limit
of $N=10^{6}-10^{10}$, as is the case of collective SC, becomes intractable.
We shall see how an analogous simplification can be carried out for the
second-excitation manifold.

The intensities of spectroscopic signals depend on the off-diagonal matrix
elements of the collective dipole operator
$\hat{\bm{\mu}}=\sum_{i=1}^{N}{\hat{\mu}_{i}}$, and the photon mode creation
operator $\hat{a}_{0}^{\dagger}$. For transitions between the ground state,
$\ket{0}$, and the eigenstates in the first excitation manifold we have

$\displaystyle\bra{1_{\pm}}{\hat{\bm{\mu}}}$
$\displaystyle\ket{0}=\sqrt{N}h_{\mp}\mu_{10},$ (15)
$\displaystyle\bra{1_{\textrm{B}(k)}}{\hat{\bm{\mu}}}$
$\displaystyle\ket{0}=0,$ (16)
$\displaystyle\bra{1_{\pm}}{\hat{a}_{0}^{\dagger}}$ $\displaystyle\ket{0}=\pm
h_{\pm},$ (17) and
$\displaystyle\bra{1_{\textrm{B}(k)}}{\hat{a}_{0}^{\dagger}}$
$\displaystyle\ket{0}=0,$ (18)

where $\mu_{10}=\matrixelement{1}{\hat{\mu}}{0}$ is the dipole moment for the
$0\to 1$ transition of the bare emitters. Equation 16 can be understood under
the consideration that the operator $\hat{\bm{\mu}}$ is totally-symmetric, and
the states $\ket{1_{\textrm{B}(k)}}$ carry an orthogonal irrep.

##  III The second excitation manifold

The states with two quanta can be of the form

$\ket{1_{i}1_{j}}=\quantity(\hat{\sigma}_{j\neq
i}^{(0)})^{\dagger}\ket{1_{i}},$ (19)

in which the $i$th and the $j$th particle are both in their first excited
state, or

$\displaystyle\ket{2_{0}}=$
$\displaystyle\frac{\hat{a}_{0}^{\dagger}}{\sqrt{2}}\ket{1_{0}},$ (20a) and
$\displaystyle\ket{2_{i>0}}=$
$\displaystyle\quantity(\hat{\sigma}_{i}^{(1)})^{\dagger}\ket{1_{i}},$ (20b)

where the $i$th particle is in its second excited state. As in the singly
excited manifold, all particles not explicitly indicated are in their
respective ground state. Notice that a brute-force numerical diagonalization
of the second-excitation manifold Hamiltonian $\hat{H}_{2}$ would require
diagonalization of ${N+2\choose 2}$-dimensional matrices, that is unattainable
for the large number $N$ of emitters which concern us in the context of
collective SC. Instead, as in the previous section, we will exploit the
permutational symmetry arising from the consideration of identical emitters
coupled to the cavity.

Just as before, a symmetry-adapted basis can be defined such that there are
wavefunctions carrying the totally-symmetric irrep: $\ket{2_{0}}$,

$\begin{split}\ket{1_{0}1_{\textrm{A}}}=&\hat{a}_{0}^{\dagger}\ket{1_{\textrm{A}}}\\\
=&\frac{\hat{J}_{+}^{(0)}}{\sqrt{N}}\ket{1_{0}},\end{split}$ (21)
$\begin{split}\ket{1_{\textrm{A}}^{2}}=&\frac{\hat{J}_{+}^{(0)}\ket{1_{\textrm{A}}}}{\sqrt{2\quantity(N-1)}}\\\
=&\sqrt{\frac{2}{N(N-1)}}\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\ket{1_{i}1_{j}},\end{split}$
(22)

and

$\begin{split}\ket{2_{\textrm{A}}}=&\hat{J}_{+}^{(1)}\ket{1_{\textrm{A}}}\\\
=&\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{2_{i}}.\end{split}$ (23)

There are also wavefunctions carrying the standard irrep:

$\ket{1_{0}1_{\textrm{B}(k)}}=\hat{a}_{0}^{\dagger}\ket{1_{\textrm{B}(k)}},$
(24)
$\begin{split}\ket{1_{\textrm{B}(k)}^{2}}=&\frac{\hat{J}_{+}^{(0)}}{\sqrt{N-2}}\ket{1_{\textrm{B}(k)}}\\\
=&\sum_{m=1}^{N-1}\sum_{n=m+1}^{N}c_{mn}^{(k)}\ket{1_{m}1_{n}},\end{split}$
(25)

and

$\begin{split}\ket{2_{\textrm{B}(k)}}=&\hat{J}_{+}^{(1)}\ket{1_{\textrm{B}(k)}}\\\
=&\sum_{n=1}^{N}c_{n}^{(k)}\ket{2_{n}},\end{split}$ (26)

where the coefficients $c_{n}^{(k)}$ are the same as in eqs. 7 and 8, and the
coefficients $c_{mn}^{(k)}=\quantity(c_{m}^{(k)}+c_{n}^{(k)})/\sqrt{N-2}$
fulfill

$\displaystyle\sum_{m=1}^{N-1}\sum_{n=m+1}^{N}c_{mn}^{(k)}=$ $\displaystyle
0,$ (27a) and
$\displaystyle\sum_{m=1}^{N-1}\sum_{n=m+1}^{N}c_{mn}^{(k)}c_{mn}^{(k^{\prime})*}=$
$\displaystyle\delta_{kk^{\prime}}.$ (27b)

In the Fourier basis, these coefficients are

$c_{mn}^{(k)}=\frac{2\exp\left[\pi
i(k-1)\frac{m+n}{N}\right]}{\sqrt{N(N-2)}}\cos\left(\pi(k-1)\frac{m-n}{N}\right),$
(28)

while in the Schur-Weyl basis,

$c_{mn}^{(k)}=\frac{\alpha_{mn}^{(k)}}{\sqrt{(N-2)k(k-1)}},$ (29)

where

$\alpha_{mn}^{(k)}=\begin{cases}-2&\text{for }1\leq m<k,\enskip 1<n<k\\\
k-2&\text{for }1\leq m<k,\enskip n=k\\\ -1&\text{for }1\leq m<k,\enskip
k<n\leq N\\\ k-1&\text{if }m=k,\enskip k<n\leq N\\\
0&\text{otherwise}\end{cases}.$ (30)

Finally, there are wavefunctions,

$\ket{1_{\textrm{C}(k\ell)}^{2}}=\sum_{m=1}^{N-1}\sum_{n=m+1}^{N}c_{mn}^{(k\ell)}\ket{1_{m}1_{n}},$
(31)

carrying the $N(N-3)/2$-dimensional irrep that is orthogonal to both the
totally-symmetric and the standard irreps. Their coefficients fulfill

$\displaystyle\sum_{m=1}^{N-1}\sum_{n=m+1}^{N}c_{mn}^{(k\ell)}=$
$\displaystyle 0,$ (32a)
$\displaystyle\sum_{m=1}^{N-1}\sum_{n=m+1}^{N}c_{mn}^{(k\ell)}c_{mn}^{(k^{\prime})*}=$
$\displaystyle 0,$ (32b) and
$\displaystyle\sum_{m=1}^{N-1}\sum_{n=m+1}^{N}c_{mn}^{(k\ell)}c_{mn}^{(k^{\prime}\ell^{\prime})*}=$
$\displaystyle\delta_{kk^{\prime}}\delta_{\ell\ell^{\prime}}.$ (32c)

These coefficients, in the Fourier basis, have the general form

$c_{mn}^{(k\ell)}=\sum_{q}\sum_{K(q)}\psi(K_{k\ell},q_{k\ell})\textrm{e}^{\pi
iK_{k\ell}\frac{m+n}{N}}\cos\left(\pi q_{kl}\frac{m-n}{N}\right),$ (33)

where $K_{k\ell}$ labels a center-of-mass wave number, $q_{k\ell}$ identifies
a relative wave number, and $\psi(q_{k\ell})$ is a coefficient resulting from
symmetrization. The structure of this basis is discussed in greater detail in
appendix B. In the Schur-Weyl basis, the coefficients are given by

$c_{mn}^{(k\ell)}=-\frac{\alpha_{mn}^{(k\ell)}}{\sqrt{k(k-1)(\ell-2)(\ell-3)}},$
(34)

with

$\alpha_{mn}^{(k\ell)}=\begin{cases}-2&\text{for }1\leq m<k,1<n<k\\\
k-2&\text{for }1\leq m<k,n=k\\\ -1&\text{for }1\leq m<k,k<n<\ell\\\
\ell-3&\text{for }1\leq m<k,n=\ell\\\ k-1&\text{for }m=k,k<n<\ell\\\
-(k-1)(\ell-3)&\text{if }m=k,n=\ell\\\ 0&\text{otherwise}\end{cases}.$ (35)

We note that the involvement of $\hat{J}_{+}^{(1)}$ in eqs. 23 and 26 implies
that the SU(2) algebra with operators $\hat{J}_{x}^{(0)}$ is insufficient to
characterize the states with two excitations in the same emitter.[62] These
states need labels from the SU(3) algebra for their correct identification. In
table 1, we present these labels, including the eigenvalues, $y$, of the
operator

$\hat{Y}=\frac{2}{3}\quantity(\quantity[\hat{J}_{+}^{(1)},\hat{J}_{-}^{(1)}]-\hat{J}_{0}^{(0)}),$
(36)

usually known as hypercharge in the literature of particle physics.[63] We
also require labels for the SU(3) irrep or multiplet; the methods to identify
them can be found elsewhere.[64] As we can see, $J_{0}$, is useful to identify
the symmetry of states with at most one excitation per emitter, but is not
capable of classifying the remaining states. On the other hand, the irreps of
SU(3) thoroughly fulfill this job. However, keeping track of that many quantum
numbers is cumbersome for the purposes of this work, and the A, B, C scheme
here introduced conveys the important information in a more condensed way.

Up to this point, we know that the matrix element

$\matrixelement{\Psi_{\Gamma}}{\hat{H}_{2}}{\Psi_{\Gamma^{\prime}}}=0,$ (37)

if $\Gamma\neq\Gamma^{\prime}$, where $\hat{H}_{2}$ is the Hamiltonian with
$n_{\textrm{exc}}=2$, and $\ket{\Psi_{\Gamma}}$ is a symmetrized state
carrying the irrep $\Gamma$ which can be either A, B or C.

The form of the Hamiltonian in the subspaces carrying the irreps A and B
depends on the considered spectrum of the emitters, and will be discussed in
the following sections. On the other hand, the states
$\ket{1_{\textrm{C}(k\ell)}^{2}}$ do not couple to the EM mode in any capacity
nor among themselves, and are therefore eigenstates of $\hat{H}_{2}$ in every
model discussed below, i.e.,

$\matrixelement{1_{\textrm{C}(k\ell)}^{2}}{\hat{H}_{2}}{1_{\textrm{C}(k^{\prime}\ell^{\prime})}^{2}}=2\omega_{10}\delta_{kk^{\prime}}\delta_{\ell\ell^{\prime}};$
(38)

as such, they will not be discussed in the next subsections.

Figure 1 illustrates the distribution of every spectral configuration among
the irreps. To be specific, the Hamiltonian matrix in the second excitation
manifold decomposes according to

$\mathbf{H}_{2}=\mathbf{H}_{2}^{(\textrm{A})}\oplus\quantity(\mathbf{1}_{N-1}\otimes\mathbf{H}_{2}^{(\textrm{B})})\oplus\quantity(\mathbf{1}_{N\quantity(N-3)/2}\otimes\mathbf{H}_{2}^{(\textrm{C})}),$
(39)

where the dimensions of $\mathbf{H}_{2}^{(\textrm{A})}$ and
$\mathbf{H}_{2}^{(\textrm{B})}$ depend upon the particular energetic structure
of the emitters, and $\mathbf{H}_{2}^{(\textrm{C})}=2\omega_{10}$.

Figure 1: Distribution of states in the doubly excited manifold among irreps.

###  III.1 Tavis-Cummings

If the emitters are well-approximated by two-level systems, the matrix
elements of the Hamiltonian for the doubly excited manifold are

$\displaystyle\bra{2_{0}}{\hat{H}_{2}}$
$\displaystyle\ket{2_{0}}=2\omega_{0},$ (40a)
$\displaystyle\bra{2_{0}}{\hat{H}_{2}}$
$\displaystyle\ket{1_{i}1_{j}}=\sqrt{2}g_{01}\delta_{i0},$ (40b)
$\displaystyle\bra{1_{i}1_{j}}{\hat{H}_{2}}$
$\displaystyle\ket{2_{0}}=\sqrt{2}g_{10}\delta_{i0},$ (40c) and
$\displaystyle\bra{1_{i}1_{j}}{\hat{H}_{2}}$
$\displaystyle\ket{1_{i^{\prime}}1_{j^{\prime}}}=\begin{cases}\omega_{0}+\omega_{10}&\text{if
}i=i^{\prime}=0,\enskip j=j^{\prime}\\\ 2\omega_{10}&\text{if
}i=i^{\prime}\neq 0,\enskip j=j^{\prime}\neq 0\\\ g_{01}&\text{if }i=0\neq
i^{\prime},\enskip j=j^{\prime}\\\ g_{10}&\text{if }i\neq 0=i^{\prime},\enskip
j=j^{\prime}\\\ 0&\text{otherwise}\end{cases}.$ (40d)

In turn, the matrix in the totally-symmetric subspace, i.e., with the
symmetrized basis
$\quantity{\ket{2_{0}},\ket{1_{0}1_{\textrm{A}}},\ket{1_{\textrm{A}}^{2}}}$,
is

$\mathbf{H}_{2,\textrm{TC}}^{(\textrm{A})}=\begin{pmatrix}2\omega_{0}&\sqrt{2N}g_{10}&0\\\
\sqrt{2N}g_{01}&\omega_{0}+\omega_{10}&\sqrt{2\quantity(N-1)}g_{10}\\\
0&\sqrt{2\quantity(N-1)}g_{01}&2\omega_{10}.\end{pmatrix}$ (41)

For states carrying the standard irrep,
$\quantity{\ket{1_{0}1_{\textrm{B}(k)}},\ket{1_{\textrm{B}(k)}^{2}}}$, the
corresponding Hamiltonian matrix is

$\mathbf{H}_{2,\textrm{TC}}^{(\textrm{B})}=\begin{pmatrix}\omega_{0}+\omega_{10}&\sqrt{N-2}g_{10}\\\
\sqrt{N-2}g_{01}&2\omega_{10}.\end{pmatrix}$ (42)

Notice that, as opposed to the singly excited manifold, the sub-space carrying
the standard irrep of the doubly excited Hamiltonian is no longer diagonal in
the symmetry-adapted basis.

Diagonalization of $\mathbf{H}_{2,\textrm{TC}}^{(\textrm{A})}$ gives rise to
the so-called bipolaritonic states[65, 38, 42, 49]

$\displaystyle{\hat{H}_{2}}$
$\displaystyle\ket{2_{\pm}}=\omega_{2\pm}\ket{2_{\pm}},$ (43a) and
$\displaystyle{\hat{H}_{2}}$
$\displaystyle\ket{1_{+}1_{-}}=\omega_{1_{+}1_{-}}\ket{1_{+}1_{-}}.$ (43b)

Eigenstates and eigenfrequencies of this $3\times 3$ matrix can be easily
obtained numerically. Analytical expressions, however, are cumbersome (see
appendix A) unless the reasonable approximation $N\gg 1$ is considered, in
which case:

$\displaystyle\ket{2_{\pm}}_{\textrm{TC}}=$ $\displaystyle
h_{\pm}^{2}\ket{2_{0}}\pm\sqrt{2}h_{+}h_{-}\ket{1_{0}1_{\textrm{A}}}+h_{\mp}^{2}\ket{1_{\textrm{A}}^{2}},$
(44) and $\displaystyle\ket{1_{+}1_{-}}_{\textrm{TC}}=$
$\displaystyle\sqrt{2}h_{+}h_{-}\quantity(-\ket{2_{0}}+\ket{1_{\textrm{A}}^{2}})$
$\displaystyle+\quantity(h_{+}^{2}-h_{-}^{2})\ket{1_{0}1_{\textrm{A}}}.$ (45)

The eigenvectors of $\mathbf{H}_{2,\textrm{TC}}^{(\textrm{B})}$ give the
coefficients for the oftentimes referred to as “upper(lower) polariton-dark
states”: [38, 42]

$\ket{1_{\pm}1_{\textrm{B}(k)}}_{\textrm{TC}}=\pm
h_{\pm}\ket{1_{0}1_{\textrm{B}(k)}}+h_{\mp}\ket{1_{\textrm{B}(k)}^{2}},$ (46)

which fulfill

$\hat{H}_{2}\ket{1_{\pm}1_{\textrm{B}(k)}}_{\textrm{TC}}=\omega_{1_{\pm}1_{\textrm{B}}}\ket{1_{\pm}1_{\textrm{B}(k)}}_{\textrm{TC}}.$
(47)

The eigenfrequencies are tabulated in table 2, while the matrix elements of
transition operators between the eigenstates in the first excitation manifold
and those in eqs. 44 to 46 are included in tables 3, 4, 5 and 6. The
analytical expressions for the solutions with arbitrary $N$ can be found in
Appendix A.

###  III.2 Harmonic limit

When two levels are not enough to depict the structure of the emitters, the
simplest way to add complexity to the system is by considering the emitters as
harmonic oscillators, i.e., with evenly-spaced levels in their energy
spectrum. This assumption is valid for modes related to bonds with high
dissociation energies. In this scenario, the new light-matter Hamiltonian
incorporates additional matrix elements given by

$\displaystyle\bra{2_{i>0}}{\hat{H}_{2}}$
$\displaystyle\ket{2_{j}}=2\omega_{10}\delta_{ij},$ (48a)
$\displaystyle\bra{1_{0}1_{i}}{\hat{H}_{2}}$
$\displaystyle\ket{2_{j}}=\sqrt{2}g_{01}\delta_{ij},$ (48b) and
$\displaystyle\bra{2_{i}}{\hat{H}_{2}}$
$\displaystyle\ket{1_{0}1_{j}}=\sqrt{2}g_{10}\delta_{ij}.$ (48c)

Notice that the harmonic approximation assumes $\omega_{12}=\omega_{10}$ and
$g_{12}=\sqrt{2}g_{01}$.

With the new symmetrized contributions to the space of functions, namely
$\ket{2_{\textrm{A}}}$ for $\mathbf{H}_{2}^{(\textrm{A})}$, and
$\quantity{\ket{2_{\textrm{B}(k)}}}$ for $\mathbf{H}_{2}^{(\textrm{B})}$, the
expanded Hamiltonian matrices are

$\mathbf{H}_{2,\textrm{HO}}^{(\textrm{A})}=\\\
\begin{pmatrix}2\omega_{0}&\sqrt{2N}g_{10}&0&0\\\
\sqrt{2N}g_{01}&\omega_{0}+\omega_{10}&\sqrt{2\quantity(N-1)}g_{10}&\sqrt{2}g_{10}\\\
0&\sqrt{2\quantity(N-1)}g_{01}&2\omega_{10}&0\\\
0&\sqrt{2}g_{01}&0&2\omega_{10}\end{pmatrix},$ (49)

for the totally-symmetric irrep, and

$\mathbf{H}_{2,\textrm{HO}}^{(\textrm{B})}=\begin{pmatrix}\omega_{0}+\omega_{10}&\sqrt{N-2}g_{10}&g_{21}\\\
\sqrt{N-2}g_{01}&2\omega_{10}&0\\\ g_{12}&0&2\omega_{10}\end{pmatrix},$ (50)

for the standard irrep.

While the separation of the Hamiltonian accomplished with this approach
significantly facilitates the numerical diagonalization, a simpler approach to
obtain analytical results can be achieved by noticing that a harmonic bi-
linear Hamiltonian can always be written as a sum of normal harmonic modes. In
other words, the Hamiltonian of the system becomes

$\begin{split}\hat{H}_{\textrm{HO}}=&\omega_{0}\hat{a}_{0}^{\dagger}\hat{a}_{0}+\omega_{10}\sum_{i=1}^{N}\hat{a}_{i}^{\dagger}\hat{a}_{i}\\\
&+\sum_{i=1}^{N}\quantity(g_{10}\hat{a}_{i}^{\dagger}\hat{a}_{0}+g_{01}\hat{a}_{0}^{\dagger}\hat{a}_{i}),\end{split}$
(51)

where $\hat{a}_{i}=\sum_{v=0}^{\infty}\sqrt{v+1}\hat{\sigma}_{i}^{(v)}$ is the
bosonic annihilation operator acting on the $i$th emitter. A set of
symmetrized operators can be defined such that the Hamiltonian becomes

$\begin{split}\hat{H}_{\textrm{HO}}=&\omega_{0}\hat{a}_{0}^{\dagger}\hat{a}_{0}+\omega_{10}\hat{a}_{\textrm{A}}^{\dagger}\hat{a}_{\textrm{A}}+\sqrt{N}\quantity(g_{10}\hat{a}_{\textrm{A}}^{\dagger}\hat{a}_{0}+g_{01}\hat{a}_{0}^{\dagger}\hat{a}_{\textrm{A}})\\\
&+\omega_{10}\sum_{k=2}^{N}\hat{a}_{\textrm{B}(k)}^{\dagger}\hat{a}_{\textrm{B}(k)},\end{split}$
(52)

with creation operators
$\hat{a}_{\textrm{A}}^{\dagger}=\sum_{v=0}^{\infty}\sqrt{v+1}\hat{J}_{+}^{(v)}/\sqrt{N}$,
and $a_{\textrm{B}(k)}^{\dagger}=\sum_{n=1}^{N}c_{n}^{(k)}a_{i}^{\dagger}$,
where the coefficients $c_{n}^{(k)}$ are the same as in eq. 6. Furthermore,
defining the polaritonic modes through the creation operators
$\hat{a}_{\pm}^{\dagger}=\pm
h_{\pm}\hat{a}_{0}^{\dagger}+h_{\mp}\hat{a}_{\textrm{A}}^{\dagger}$, allows to
rewrite the Hamiltonian as

$\hat{H}_{\textrm{HO}}=\omega_{+}\hat{a}_{+}^{\dagger}\hat{a}_{+}+\omega_{-}\hat{a}_{-}^{\dagger}\hat{a}_{-}+\omega_{10}\sum_{k=2}^{N}\hat{a}_{\textrm{B}(k)}^{\dagger}\hat{a}_{\textrm{B}(k)}.$
(53)

Since the modes defined by this Hamiltonian are independent, any combination
of creation operators applied to the groundstate will generate an eigenstate.
Consequently, for the totally-symmetric subspace,

$\displaystyle\ket{2_{\pm}}_{\textrm{HO}}=$
$\displaystyle\frac{\quantity(\hat{a}_{\pm}^{\dagger})^{2}}{\sqrt{2}}\ket{0}$
$\displaystyle=$ $\displaystyle
h_{\pm}^{2}\ket{2_{0}}\pm\sqrt{2}h_{+}h_{-}\ket{1_{0}1_{\textrm{A}}}$ (54a)
$\displaystyle+h_{\mp}^{2}\quantity(\frac{\sqrt{N-1}\ket{1_{\textrm{A}}^{2}}+\ket{2_{\textrm{A}}}}{\sqrt{N}}),$
$\displaystyle\ket{1_{+}1_{-}}_{\textrm{HO}}=$
$\displaystyle\hat{a}_{\pm}^{\dagger}\hat{a}_{\mp}^{\dagger}\ket{0}$
$\displaystyle=$
$\displaystyle-\sqrt{2}h_{+}h_{-}\ket{2_{0}}+\quantity(h_{+}^{2}-h_{-}^{2})\ket{1_{0}1_{\textrm{A}}}$
(54b)
$\displaystyle+\sqrt{2}h_{+}h_{-}\quantity(\frac{\sqrt{N-1}\ket{1_{\textrm{A}}^{2}}+\ket{2_{\textrm{A}}}}{\sqrt{N}})$
and $\displaystyle\ket{2_{D\textrm{A}}}_{\textrm{HO}}=$
$\displaystyle\frac{1}{\sqrt{N-1}}\sum_{k=2}^{N}\frac{\quantity(\hat{a}_{\textrm{B}(k)}^{\dagger})^{2}}{\sqrt{2}}\ket{0}$
$\displaystyle=$
$\displaystyle\frac{-\ket{1_{\textrm{A}}^{2}}+\sqrt{N-1}\ket{2_{\textrm{A}}}}{\sqrt{N}}.$
(54c)

The eigenstates carrying the standard representation are

$\displaystyle\ket{1_{\pm}1_{\textrm{B}(k)}}_{\textrm{HO}}=$
$\displaystyle\hat{a}_{\pm}^{\dagger}\hat{a}_{\textrm{B}(k)}\ket{0}=\hat{a}_{\textrm{B}(k)}^{\dagger}\hat{a}_{\pm}\ket{0}$
$\displaystyle=$ $\displaystyle\pm h_{\pm}\ket{1_{0}1_{\textrm{B}(k)}}$ (55a)
$\displaystyle+h_{\mp}\quantity(\frac{\sqrt{N-2}\ket{1_{\textrm{B}(k)}^{2}}+\sqrt{2}\ket{2_{\textrm{B}(k)}}}{\sqrt{N}})$
and $\displaystyle\ket{2_{D\textrm{B}(k)}}_{\textrm{HO}}=$
$\displaystyle\frac{-\sqrt{2}\ket{1_{\textrm{B}(k)}^{2}}+\sqrt{N-2}\ket{2_{\textrm{B}(k)}}}{\sqrt{N}}.$
(55b)

The states $\ket{2_{D\textrm{A}}}$ and $\ket{2_{D\textrm{B}(k)}}$ are purely
molecular and account for the addition of states of the form $\ket{2_{i}}$ to
the Hilbert space.

The emerging eigenvalues are on display in table 2, and the transition
intensities between first and second excitation manifolds for this model are
summarized in tables 3, 4, 5 and 6. Figure 2 compares the behavior of the
energy spectrum for the TC model and the harmonic approximation as functions
of the number of coupled emitters and detuning. As it can be seen, the TC
model and the harmonic approximation yield indistinguishable results as the
number of emitters increases. To explain this observation, we remark that the
interaction between the cavity and the $1\to 2$ transition in the emitters is
not extensive, as opposed to the $0\to 1$ transition. Therefore, although the
overall light-matter coupling in the multi-level system is stronger than in
the TC model, the contribution from the additional level gets diluted as the
ensemble grows. We emphasize again the emergence of states in the multi-level
case, $\ket{2_{D\textrm{A}}}$ and $\ket{2_{D\textrm{B}}}$, whose lack of
photonic character makes them impervious to detuning.

Figure 2: Energy spectrum of the doubly excited manifold in the Tavis-Cummings
(TC) and the harmonic oscillator (HO) models. Spectra of totally-symmetric
states (a) and states carrying the standard irrep (b) as a function of the
number of emitters, $N$, for a cavity in resonance with the $0\to 1$
vibrational transition. Spectra of totally-symmetric states (c) and states
carrying the standard irrep (d) as a function of the detuning between the
cavity and the $0\to 1$ transition, for $N=10^{9}$. All plots were calculated
with constant collective coupling amplitude of $7\omega_{10}/100$. Notice that
(b) starts at $\log(N)>0$ because there are no states carrying the standard
representation for $N=1$. Also, notice that the $2_{D\mathrm{A}}$ and
$2_{D\mathrm{B}}$ states are only defined for the HO model, as the TC model
does not afford states with two excitations in the same emitter. Finally, the
lack of photonic character for $2_{D\mathrm{A}}$ and $2_{D\mathrm{B}}$ implies
that they remain dispersionless as a function of detuning $\Delta$.

###  III.3 General (anharmonic) case.

More realistically, emitters deviate from the behavior of harmonic
oscillators. In particular, for molecular vibrations, there is a mechanical
anharmonicity reflecting that vibrational energy levels are not evenly spaced,
i.e.,

$\omega_{12}=\omega_{10}\left(1-\chi\right),$ (56)

where $\chi$ is the mechanical anharmonicity constant. Additionally, an
electrical anharmonicity constant, $\gamma$, can also be defined stemming from

$g_{12}=\sqrt{2}g_{01}\left(1+\gamma\right).$ (57)

These anharmonicities are generally present in electronic transtions as well.
Notice also that the TC model effectively corresponds to the anharmonic case
with $\chi=0$ and $\gamma=-1$.

In terms of these parameters, the matrix elements of the light-matter
Hamiltonian need to be updated to

$\displaystyle\bra{2_{i>0}}{\hat{H}_{2}}$
$\displaystyle\ket{2_{j}}=\quantity(\omega_{10}+\omega_{21})\delta_{ij},$
(58a) $\displaystyle\bra{1_{0}1_{i}}{\hat{H}_{2}}$
$\displaystyle\ket{2_{j}}=g_{12}\delta_{ij},$ (58b)
$\displaystyle\bra{2_{i}}{\hat{H}_{2}}$
$\displaystyle\ket{1_{0}1_{j}}=g_{21}\delta_{ij}.$ (58c)

In the symmetry-adapted basis they become

$\displaystyle\bra{2_{\textrm{A}}}{\hat{H}_{2}}$
$\displaystyle\ket{2_{\textrm{A}}}=\omega_{10}+\omega_{21},$ (59a)
$\displaystyle\bra{2_{\textrm{A}}}{\hat{H}_{2}}$
$\displaystyle\ket{1_{0}1_{\textrm{A}}}=g_{12},$ (59b) and
$\displaystyle\bra{1_{0}1_{\textrm{A}}}{\hat{H}_{2}}$
$\displaystyle\ket{2_{\textrm{A}}}=g_{21},$ (59c)

for the totally-symmetric irrep, and

$\displaystyle\bra{2_{\textrm{B}(k)}}{\hat{H}_{2}}$
$\displaystyle\ket{2_{\textrm{B}(k^{\prime})}}=\quantity(\omega_{10}+\omega_{21})\delta_{kk^{\prime}},$
(60a) $\displaystyle\bra{2_{\textrm{B}(k)}}{\hat{H}_{2}}$
$\displaystyle\ket{1_{0}1_{\textrm{B}(k^{\prime})}}=g_{12}\delta_{kk^{\prime}},$
(60b) and $\displaystyle\bra{1_{0}1_{\textrm{B}(k)}}{\hat{H}_{2}}$
$\displaystyle\ket{2_{\textrm{B}(k^{\prime})}}=g_{21}\delta_{kk^{\prime}},$
(60c)

for the standard irrep. The effects of anharmonicities can thus be introduced
as corrections to the harmonic Hamiltonian matrices in eqs. 49 and 50.
Explicitly, the full Hamiltonian matrices, carrying the irrep $\Gamma$, in the
second-excitation manifold read

$\mathbf{H}_{2}^{(\Gamma)}=\mathbf{H}_{2,\textrm{HO}}^{(\Gamma)}+{\mathbf{H}^{\prime}_{2}}^{(\Gamma)},$
(61)

where

${\mathbf{H}^{\prime}_{2}}^{(\textrm{A})}=\begin{pmatrix}0&0&0&0\\\
0&0&0&\gamma\sqrt{2}g_{10}\\\ 0&0&0&0\\\
0&\gamma\sqrt{2}g_{01}&0&-\chi\omega_{10}\end{pmatrix}$ (62)

and

${\mathbf{H}^{\prime}_{2}}^{(\textrm{B})}=\begin{pmatrix}0&0&\gamma\sqrt{2}g_{10}\\\
0&0&0\\\ \gamma\sqrt{2}g_{01}&0&-\chi\omega_{10}\end{pmatrix},$ (63)

are the anharmonic corrections with $\Gamma=\textrm{A}$ and B, respectively.
The interaction between levels according to their symmetry is summarized in
fig. 3.

Figure 3: Interaction scheme among levels in the singly and doubly excited
manifolds. States carrying the totally-symmetric (A) irreducible
representation are non-degenerate, while states with the standard (B) and C
irreps are $(N-1)$-fold and $N(N-3)/2$-fold degenerate, respectively. States
with labels $2_{\textrm{A}}$ and $2_{\textrm{B}}$ do not appear in the Tavis-
Cummings model.

To analyze the effect of anharmonicity in the energy spectra we consider,
without loss of generality, the emitters as Morse oscillators, i.e., the mass-
normalized normal-mode elongation, $x$, is subject to a potential energy of
the form[66]

$V(x)=\hbar\omega_{10}\frac{\quantity(1+\chi)^{2}}{2\chi}\quantity[1-\exp(-\sqrt{\frac{\omega_{10}\chi}{\hbar}}x)]^{2}.$
(64)

For this system, the electrical and mechanical anharmonicities are related
through[67]

$\gamma=\frac{1}{2}\sqrt{\frac{(1+\chi)(16-\chi^{2})}{2(2+\chi)}}-1.$ (65)

Figure 4 compares the exact energy spectra of totally-symmetric states as a
function of anharmonicity for several values of $N$ obtained through numerical
diagonalization of $\mathbf{H}_{2}^{(\textrm{A})}$.

Figure 4: Exact energy spectrum of doubly-excited totally-symmetric
eigenstates as a function of anharmonicity, $\chi$, for several numbers of
Morse emitters. Resonant case ($\omega_{0}=\omega_{10}$). The labeling
identifies the frequencies of the adiabatic eigenstates as a function of
$\chi/\omega_{10}$. Since the states $\ket{2_{+}}$ and $\ket{1_{+}1_{-}}$ are
well-separated in energy, they have been labeled to reflect that they are very
similar to the corresponding harmonic eigenstates. The $\ket{2_{X}}$ and
$\ket{2_{Y}}$ adiabatic eigenstates show avoided crossings resulting from
anharmonic mixing of the harmonic states $\ket{2_{D\textrm{A}}}$ and
$\ket{2_{-}}$. In all cases, $\ket{2_{X}}\approx\ket{2_{D\textrm{A}}}$ and
$\ket{2_{Y}}\approx\ket{2_{-}}$ for $\chi\to 0$, and
$\ket{2_{X}}\approx\ket{2_{-}}$ and $\ket{2_{Y}}\approx\ket{2_{D\textrm{A}}}$
for large $\chi$. The photon content is evaluated as the expectation value of
the number operator $\hat{a}_{0}^{\dagger}\hat{a}_{0}$. Notice that, for large
$N$, the avoided crossings giving rise to the $\ket{2_{X}}$ and $\ket{2_{Y}}$
states are essentially crossings between the very anharmonic state
$\ket{2_{D\mathrm{A}}}$ and the bipolariton $\ket{2_{-}\\}}$. These crossings
are loci of two-photon absorption enhancement.

As can be seen, the harmonic states $\ket{2_{-}}$ and $\ket{2_{D\textrm{A}}}$
interact anharmonically with a coupling strength determined by the single-
molecule light-matter coupling constant $g_{21}$. Therefore, it is possible to
find a set of parameters in the Hamiltonian for which these states are near-
resonant for large $N$, which is the limit that concerns us for collective SC.
Under this condition, the absorption of two photons with the frequency of the
$\ket{0}\to\ket{1_{-}}$ transition will experience an enhancement, as
discussed in refs. 40, 49. This phenomenon can be understood as follows: the
bipolariton $\ket{2_{-}}$ provides an optical window to funnel energy
efficiently into the anharmonic $\ket{2_{D\textrm{A}}}$ state.[48, 49] The
effect can be observed when $\omega_{2-}\approx 2\omega_{10}-\chi$. At light-
matter resonance with the fundamental transition ($\omega_{0}=\omega_{10}$),
this condition translates into the more intuitive $\Omega_{10}\approx\chi$,
i.e., the Rabi splitting must be tuned to to match the anharmonic shift. Note
that the bipolariton $\ket{2_{-}}$ is the only possible optical window to
resonate with $\ket{2_{D\textrm{A}}}$ when $\omega_{0}=\omega_{10}$. However,
if the photon is negatively detuned from the fundamental transition,
$\omega_{0}<\omega_{10}$, a second bipolariton $\ket{1_{+}1_{-}}$ can also be
in resonance with $\ket{2_{D\textrm{A}}}$, as shown in fig. 5. This effect was
not discussed in Ref. 49, and should provide an extra tuning parameter to
enhance multiphoton absorption processes in polaritonic systems. The
conditions for enhanced two-photon absorption (TPA) are better illustrated in
the correlation diagrams of fig. 6.

Figure 5: Analogous plot to fig. 4, panel $N=10^{6}$, except that the cavity
is negatively detuned ($\omega_{0}=0.95\omega_{10}$). The eigenstate
$\ket{2_{+}}$ is well-separated in energy and has been labeled to reflect that
it is very similar to the corresponding harmonic eigenstate. The
$\ket{2_{R}}$, $\ket{2_{S}}$, $\ket{2_{T}}$ adiabatic eigenstates show weakly
(because of large $N$) avoided crossings resulting from anharmonic mixing of
harmonic states $\ket{2_{D\textrm{A}}}$, $\ket{1_{+}1_{-}}$, and
$\ket{2_{-}}$. Three anharmonicity regimes can be recognized: small $\chi$
($\ket{2_{R}}\approx\ket{2_{D\textrm{A}}}$,
$\ket{2_{S}}\approx\ket{1_{+}1_{-}}$, and $\ket{2_{T}}\approx\ket{2_{-}}$),
intermediate $\chi$ ($\ket{2_{R}}\approx\ket{1_{+}1_{-}}$,
$\ket{2_{S}}\approx\ket{2_{D\textrm{A}}}$, and
$\ket{2_{T}}\approx\ket{2_{-}}$), and large $\chi$
($\ket{2_{R}}\approx\ket{1_{+}1_{-}}$, $\ket{2_{S}}\approx\ket{2_{-}}$, and
$\ket{2_{T}}\approx\ket{2_{D\textrm{A}}}$). The avoided crossings occur at
resonances between the very anharmonic $\ket{2_{D\textrm{A}}}$ with the
bipolaritons $\ket{1_{+}1_{-}}$ and $\ket{2_{-}}$, the latter of which
features a significant photonic content. These crossings are the loci of two-
photon absorption enhancement. Notice that, owing to the negative detuning,
there are two (rather than just one) such crossings, as opposed to the
resonant case shown in fig. 4. Figure 6: Correlation diagrams of polaritonic
states under the conditions that afford enhanced TPA. (a) If the photon mode
is resonant with the $0\to 1$ molecular transition (i.e.
$\omega_{0}=\omega_{10}$), enhancement of TPA occurs when
$\omega_{2-}=\omega_{10}+\omega_{21}$ due to anharmonic coupling between the
bipolariton $\ket{2_{-}}$ and the anharmonic state $\ket{2_{D\textrm{A}}}$
(red dotted line). (b) If the photon mode is negatively detuned from the $0\to
1$ molecular transition (i.e. $\omega_{0}<\omega_{10}$), enhancement of TPA
can also occur when $\omega_{1+1-}=\omega_{10}+\omega_{21}$ due to anharmonic
coupling between the bipolariton $\ket{1_{+}1_{-}}$ and the anharmonic state
$\ket{2_{D\textrm{A}}}$ (red dotted line).

The calculations above show that, as opposed to the case of linear response,
simulations with small $N$ are not a reliable representation of the nonlinear
optics of systems under collective SC for anharmonicities up to one order of
magnitude larger than the Rabi splitting, even if the intensity of the total
coupling, $\sqrt{N}g_{10}$, is fixed to match experimental values of
$\Omega_{10}$. Thus, we believe that the interpretation of experimental
spectra as reported in refs. 28, 43 using a single-molecule model is worth
revisiting.

While the numerical results above shed much light on the physics of the
problem for all ranges of anharmonicity, it is iluminating to carry out a
perturbation theory analysis for small values of $\chi$ to obtain closed
approximate expressions for eigenvalues and eigenvectors of $\hat{H}_{2}$. To
zeroth order, we have $\ket{\Psi}^{(0)}=\ket{\Psi}_{\textrm{HO}}$. The first-
order corrections to the eigenfrequencies of the HO limit are included in
table 2. The $O(\chi,\gamma)$ corrections to the doubly excited polaritonic
modes are computed to be

$\begin{pmatrix}\bra{2_{0}}\\\ \bra{1_{0}1_{\textrm{A}}}\\\
\bra{1_{\textrm{A}}^{2}}\\\
\bra{2_{\textrm{A}}}\end{pmatrix}\ket{2_{\pm}}^{(1)}=\frac{h_{\mp}^{2}}{2N\Omega_{10}}\begin{pmatrix}{{h_{+}}^{2}{h_{-}}^{2}\left[{\pm
3\chi{\omega_{10}}+2\left({{h_{\mp}}^{2}-5{h_{\pm}}^{2}}\right)\gamma{\Omega_{10}}}\right]}\\\
{\sqrt{2}{h_{+}}{h_{-}}\left[{\left({2{h_{\mp}}^{2}-{h_{\pm}}^{2}}\right)\chi{\omega_{10}}\pm
2\left({{h_{\mp}}^{2}+2{h_{\pm}}^{2}-6{h_{+}}^{2}{h_{-}}^{2}}\right)\gamma{\Omega_{10}}}\right]}\\\
{{h_{\mp}}^{2}\sqrt{\frac{{N-1}}{N}}\left[{\left({\frac{{{h_{-}}^{2}-{h_{+}}^{2}\pm
3{h_{+}}^{2}{h_{-}}^{2}}}{{{h_{\pm}}^{2}}}}\right)\chi{\omega_{10}}+2\left({2{h_{\pm}}^{2}-{h_{\mp}}^{2}-6{h_{+}}^{2}{h_{-}}^{2}}\right)\gamma{\Omega_{10}}}\right]}\\\
{\frac{1}{{\sqrt{N}}}\left\\{{\mp\frac{{N+{h_{\mp}}^{2}\left({{h_{\pm}}^{2}-{h_{\mp}}^{2}-3{h_{+}}^{2}{h_{-}}^{2}}\right)}}{{{h_{\pm}}^{2}}}\chi{\omega_{10}}+2\left[{N+2{h_{\mp}}^{2}\left({2{h_{\pm}}^{2}-{h_{\mp}}^{2}-6{h_{+}}^{2}{h_{-}}^{2}}\right)}\right]\gamma{\Omega_{10}}}\right\\}}\end{pmatrix}.$
(66)

For non-zero detuning, the corrections to the remaining totally-symmetric
eigenstates are

$\displaystyle\begin{pmatrix}\bra{2_{0}}\\\ \bra{1_{0}1_{\textrm{A}}}\\\
\bra{1_{\textrm{A}}^{2}}\\\
\bra{2_{\textrm{A}}}\end{pmatrix}\ket{1_{+}1_{-}}^{(1)}=$
$\displaystyle\frac{2h_{+}h_{-}}{N\Omega_{10}}\begin{pmatrix}{-\sqrt{2}{h_{+}}^{2}{h_{-}}^{2}\gamma{\Omega_{10}}}\\\
{{h_{+}}{h_{-}}\left[{\chi{\omega_{10}}-\left({{h_{+}}^{2}-{h_{-}}^{2}}\right)\gamma{\Omega_{10}}}\right]}\\\
{\sqrt{\frac{{2\left({N-1}\right)}}{N}}{h_{+}}^{2}{h_{-}}^{2}\left({\frac{{2}}{{{h_{+}}^{2}-{h_{-}}^{2}}}\chi{\omega_{10}}-3\gamma{\Omega_{10}}}\right)}\\\
{-\frac{1}{{\sqrt{2N}}}\left[{\frac{{N-4{h_{+}}^{2}{h_{-}}^{2}}}{{{h_{+}}^{2}-{h_{-}}^{2}}}\chi{\omega_{10}}+\left({N-6{h_{+}}^{2}{h_{-}}^{2}}\right)\gamma{\Omega_{10}}}\right]}\end{pmatrix}$
(67) and $\displaystyle\begin{pmatrix}\bra{2_{0}}\\\
\bra{1_{0}1_{\textrm{A}}}\\\ \bra{1_{\textrm{A}}^{2}}\\\
\bra{2_{\textrm{A}}}\end{pmatrix}\ket{2_{D\textrm{A}}}^{(1)}=$
$\displaystyle\frac{\sqrt{N-1}}{2N\Omega_{10}}\begin{pmatrix}{\frac{1}{{{h_{-}}^{2}-{h_{+}}^{2}}}\chi{\omega_{10}}}\\\
{\frac{{\sqrt{2}}}{{{h_{+}}{h_{-}}}}\chi{\omega_{10}}}\\\
{\sqrt{\frac{{N-1}}{N}}\left({\frac{{1-5{h_{-}}^{2}{h_{+}}^{2}}}{{{h_{+}}^{2}{h_{-}}^{2}\left({{h_{-}}^{2}-{h_{+}}^{2}}\right)}}\chi{\omega_{10}}-2\gamma{\Omega_{10}}}\right)}\\\
{\frac{1}{{\sqrt{N}}}\left({\frac{{1-5{h_{-}}^{2}{h_{+}}^{2}}}{{{h_{+}}^{2}{h_{-}}^{2}\left({{h_{-}}^{2}-{h_{+}}^{2}}\right)}}\chi{\omega_{10}}-2\gamma{\Omega_{10}}}\right)}\end{pmatrix}.$
(68)

Under resonant conditions, these two states are degenerate in the harmonic
limit. After lifting the degeneracy, to zeroth order, these states become

$\displaystyle\ket{1_{+}1_{-}}^{(0)}=$
$\displaystyle\frac{\sqrt{N-1}\ket{2_{0}}-\sqrt{N}\ket{1_{\textrm{A}}^{2}}}{\sqrt{2N-1}}$
(69) $\displaystyle\ket{2_{D\textrm{A}}}^{(0)}=$
$\displaystyle-\frac{\sqrt{N}\ket{2_{0}}+\sqrt{N-1}\ket{1_{\textrm{A}}^{2}}}{\sqrt{2N\quantity(2N-1)}}$
$\displaystyle+\sqrt{\frac{2N-1}{2N}}\ket{2_{\textrm{A}}}.$ (70)

After this procedure, there are no first-order corrections for the state
$\ket{1_{+}1_{-}}$; on the other hand,

$\begin{split}\ket{2_{D\textrm{A}}}^{(1)}=&\frac{\sqrt{2N-1}}{2N\Omega_{10}}\left(-\frac{\gamma\Omega_{10}}{\sqrt{2}}\ket{2_{0}}+\chi\omega_{10}\ket{1_{0}1_{\textrm{A}}}\right.\\\
&\left.-\gamma\Omega_{10}\frac{\sqrt{N-1}\ket{1_{\textrm{A}}^{2}}+\ket{2_{\textrm{A}}}}{\sqrt{2N}}\right).\end{split}$
(71)

For states carrying the standard irrep, the corrections of first-order in
anharmonicity are

$\displaystyle\begin{pmatrix}\bra{1_{0}1_{\textrm{B}(k)}}\\\
\bra{1_{\textrm{B}(k)}^{2}}\\\
\bra{2_{\textrm{B}(k)}}\end{pmatrix}\ket{1_{\pm}1_{\textrm{B}(k)}}^{(1)}=$
$\displaystyle\frac{\sqrt{2}h_{\mp}}{N\Omega_{10}}\begin{pmatrix}{\sqrt{2}{h_{+}}{h_{-}}\left[{\chi{\omega_{10}}-\left({{h_{+}}^{2}-{h_{-}}^{2}}\right)\gamma{\Omega_{10}}}\right]}\\\
{{h_{\mp}}^{2}\sqrt{\frac{{2\left({N-2}\right)}}{N}}\left[{\pm\left({\frac{{2{h_{\pm}}^{2}+{h_{\mp}}^{2}}}{{{h_{\pm}}^{2}}}}\right)\chi{\omega_{10}}-\left({3{h_{\pm}}^{2}+{h_{\mp}}^{2}}\right)\gamma{\Omega_{10}}}\right]}\\\
{\frac{1}{{\sqrt{N}}}\left\\{{\mp\left[{\frac{{N-2{h_{\mp}}^{2}\left({3{h_{\pm}}^{2}+{h_{\mp}}^{2}}\right)}}{{{h_{\pm}}^{2}}}}\right]\chi{\omega_{10}}+\left[{N-2{h_{\mp}}^{2}\left({3{h_{\pm}}^{2}+{h_{\mp}}^{2}}\right)}\right]\gamma{\Omega_{10}}}\right\\}}\end{pmatrix},$
(72) and $\displaystyle\begin{pmatrix}\bra{1_{0}1_{\textrm{B}(k)}}\\\
\bra{1_{\textrm{B}(k)}^{2}}\\\
\bra{2_{\textrm{B}(k)}}\end{pmatrix}\ket{2_{D\textrm{B}(k)}}^{(1)}=$
$\displaystyle\frac{\sqrt{2(N-2)}}{h_{+}h_{-}N\Omega_{10}}\begin{pmatrix}{\chi{\omega_{10}}}\\\
{-\sqrt{\frac{{N-2}}{N}}\left({\frac{{{h_{+}}^{2}-{h_{-}}^{2}}}{{{h_{+}}{h_{-}}}}\chi{\omega_{10}}+{h_{+}}{h_{-}}\gamma{\Omega_{10}}}\right)}\\\
{-\sqrt{\frac{2}{N}}\left({\frac{{{h_{+}}^{2}-{h_{-}}^{2}}}{{{h_{+}}{h_{-}}}}\chi{\omega_{10}}+{h_{+}}{h_{-}}\gamma{\Omega_{10}}}\right)}\end{pmatrix}.$
(73)

Table 2: Eigenfrequencies in the doubly excited manifold. The Tavis-Cummings (TC) model is effectively an anharmonic system with $\chi=0$ and $\gamma=-1$; therefore, when $N\gg 1$, the TC model yields the same results as the Harmonic limit. Eigenmode | Harmonic limit |  $O(\chi,\gamma)$ anharmonic corrections
---|---|---
$2_{\pm}$ | $2\omega_{\pm}$ | ${\displaystyle-\frac{h_{\mp}^{4}}{N}\quantity(\chi\omega_{10}\mp 4h_{\pm}^{2}\gamma\Omega_{10})}$
$1_{+}1_{-}$ | $\omega_{0}+\omega_{10}$ | ${\displaystyle-\frac{2h_{+}^{2}h_{-}^{2}}{N}\quantity[\quantity(1-\delta_{h_{+}h_{-}})\chi\omega_{10}-2\quantity(h_{+}^{2}-h_{-}^{2})\gamma\Omega_{10}]}$
$2_{D\textrm{A}}$ | $2\omega_{10}$ | ${\displaystyle-\frac{2N-2+\delta_{h_{+}h_{-}}}{2N}\chi\omega_{10}}$
$1_{\pm}1_{\textrm{B}}$ | $\omega_{10}+\omega_{\pm}$ | ${\displaystyle-\frac{2h_{\mp}^{2}}{N}\quantity(\chi\omega_{10}\mp 2h_{\pm}^{2}\gamma\Omega_{10})}$
$2_{D\textrm{B}}$ | $2\omega_{10}$ | ${\displaystyle-\frac{N-2}{N}\chi\omega_{10}}$
$1_{\textrm{C}}^{2}$ | $2\omega_{10}$ | 0
Table 3: Dipolar transition intensities between eigenstates carrying the totally symmetric irrep in the first and second excited manifolds. $\ket{\Psi}$ | $\matrixelement{\Psi}{\hat{\bm{\mu}}/\mu_{10}}{\Phi}$
---|---
| $\ket{\Phi}$
|  | $\ket{1_{+}}$ | $\ket{1_{-}}$
Tavis-Cummings
$\ket{2_{+}}$ |  | $\sqrt{2}h_{-}\quantity(h_{+}^{2}\sqrt{N}+h_{-}^{2}\sqrt{N-1})$ | $\sqrt{2}h_{+}h_{-}^{2}\quantity(\sqrt{N-1}-\sqrt{N})$
$\ket{2_{-}}$ |  | $\sqrt{2}h_{+}^{2}h_{-}\quantity(\sqrt{N-1}-\sqrt{N})$ | $\sqrt{2}h_{+}\quantity(h_{+}^{2}\sqrt{N-1}+\sqrt{2}h_{-}^{2})$
$\ket{1_{+}1_{-}}$ |  | $h_{+}\quantity[h_{+}^{2}\sqrt{N}+h_{-}^{2}\quantity(2\sqrt{N-1}-\sqrt{N})]$ | $h_{-}\quantity[h_{+}^{2}\quantity(2\sqrt{N-1}-\sqrt{N})+h_{-}^{2}\sqrt{N}]$
Harmonic approximation
$\ket{2_{+}}$ |  | $\sqrt{2N}h_{-}$ | 0
$\ket{2_{-}}$ |  | 0 | $\sqrt{2N}h_{+}$
$\ket{1_{+}1_{-}}$ |  | $\sqrt{N}h_{+}$ | $\sqrt{N}h_{-}$
$\ket{2_{D\textrm{A}}}$ |  | 0 | 0
$O(\chi,\gamma)$ anharmonic corrections
$\ket{2_{+}}$ |  | $\sqrt{\frac{2}{N}}h_{+}^{2}h_{-}^{3}\quantity[-\frac{\chi\omega_{10}}{\Omega_{10}}+\quantity(3h_{+}-h_{-}^{2})\gamma]$ | $-\frac{h_{+}h_{-}^{2}}{\sqrt{2N}}\quantity[\quantity(h_{+}^{2}+2h_{-}^{2})\frac{\chi\omega_{10}}{\Omega_{10}}-2\quantity(h_{+}^{2}-h_{-}^{2}-2h_{+}^{2}h_{-}^{2})\gamma]$
$\ket{2_{-}}$ |  | $-\frac{h_{+}^{2}h_{-}}{\sqrt{2N}}\quantity[\quantity(2h_{+}^{2}+h_{-}^{2})\frac{\chi\omega_{10}}{\Omega_{10}}-2\quantity(h_{+}^{2}-h_{-}^{2}-2h_{+}^{2}h_{-}^{2})\gamma]$ | $\sqrt{\frac{2}{N}}h_{+}^{3}h_{-}^{2}\quantity[\frac{\chi\omega_{10}}{\Omega_{10}}-\quantity(h_{+}-3h_{-}^{2})\gamma]$
$\ket{1_{+}1_{-}}$ | $\Delta\neq 0$ | $\frac{2h_{+}h_{-}^{4}}{\sqrt{N}}\quantity[\frac{\chi\omega_{10}}{\Omega_{10}}-\quantity(3h_{+}^{2}-h_{-}^{2})\gamma]$ | $\frac{2h_{+}^{4}h_{-}}{\sqrt{N}}\quantity[-\frac{\chi\omega_{10}}{\Omega_{10}}+\quantity(h_{+}^{2}-3h_{-}^{2})\gamma]$
| $\Delta=0$ | $-\sqrt{\frac{N\quantity(N-1)}{2N-1}}$ | $-\sqrt{\frac{N\quantity(N-1)}{2N-1}}$
$\ket{2_{D\textrm{A}}}$ | $\Delta\neq 0$ | $\sqrt{\frac{N-1}{2N}}h_{-}\quantity[\frac{2h_{+}^{2}-h_{-}^{2}}{h_{+}^{2}\quantity(h_{+}^{2}-h_{-}^{2})}\frac{\chi\omega_{10}}{\Omega_{10}}-2\gamma]$ | $\sqrt{\frac{N-1}{2N}}h_{+}\quantity[-\frac{h_{+}^{2}-2h_{-}^{2}}{h_{-}^{2}\quantity(h_{+}^{2}-h_{-}^{2})}\frac{\chi\omega_{10}}{\Omega_{10}}-2\gamma]$
| $\Delta=0$ | $\sqrt{\frac{N}{2\quantity(2N-1)}}+\frac{1}{2}\sqrt{\frac{2N-1}{2N}}\quantity(\frac{\chi\omega_{10}}{\Omega_{10}}-\gamma)$ | $\sqrt{\frac{N}{2\quantity(2N-1)}}-\frac{1}{2}\sqrt{\frac{2N-1}{2N}}\quantity(\frac{\chi\omega_{10}}{\Omega_{10}}+\gamma)$
Table 4: Dipolar transition intensities between eigenstates carrying the standard irrep in the first and second excited manifolds. Model | $\matrixelement{1_{\pm}1_{\textrm{B}}}{\hat{\bm{\mu}}/\mu_{10}}{1_{\textrm{B}}}$ | $\matrixelement{2_{D\textrm{B}}}{\hat{\bm{\mu}}/\mu_{10}}{1_{\textrm{B}}}$
---|---|---
Tavis-Cummings | $h_{\mp}\sqrt{N-2}$ | -
Harmonic approximation | $\sqrt{N}$ | 0
$O(\chi,\gamma)$ anharmonic corrections | $\mp\frac{2h_{\pm}h_{+}h_{-}}{\sqrt{N}}\quantity[\frac{\chi\omega_{10}}{\Omega_{10}}-\quantity(h_{+}^{2}-h_{-}^{2})\gamma]$ | $-\sqrt{\frac{2\quantity(N-2)}{N}}\quantity(\frac{h_{+}^{2}-h_{-}^{2}}{h_{+}^{2}h_{-}^{2}}\frac{\chi\omega_{10}}{\Omega_{10}}+\gamma)$
Table 5: Photon-induced transition intensities between eigenstates carrying the totally-symmetric irrep in the first and second excited manifolds. $\ket{\Psi}$ | $\matrixelement{\Psi}{\hat{a}_{0}}{\Phi}$
---|---
| $\ket{\Phi}$
|  | $\ket{1_{+}}$ | $\ket{1_{-}}$
Tavis-Cummings
$\ket{2_{+}}$ |  | $\sqrt{2}h_{+}$ | 0
$\ket{2_{-}}$ |  | 0 | $-\sqrt{2}h_{-}$
$\ket{1_{+}1_{-}}$ |  | $-h_{-}$ | $h_{+}$
Harmonic approximation
$\ket{2_{+}}$ |  | $\sqrt{2}h_{+}$ | 0
$\ket{2_{-}}$ |  | 0 | $-\sqrt{2}h_{-}$
$\ket{1_{+}1_{-}}$ |  | $-h_{-}$ | $h_{+}$
$\ket{2_{D\textrm{A}}}$ |  | 0 | 0
$O(\chi,\gamma)$ anharmonic corrections
$\ket{2_{+}}$ |  | $\frac{\sqrt{2}h_{+}h_{-}^{4}}{N}\quantity[\frac{\chi\omega_{10}}{\Omega_{10}}-\quantity(3h_{+}-h_{-}^{2})\gamma]$ | $\frac{h_{+}^{2}h_{-}^{3}}{\sqrt{2}N}\quantity(-\frac{\chi\omega_{10}}{\Omega_{10}}+4h_{+}^{2}\gamma)$
$\ket{2_{-}}$ |  | $-\frac{h_{+}^{3}h_{-}^{2}}{\sqrt{2}N}\quantity(\frac{\chi\omega_{10}}{\Omega_{10}}+4h_{-}^{2}\gamma)$ | $\frac{\sqrt{2}h_{+}^{4}h_{-}}{N}\quantity[\frac{\chi\omega_{10}}{\Omega_{10}}-\quantity(h_{+}-3h_{-}^{2})\gamma]$
$\ket{1_{+}1_{-}}$ | $\Delta\neq 0$ | $\frac{2h_{+}^{2}h_{-}^{3}}{{N}}\quantity[\frac{\chi\omega_{10}}{\Omega_{10}}-\quantity(3h_{+}^{2}-h_{-}^{2})\gamma]$ | $\frac{2h_{+}^{3}h_{-}^{2}}{{N}}\quantity[\frac{\chi\omega_{10}}{\Omega_{10}}-\quantity(h_{+}^{2}-3h_{-}^{2})\gamma]$
| $\Delta=0$ | $-\sqrt{\frac{N-1}{2N-1}}$ | $-\sqrt{\frac{N-1}{2N-1}}$
$\ket{2_{D\textrm{A}}}$ | $\Delta\neq 0$ | $-\frac{1}{N}\sqrt{\frac{N-1}{2}}\frac{h_{-}^{2}}{h_{+}\quantity(h_{+}^{2}-h_{-}^{2})}\frac{\chi\omega_{10}}{\Omega_{10}}$ | $\frac{1}{N}\sqrt{\frac{N-1}{2}}\frac{h_{+}^{2}}{h_{-}\quantity(h_{+}^{2}-h_{-}^{2})}\frac{\chi\omega_{10}}{\Omega_{10}}$
| $\Delta=0$ | $-\frac{1}{\sqrt{2\quantity(2N-1)}}+\frac{1}{2N}\sqrt{\frac{2N-1}{2}}\quantity(\frac{\chi\omega_{10}}{\Omega_{10}}-\gamma)$ | $\frac{1}{\sqrt{2\quantity(2N-1)}}+\frac{1}{2N}\sqrt{\frac{2N-1}{2}}\quantity(\frac{\chi\omega_{10}}{\Omega_{10}}+\gamma)$
Table 6: Photon-induced transition intensities between eigenstates carrying the standard irrep in the first and second excited manifolds. Model | $\matrixelement{1_{\pm}1_{\textrm{B}}}{\hat{a}_{0}}{1_{\textrm{B}}}$ | $\matrixelement{2_{D\textrm{B}}}{\hat{a}_{0}}{1_{\textrm{B}}}$
---|---|---
Tavis-Cummings | $\pm h_{\pm}$ | -
Harmonic approximation | $\pm h_{\pm}$ | 0
$O(\chi,\gamma)$ anharmonic corrections | $\mp\frac{2h_{\pm}h_{+}h_{-}}{\sqrt{N}}\quantity[\frac{\chi\omega_{10}}{\Omega_{10}}-\quantity(h_{+}^{2}-h_{-}^{2})\gamma]$ | $-\sqrt{\frac{2\quantity(N-2)}{N}}\quantity(\frac{h_{+}^{2}-h_{-}^{2}}{h_{+}^{2}h_{-}^{2}}\frac{\chi\omega_{10}}{\Omega_{10}}+\gamma)$

The main lesson from these perturbative correction calculations is that the
anharmonicity-induced energy shifts should be undetectable as they scale as
$1/N\to 0$ when $N\gg 1$, as is the case of collective SC. Hence, calculations
using one or small number of emitters cannot accurately model the nonlinear
response of vibrational polariton systems under those conditions.

##  IV Conclusions

In this work, we have found the exact eigenspectrum of the doubly excited
manifold of an ensemble of $N$ anharmonic oscillators under collective VSC. We
provided a group-theoretical formalism to solve this many-body problem,
reducing it to small matrix diagonalizations (the largest of which involves a
$4\times 4$ matrix). This procedure is a significant simplification from the
brute-force numerical approach, which involves diagonalization of
astronomically large ${N+2\choose 2}$-dimensional matrices where
$N=10^{6}-10^{10}$. We provide compact expressions and tables that should
serve as a concise reference for future work involving the nonlinear
spectroscopy of molecular polariton systems.

Through numerically exact examples and analytical studies using perturbation
theory, we have contrasted this model to the standard TC and HO models, and
have shown that the additionally available anharmonic transitions per molecule
can give rise to new phenomena, such as enhancement of two-photon absorption
cross-sections owing to new resonances between bipolariton states and
anharmonic two-quanta vibrational states. We have demonstrated that there are
at most two such bipolariton states that fullfill those resonant conditions.
These conclusions, while studied specifically for vibrational SC conditions,
should have analogues in other frequency ranges, such as in the UV-visible
regime, where anharmonic shifts in electronic transitions are ubiquitous.

We have provided with a summary of analytical expresions for eigenfrequencies
and matrix elements of typical transition operators. These might serve as a
reference guide for future theoretical and experimental explorations of the
doubly excited manifold. We highlight the finding that the deviation from the
harmonic behavior is negliglible for large $N$. This statement seems to
contradict reports of strong optical nonlinearities in systems under
collective SC.[68, 39] The resolution to this apparent conundrum lies in
recognizing the essential role of anharmonic dissipative processes in the
creation of reservoirs populated with dark states.[44] This effect is not
considered in the present formalism as it ignores dissipation, and focuses on
the regime where multiphoton absorption processes occur faster than relaxation
into dark states.

Finally, our work also shows that calculations involving few emitters and
higher excitation manifolds (i.e. nonlinear response calculations) cannot
provide a satisfactory approximation to the nonlinear optics of the collective
SC regime even if the single-molecule dipole is artificially increased to keep
the overall light-matter coupling constant. We hope that these findings inform
future experimental endeavors and the theoretical considerations regarding
polaritonic chemical behavior and non-linear response.

##  Appendix A Exact expressions for TC solutions.

Diagonalization of the TC Hamiltonian carrying the totally-symmetric irrep,
eq. 41, entails solving a cubic characteristic polynomial. Although it is
possible to obtain analytical expressions, this task is rather complicated and
a numerical approach can be easily implemented if exact solutions are needed.
Nonetheless, we list the eigenfrequencies for the corresponding bipolaritons
as a reference:

$\expectationvalue{\hat{H}_{2}}{\Psi}_{\textrm{TC}}=\omega_{0}+\omega_{10}+f_{\Psi}(\rho)\Omega_{10},$
(74)

where

$\displaystyle f_{2+}(\rho)=$ $\displaystyle 2\rho^{\prime},$ (75a)
$\displaystyle f_{2-}(\rho)=$
$\displaystyle-\rho^{\prime}-\sqrt{3}\rho^{\prime\prime},$ (75b) and
$\displaystyle f_{1+1-}(\rho)=$
$\displaystyle-\rho^{\prime}+\sqrt{3}\rho^{\prime\prime}.$ (75c)

In eq. 75, $\rho^{\prime}$ and $\rho^{\prime\prime}$ refer to the real and
imaginary parts, respectively, of the quantity $\rho$, which fulfills

$\rho^{3}=q^{3}+\textrm{i}\sqrt{p^{6}-q^{6}},$ (76a) where the coefficients
$p$ and $q$ are such that $\displaystyle
p^{2}=\frac{N-2h_{+}^{2}h_{-}^{2}}{3N}=\frac{1}{3}\quantity(1-\frac{2g_{10}^{2}}{\Omega_{10}^{2}}),$
(76b) and $\displaystyle
q^{3}=\frac{h_{+}^{2}h_{-}^{2}}{N}\quantity(h_{+}^{2}-h_{-}^{2})=\frac{g_{10}^{2}}{\Omega_{10}^{3}}\Delta.$
(76c)

Finally, the eigenvectors are given by

$\ket{\Psi}=\mathcal{N}_{\Psi}\begin{pmatrix}2\sqrt{N-1}h_{+}^{2}h_{-}^{2}\\\
\sqrt{2\quantity(N-1)}h_{+}h_{-}\varphi_{\Psi}(\rho)\\\
\sqrt{N}\quantity(f_{\Psi}(\rho)\varphi_{\Psi}(\rho)-2h_{+}^{2}h_{-}^{2})\end{pmatrix},$
(77)

where $\varphi_{\Psi}(\rho)=f_{\Psi}(\rho)-h_{+}^{2}+h_{-}^{2}$, and
$\mathcal{N}_{\Psi}$ is a normalization coefficient.

##  Appendix B Fourier basis with two excitations.

In this section, we discuss the dark states that emerge in the doubly excited
manifold beyond of what was covered on the main body of the manuscript. In the
harmonic limit, a transparent way to generate dark states with meaningful
labels is through the application of two creation operators related to the
eigenmodes. Explicitly, these states are

$\displaystyle\ket{2_{\gamma D(k)}}$
$\displaystyle=\frac{\hat{a}_{\textrm{B}(k)}^{\dagger}}{\sqrt{2}}\ket{1_{\textrm{B}(k)}},$
(78a) and $\displaystyle\ket{1_{\gamma D(k\ell)}^{2}}$
$\displaystyle=\hat{a}_{\textrm{B}(k)}^{\dagger}\ket{1_{\textrm{B}(\ell)}}=\hat{a}_{\textrm{B}(\ell)}^{\dagger}\ket{1_{\textrm{B}(k)}}.$
(78b)

The labels $k$ and $\ell$ clearly indicate either the emitters in which the
wavefunction is mostly localized, if the Schur-Weyl basis is used, or the wave
numbers assigned to the eigenfunction, in the case of the Fourier basis.
However, these functions no longer carry a defined irrep, but instead the
$N\quantity(N-1)/2$-dimensional reducible representation
$\gamma=\textrm{A}+\textrm{B}+\textrm{C}$ (not to be confused with the
electric anharmonicity parameter).

The generation of states carrying the standard irrep has been presented in
section III, as well as the explicit form of the C-symmetric states in the
Schur-Weyl basis, where

$\alpha_{mn}^{(k\ell)}=\alpha_{m}^{(k)}\alpha_{n}^{(\ell)}+\alpha_{n}^{(k)}\alpha_{m}^{(\ell)}-2\quantity(\alpha_{m}^{(k)}\delta_{n\ell}+\alpha_{n}^{(k)}\delta_{m\ell}).$
(79)

In the case of the Fourier basis, we illustrate the application of eq. 33 for
a scenario with four emitters. After removing the functions with symmetries A
and B, for $N=4$, we have $N\quantity(N-3)/2=2$, and it can be checked, with
the Schur-Weyl basis, that the allowed values for $k$ and $\ell$ are
$\\{2,3\\}$ and $4$, respectively. For $N=4$, we can define the wave vectors
proportional to $0\leq w\leq 3$. This fact gives rise to the combinations of
center of mass, $K$, and relative wave number, $q$, shown in table 7.

Table 7: Coefficients of center of mass and relative wave number for wavefunctions with $N=4$. | $w$ | $w^{\prime}$ | $K$ | $q$ |
---|---|---|---|---|---
| 0 | 0 | 0 | 0 |
| 0 | 1 | 1 | 1 |
| 0 | 2 | 2 | 2 |
| 0 | 3 | 3 | 3 |
| 1 | 1 | 2 | 0 |
| 1 | 2 | 3 | 1 |
| 1 | 3 | 4 | 2 |
| 2 | 2 | 4 | 0 |
| 2 | 3 | 5 | 1 |
| 3 | 3 | 6 | 0 |

Plugging these into eq. 33, we get

$c_{mn}^{(k,4)}=\psi(0,0)+\psi(2,0)\textrm{e}^{\frac{\pi
i}{2}\quantity(m+n)}+\psi(4,0)\textrm{e}^{\pi
i\quantity(m+n)}+\psi(6,0)\textrm{e}^{\frac{3\pi i}{2}\quantity(m+n)}\\\
+\quantity[\psi(1,1)\textrm{e}^{\frac{\pi
i}{4}\quantity(m+n)}+\psi(3,1)\textrm{e}^{\frac{3\pi
i}{4}\quantity(m+n)}+\psi(5,1)\textrm{e}^{\frac{5\pi
i}{4}\quantity(m+n)}]\cos\quantity(\pi\frac{m-n}{4})\\\
+\quantity[\psi(2,2)\textrm{e}^{\frac{\pi
i}{2}\quantity(m+n)}+\psi(4,2)\textrm{e}^{\pi
i\quantity(m+n)}]\cos\quantity(\pi\frac{m-n}{2})+\psi(3,3)\textrm{e}^{\frac{3\pi
i}{4}\quantity(m+n)}\cos\quantity(3\pi\frac{m-n}{4}).$ (80)

With help of Gram-Schmidt orthogonalization, it is possible to conclude that
all coefficients $\psi(K_{k,4},q_{k,4})$ vanish, except for $\psi(6,0)$,
$\psi(4,0)$, $\psi(4,2)$, and $\psi(2,2)$. With this information, we can write
write

$\displaystyle c_{mn}^{(k,4)}=$
$\displaystyle\frac{1}{2\sqrt{3}}\textrm{e}^{\pi
i\quantity(m+n)}\quantity[\cos\quantity(\pi\frac{m-n}{2})-1],$ (81a) and
$\displaystyle c_{mn}^{(k^{\prime},4)}=$
$\displaystyle-\frac{1}{2}\quantity[\textrm{e}^{\frac{3\pi
i}{2}\quantity(m+n)}+\textrm{e}^{\frac{\pi
i}{2}\quantity(m+n)}\cos\quantity(\pi\frac{m-n}{2})].$ (81b)

Identifying the appropriate values for $k$ and $k^{\prime}$ is rather
cumbersome and escapes the scope of this work.

###### Acknowledgements.

JACGA thanks Stephan van den Wildenberg and Matthew Du for their useful
comments and insights. The authors acknowledge funding support from the Air
Force Office of Scientific Research award FA9550-18-1-0289.

## References

  * [1] R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Campos-Gonzalez-Angulo, and J. Yuen-Zhou,  Chem. Sci. 9, 6325 (2018).
  * [2] J. Feist, J. Galego, and F. J. Garcia-Vidal,  ACS Photonics 5, 205 (2018).
  * [3] J. Flick, N. Rivera, and P. Narang,  Nanophotonics 7, 1479 (2018).
  * [4] J. Garcia-Vidal Francisco, C. Cristiano, and W. Ebbesen Thomas,  Science 373, eabd0336 (2021).
  * [5] D. N. Basov, A. Asenjo-Garcia, P. J. Schuck, X. Zhu, and A. Rubio,  Nanophotonics 10, 549 (2021).
  * [6] J. D. Plumhof, T. Stöferle, L. Mai, U. Scherf, and R. F. Mahrt,  Nature Materials 13, 247 (2014).
  * [7] J. Keeling and S. Kéna-Cohen,  Annu. Rev. Phys. Chem. 71, 435 (2020).
  * [8] M. A. Zeb, P. G. Kirton, and J. Keeling,  arXiv preprint arXiv:2004.09790 (2020).
  * [9] S. Pannir-Sivajothi, J. A. Campos-Gonzalez-Angulo, L. A. Martínez-Martínez, S. Sinha, and J. Yuen-Zhou,  arXiv preprint arXiv:2106.12156 (2021).
  * [10] K. Georgiou et al.,  ACS Photonics 5, 258 (2018).
  * [11] L. A. Martínez-Martínez, E. Eizner, S. Kéna-Cohen, and J. Yuen-Zhou,  The Journal of Chemical Physics 151, 054106 (2019).
  * [12] T. W. Ebbesen,  Acc. Chem. Res. 49, 2403 (2016).
  * [13] J. Yuen-Zhou and V. M. Menon,  Proceedings of the National Academy of Sciences 116, 5214 (2019).
  * [14] F. Herrera and J. Owrutsky,  J. Chem. Phys. 152, 100902 (2020).
  * [15] A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf,  Phys. Rev. A 69, 062320 (2004).
  * [16] J. Peng et al.,  Phys. Rev. Lett. 127, 043604 (2021).
  * [17] V. Agranovich,  Optika i Spektroskopiya 3, 84 (1957).
  * [18] J. J. Hopfield,  PR 112, 1555 (1958).
  * [19] C. Weisbuch, M. Nishioka, A. Ishikawa, and Y. Arakawa,  Phys. Rev. Lett. 69, 3314 (1992).
  * [20] D. G. Lidzey et al.,  Nature 395, 53 (1998).
  * [21] I. I. Rabi,  Phys. Rev. 51, 652 (1937).
  * [22] E. T. Jaynes and F. W. Cummings,  Proceedings of the IEEE 51, 89 (1963).
  * [23] R. H. Dicke,  Phys. Rev. 93, 99 (1954).
  * [24] M. Tavis and F. W. Cummings,  Phys. Rev. 170, 379 (1968).
  * [25] V. F. Crum, S. R. Casey, and J. R. Sparks,  Phys. Chem. Chem. Phys. 20, 850 (2018).
  * [26] K. Hirai, J. A. Hutchison, and H. Uji-i,  ChemPlusChem 85, 1981 (2020).
  * [27] J. D. Erwin, Y. Wang, R. C. Bradley, and J. V. Coe,  J. Phys. Chem. B 125, 8472 (2021).
  * [28] A. B. Grafton et al.,  Nature Communications 12, 214 (2021).
  * [29] K. Nagarajan, A. Thomas, and T. W. Ebbesen,  J. Am. Chem. Soc. 143, 16877 (2021).
  * [30] D. S. Wang and S. F. Yelin,  ACS Photonics 8, 2818 (2021).
  * [31] B. Xiang and W. Xiong,  The Journal of Chemical Physics 155, 050901 (2021).
  * [32] J. Yuen-Zhou, J. A. Campos-González-Angulo, R. F. Ribeiro, and M. Du,  Vibropolaritonic chemistry: theoretical perspectives,  in Metamaterials, Metadevices, and Metasystems 2021, edited by N. Engheta, M. A. Noginov, and N. I. Zheludev, volume 11795, pages 6 – 17, International Society for Optics and Photonics, SPIE, 2021.
  * [33] F. J. Hernández and F. Herrera,  The Journal of Chemical Physics 151, 144116 (2019).
  * [34] J. F. Triana, F. J. Hernández, and F. Herrera,  The Journal of Chemical Physics 152, 234111 (2020).
  * [35] E. W. Fischer and P. Saalfrank,  The Journal of Chemical Physics 154, 104311 (2021).
  * [36] C. Schäfer, J. Flick, E. Ronca, P. Narang, and A. Rubio,  arXiv preprint arXiv:2104.12429 (2021).
  * [37] D. S. Wang, T. Neuman, S. F. Yelin, and J. Flick,  arXiv preprint arXiv:2109.06631 (2021).
  * [38] N. Takemura et al.,  Phys. Rev. B 92, 125415 (2015).
  * [39] B. Xiang et al.,  Proceedings of the National Academy of Sciences 115, 4845 (2018).
  * [40] B. Xiang et al.,  J. Phys. Chem. A 123, 5918 (2019).
  * [41] T. M. Autry et al.,  Phys. Rev. Lett. 125, 067403 (2020).
  * [42] C. A. DelPo et al.,  J. Phys. Chem. Lett. 11, 2667 (2020).
  * [43] R. Duan, J. N. Mastron, Y. Song, and K. J. Kubarych,  J. Phys. Chem. Lett. 12, 11406 (2021).
  * [44] R. F. Ribeiro et al.,  J. Phys. Chem. Lett. 9, 3766 (2018).
  * [45] P. Saurabh and S. Mukamel,  J. Chem. Phys. 144, 124115 (2018).
  * [46] A. Debnath and A. Rubio,  Journal of Applied Physics 128, 113102 (2020).
  * [47] B. Gu and S. Mukamel,  J. Phys. Chem. Lett. 11, 8177 (2020).
  * [48] T. E. Li, A. Nitzan, and J. E. Subotnik,  The Journal of Chemical Physics 154, 094124 (2021).
  * [49] R. F. Ribeiro, J. A. Campos-Gonzalez-Angulo, N. C. Giebink, W. Xiong, and J. Yuen-Zhou,  Phys. Rev. A 103, 063111 (2021).
  * [50] F. T. Arecchi, E. Courtens, R. Gilmore, and H. Thomas,  Phys. Rev. A 6, 2211 (1972).
  * [51] R. Gilmore,  Ann. Phys. 74, 391 (1972).
  * [52] M. Gegg and M. Richter,  New J. Phys. 18, 043037 (2016).
  * [53] N. Shammah, S. Ahmed, N. Lambert, S. De Liberato, and F. Nori,  Phys. Rev. A 98, 063815 (2018).
  * [54] Y.-H. Lee, J. Links, and Y.-Z. Zhang,  Nonlinearity 24, 1975 (2011).
  * [55] T. Skrypnyk,  J. Math. Phys. 56, 023511 (2015).
  * [56] J. A. Campos-Gonzalez-Angulo, R. F. Ribeiro, and J. Yuen-Zhou,  23, 063081 (2021).
  * [57] M. A. Bastarrachea-Magnani, S. Lerma-Hernández, and J. G. Hirsch,  Phys. Rev. A 89, 032101 (2014).
  * [58] E. Choreño, D. Ojeda-Guillén, and V. D. Granados,  Journal of Mathematical Physics 59, 073506 (2018).
  * [59] A. Klimov and S. Chumakov,  A Group-Theoretical Approach to Quantum Optics: Models of Atom-Field Interactions,  Wiley, 2009.
  * [60] A. Strashko and J. Keeling,  Phys Rev A 94, 023843 (2016).
  * [61] M. Scully and M. Zubairy,  Quantum Optics,  Cambridge University Press, 1997.
  * [62] S. Cordero, O. Castaños, R. López-Peña, and E. Nahmad-Achar,  J. Phys. A: Math. Theor. 46, 505302 (2013).
  * [63] A. Bohr and B. Mottelson,  Nuclear Structure (In 2 Volumes),  World Scientific Publishing Company, 1998.
  * [64] H. Lipkin,  Lie Groups for Pedestrians,  Dover Books on Physics Series, Dover Publications, 2002.
  * [65] A. L. Ivanov and H. Haug,  Phys. Rev. Lett. 74, 438 (1995).
  * [66] P. M. Morse,  Phys. Rev. 34, 57 (1929).
  * [67] E. F. de Lima and J. E. M. Hornos,  J. Phys. B: At., Mol. Opt. Phys. 38, 815 (2005).
  * [68] A. D. Dunkelberger, B. T. Spann, K. P. Fears, B. S. Simpkins, and J. C. Owrutsky,  Nature Communications 7, 13504 (2016).