Generalization of the Tavis-Cummings model for multi-level anharmonic systems

The interaction between anharmonic quantum emitters (e.g., molecular vibrations) and confined electromagnetic fields gives rise to quantum states with optical and chemical properties that are different from those of their precursors. The exploration of these properties has been typically constrained to the first excitation manifold, the harmonic approximation, ensembles of two-level systems [Tavis-Cummings (TC) model], or the anharmonic single-molecule case. The present work studies, for the first time, a collective ensemble of identical multi-level anharmonic emitters and their dipolar interaction with a photonic cavity mode. The permutational properties of the system allow identifying symmetry classified submanifolds in the energy spectrum. Notably, in this approach, the number of particles, typically in the order of several millions, becomes only a parameter from the operational standpoint, and the size of the dimension of the matrices to diagonalize is independent of it. The formalism capabilities are illustrated by showing the energy spectrum structure, up to the third excitation manifold, and the calculation of the photon contents as a permutationally invariant quantity. Emphasis is placed on (a) the collective (superradiant) scalings of light-matter couplings and the various submanifolds of dark (subradiant) states with no counterpart in the single-molecule case, as well as (b) the delocalized modes containing more than one excitation per molecule with no equivalent in the TC model. We expect these findings to be applicable in the study of non-linear spectroscopy and chemistry of polaritons.


Introduction
Optical microcavities and similar devices enable the dipolar interaction between the electromagnetic (EM) field they confine and a suitable degree of freedom of quantum emitters [1][2][3][4][5]. The coupled system, whose excitations receive the name of polaritons, is a hybrid between light and matter, and its properties have motivated extensive exploration from disciplines such as quantum optics [6], excitonic and two-dimensional materials science [7][8][9], and chemistry [10,11]. In many systems of interest, observable effects emerge only when a large number of material dipoles cooperatively interact with a single photon mode. [12,13].
Polaritons have been produced using diverse setups and materials whose description requires various theoretical frameworks. For instance, real spins, NV centers [14], qubits [15], SQUIDS [16], quantum dots [17,18], and electronic transitions in organic molecules [19] are approched as two-level systems, while low vibrational transitions in molecules can be modeled with harmonic oscillators [20,21]. The success of these two models stems from their simplicity and the fact that conditions can be found for their realization in interesting experimental settings.
The interest in strongly coupled systems in which the emitters have a multi-level anharmonic spectrum has recently increased. For instance, the role that light-matter interaction might play on chemical reactivity has been explored by studying the Rabi model with a single emitter described by a Morse oscillator [22,23]. Additionally, signatures of the non-linearities present in states with two energy quanta have been observed experimentally [24][25][26][27], studied from theory [28,29], and motivated innovative theoretical ideas for signal enhancement [30,31].
The dynamics of the interaction between many multi-level systems and several EM modes has been discussed when non-energy-conserving terms of the Hamiltonian can be neglected, i.e., under the so-called rotating wave approximation (RWA) [32][33][34]. From these studies, it is known that analytical solutions exist for the coupling of r-level systems to r−1 EM modes [35]. In the context of the problem of N molecules coupled to a single EM mode, the RWA allows separating the Hamiltonian of the system allowing numerical approaches to the computation of exact solutions. However, in this kind of systems the Hilbert space scales as the number of molecules considered; therefore, these calculations are suitable only for systems of relatively small size. Nonetheless, for the case of identical emitters, the permutational symmetry implies a high degree of degeneracy, which can be exploited to reduce the number of degrees of freedom of the Hamiltonian, thus increasing the computational capabilities of the model. Such an approach has been utilized to provide a formulation for open quantum systems that can be used to explore spectroscopic measurements and other processes involving relaxation [36,37]. There is, however, little to no word about the stationary features of such systems, which are discussed in depth in the present work.
Permutational symmetry has been exploited in the few body limit to study systems such as nuclear structure [38], quantum circuits [39], magnons [40], ultracold atoms [41], nuclear spins within molecules [42], and the one-dimensional Hubbard model [43]. In the general case, the properties of the symmetrized quantum states have also been extensively discussed [44,45]. However, to the best of our knowledge, there has not been an effort to present concrete results from the application of this algebraic approach to the problem of light-matter coupling in the collective regime beyond the case of two-level systems, which is widely known [44,46,47].
In the present work, we provide a fathomable, insightful and easily implementable procedure to simplify the time-independent Schrödinger equation for N emitters, each with r non-degenerate and non-equispaced bound states in their energy spectrum, coupled to a harmonic field mode through a bilinear and excitation-conserving interaction. We identify a correlation between the symmetries of the system and the distribution of photonic component among the emerging manifolds. Of particular interest is the realization that collective states comprising excitations of distinct natures give rise to cooperative couplings that depart from the well-known factor of √ N, which characterizes polaritons in the singly excited manifold. We also present the effects of anharmonicity, detuning, and intensity of collective coupling on the eigenenergies, and photon contents of the eigenstates. This paper is organized as follows: in section 2 we introduce the model Hamiltonian, discuss the implied approximations, and provide a layout of the algebraic structure that enables its separation. In section 3 we introduce the tools and concepts that enable the block-diagonalization, and show the form and distribution of the diagonal elements. In section 4 we formulate the off-diagonal matrix elements that result from the blockdiagonalization. In section 5 we summarize the simplification method, providing a general recipe to put it in practice. In section 6 we present examples for the separation for the lowest-lying excitation manifolds. In section 7 we explore the implications of varying parameters on properties of the eigenstates. Finally, we present the conclusions in section 8.

Description of the model
The HamiltonianĤ describes a collection of N identical emitters with r+1 states each (which could represent the lowest vibrational bound states of the ground electronic state of an anharmonic molecule) and an EM mode with frequency ω. Here,â ( †) 0 is the annihilation (creation) operator of the EM mode, ε v is the energy of the vth state of the emitter, andσ (i) v,u is the transition operator between levels u and v in the ith emitter, i.e., In (2), the eigenstates ofĤ bare (N) in (1), henceforth primitive bare-eigenstates (PBEs), are written in a collective Fock representation where the index v i indicates the number of bosonic excitations in the ith emitter or the EM mode (i = 0). Hereafter, modes with v = 0 are considered implicitly for brevity. The Hamiltonian describing the emitters, the EM mode, and their dipolar interactions under the RWA iŝ with collective ladder operators for the emitters, Here, the coupling constant g v±u,v = ω 2ǫ 0 V v ± u|ex|v is the product of the amplitude of the single-photon electric field, confined to a mode volume V, and the value of the transition dipole moment between emitter states v ± u and v. Notice, in (4), that the model only accounts for coupling between states that differ by one excitation, and therefore neglects overtone transitions. The v,v acts on the PBEs according tô thus indicating the total number of excitations of a given state; moreover, since Ĥ (N),n exc = 0, this operator corresponds to a constant of motion. Therefore, the Hamiltonian can be recast in the form in which each of the terms in the sum of the right-hand-side can be represented as block matrices spanned by PBEs of constant n exc , these blocks generate the so-called excitation manifolds. Explicitly, we havê . . .
where |0 = |0 0 0 1 . . . 0 N is the global ground state and H.c. stands for Hermitian conjugate. While (8a) represents a significant simplification that allows diagonalization of (3), it is important to remark that the matrices generated with this approach can quickly become intractable as the number of molecules increases. For instance,Ĥ nexc (N) is an operator in N +nexc nexc dimensions. Therefore, it is impractical to deal with the number of molecules required to achieve observable ensemble strong coupling, which is usually in the order of millions [48][49][50].

Permutational symmetry
Further simplification of (3) requires to identify the additional symmetries of the Hamiltonian. An inspection of (8a) reveals that PBEs with the same distribution of quanta are degenerate under the action ofĤ bare (N); therefore, it becomes convenient to introduce a label that characterizes PBEs yet avoids spurious identification due to accidental degeneracies. In principle, this characterization could be achieved with the bare energy and the collection of eigenvalues of the operatorŝ which commute withĤ bare (N); indeed, the relevance of these operators and those in (5) will be revisited in section 4. Nonetheless, a much simpler approach is the description of the distribution of quanta itself. Let us define the spectral configuration: a notation device that indicates the number of emitters nμ(v) in the vth excited state. See table 1 for usage examples. The spectral configuration labels not only PBEs, but also any bare eigenstate (BE) resulting from linear combinations of PBEs with the same distribution of quanta. The BEs can thus be partially characterized with three labels: the excitation manifold (n exc ), the spectral configuration (μ), and a basis-dependent degeneracy index (y). These states fulfillĤ bare (N)|n exc ,μ; y = ε (μ) nexc |n exc ,μ; y , where ε (μ) nexc = ω v (nexc,μ) 0 is the characteristic bare energy for states with those labels, and v (nexc,μ) 0 is the corresponding number of quanta in the EM mode. Since r v=0 nμ(v) = N for a givenμ, the spectral configuration corresponds to a partition of the total number of emitters. For reasons that will become apparent, it is convenient to write these partitions in their regular (non-increasing) form, i.e., where µ i = nμ(v i ) under the constraint that µ i ≤ µ i ′ for i > i ′ [51]. In what follows, all the empty elements, µ i = 0, will be omitted for brevity, as illustrated in table 1. Furthermore, the spectral configurations can be identified as partitions of the number of excitations distributed among the emitters, i.e. for everyμ there exists a regular partition ν(μ) such that as illustrated in table 1. Because of the latter, the number of possible spectral configurations within a given manifold is where p(n) is the number of partitions of the integer n. These numbers can be extracted from the expansion [52] m k=1 Since the emitters are considered as identical,Ĥ(N) is invariant under the action of the symmetric group of degree N, S N , whose elements are the N! possible permutations of labels of the emitters. In other words, for all permutationsπ ∈ S N , Ĥ (N),π = 0.
The permutation operators act on the BEs according tô π|n exc ,μ; y = |n exc ,μ; y ′ , which means that BEs with common n exc andμ can be used to form a representation of S N . For instance,π are the six elements of S 3 in the basis of PBEs withμ = 1 1 . The representations of S N identified with a partition µ are, in general, reducible; they receive the name of permutation modules and are denoted by M µ [53]. The Table 2: Spectral configurations, EM mode excitations and degeneracies for the first five excitation manifolds.
dimension of this representation is the number of BEs with the same µ, i.e., The permutation modules can, in turn, be decomposed in irreducible representations (irreps) according to Young's rule [54]: where S λ symbolizes the irreps, also known as Specht modules [55]. The direct sum in (20) runs over all partitions λ of N that dominate (denoted by the symbol ) over µ, i.e., they fulfill [54] For instance, consider the spectral configurations in the triply excited manifold: 0 N , The coefficients K λµ in (20) are known as Kostka numbers [56]; they indicate the number of times a permutation module contains a Specht module. While obtaining a closed analytical expression to calculate them remains an open problem [57,58], these coefficients can be found through their combinatorial interpretation: the number of semistandard Young tableaux (SSYT) of shape λ and content µ [54]. A Young diagram (YD) of size N and shape λ is a collection of N cells arranged in r left-justified rows with λ i cells on the ith row. The present work uses the English notation, which is consistent with the regular form of the partitions (see table 1 for selected examples). A SSYT of shape λ and content µ is obtained by filling in the cells of a λ-shaped YD with a collection of ordered symbols partitioned according to µ in such a manner that the rows do not decrease to the right and the columns increase to the bottom [54].
The dimension of the Specht modules corresponds to the number of standard Young tableaux (SYT) of shape λ, i.e., the number of ways in which the sequence [1, 2, . . . , N] fills a λ-shaped YD such that the entries increase across rows and columns [54].
This quantity is given by the hook-length formula [59]: where h λ (i) represents the number of cells in the hook of the ith cell in the λ-shaped YD, i.e. the number of cells that are either directly below of, directly to the right of, or the ith cell itself. For example, for λ = [6, 4, 3, 2, 1], the hook of the cell with coordinates (2, 2) is , and the hook-lengths for each cell are 10 8 6 4 2 1 7 5 3 1 The dimensions of the permutation and Specht modules for the spectral configurations found up to n exc = 4 are displayed in table 2. The irreps identify the smallest possible subspaces that remain excluded under permutations. In other words, the space of BEs with the sameμ can be split into sets of symmetry adapted linear combinations of BEs (SABEs) labeled by λ for whicĥ π|n exc ,μ, λ; y, w = |n exc ,μ, λ; y ′ , w , for allπ ∈ S n , where the label w has been added to acknowledge repetition of irreps, i.e., when K λµ > 1. It is a well-known result from Representation Theory that where C r+1 is the vector space spanned by the energy eigenstates of each molecule, and therefore (C r+1 ) ⊗N is the vector space spanned by the BEs. The symbol ⊢ reads as partition of, and V λ r+1 is a so-called weight-space; it corresponds to an irrep of the unitary group U(r + 1). This result, known as the Schur-Weyl duality [60,61], implies that the SABEs are arranged in exclusive subspaces, labeled by λ, in which they are classified according to not only their behavior under permutations, but also under unitary operators. Furthermore, the decomposition in (29) provides with the meaning of all the indices in (28).
As previously discussed, the dimension of a Specht module corresponds with the number of SYT of shape λ, and the index y was used to enumerate them; therefore the indices λ and y uniquely define a SYT. In what remains of this paper, the index y will encode the elements of the Young tableaux after removal of the top row. See some examples in table 3. On the other hand, the dimension of the irreps of U(r + 1) correspond to the Kostka numbers. To be specific, the elements of the weight space V λ r+1 can be enumerated with SSYT, and thus uniquely identified with the indicesμ, λ, and w. In the remaining of the present work, w will encode the elements of the SSYT after removal of the top row. See some examples in table 4. Figure 1: Diagram of the relations between SABEs in the doubly excited manifold. GS denotes Gram-Schmidt orthogonalization, 2 ≤ k ≤ N, 2 ≤ k ′ < ℓ, and 4 ≤ ℓ ≤ N.
To gain some insight into the meaning of the irreps, let us consider the global ground-state of the emitters, i.e., the state withμ = 0 N . The only partition that dominates over [N] is itself, and there is a unique SSYT when the shape and content correspond to the same partition; therefore, M To exemplify how these concepts unfold, let us return to the representation in (18). Starting from the state |1, 0 3 , [3]; 0, 0 , the action of the operatorâ 0Ĵ The remaining two SABEs withμ = 0 2 1 1 are and The permutation operators in this basis are where the notation was simplified by making implicit the common indices n exc = 1 and µ = 0 2 1 1 , as well as by omitting the label w. As it can be seen in (33), the permutation operators have two clear independent subspaces labeled by [3] and [2,1], respectively. The spectral configurations,μ, giving rise to partitions that are valid labels, λ, for the irreps fulfill this constraint produces the so-called partitions with weakly decreasing multiplicities [62]. The enumeration of these partitions imply that the total number of irreps in the n exc th manifold is where q(n) is the number of partitions of n with weakly decreasing multiplicities as illustrated in table 2. These numbers can be extracted from the expansion [51,52] m k=2 The main result of this section is that the Hamiltonians of the excitation manifolds can be split according toĤ whereĤ with eachĤ (λ,y) nexc encoding the energies and interactions of spectral configurations with the same symmetry. Equation (38) implies that the states obtained from diagonalizingĤ (λ) nexc are dim(S λ )-fold degenerate. In particular, for the diagonal part of the Hamiltonian, we havê Notice that dim Ĥ (λ,y) The simplification achieved with this strategy is significant since, for each excitation manifold, the problem has been reduced to the diagonalization of |λ nexc | ∼ exp(An [63], both independent of N, which is a substantial gain over the original N +nexc nexc -dimensional matrices.
A complete discussion of the structure of the Hamiltonians requires to include the couplings between spectral configurations induced by collective excitations; these will be explored in section 4.

Schur-Weyl basis.
The SABEs obtained through application of the collective excitations,Ĵ (u,v) + , and Gram-Schmidt orthogonalization have been discussed in the literature under the name of Schur-Weyl states [40,[64][65][66], or Gelfand-Tsetlin states [67][68][69][70]. Although not necessary for the block-diagonalization of the Hamiltonian, an illustration of the explicit form of the symmetrized states in terms of the PBEs might be useful for the reader. Since symmetrization only affects the portion of the Hilbert space concerning the emitters, the quantum number n exc , as well as the operators acting on the photonic mode,â ( †) 0 , are not included in the following discussion, i.e., the SABEs are denoted by |μ, λ; y, w .
If the explicit form of the SABE in terms of PBEs is known, the application of collective excitation operators gives a straightforward result. For instance, the groundstate and some of its collective excitations are where λ = [N] and y = 0 indicate the SYT 1 2 . . . N . This tableau implies that the wavefunction must be invariant under any permutation of all the labels. Since all emitters are in the same state in the collective ground-state, this condition is met by default. The linear combinations of the excited spectral configurations are thus totally symmetric.
The states with spectral configuration such that µ = λ cannot be generated through excitation operators, and thus require a orthogonalization strategy. Since these states are highly degenerate, the change of basis is not unique; however, there is an approach that highlights the meaning of the label y. For a 1D array (i 1 , i 2 , . . . , i n ), let's define the vandermonde matrixV(i 1 , i 2 , . . . , i n ) with elements A given SYT can be associated to Young operators of the form where col j (A) extracts the jth column of array A, and y ′ ψ is any permutation of the elements in each row of y that produces downwards increasing columns. For instance, the SYT generates the arrays And the corresponding operator iŝ If µ = λ, the wavefunction |μ, λ; y, w is proportional toŶ (λ, y)|0 N , [N]; 0, 0 . Notice that the constraint µ = λ uniquely defines a SSYT, and therefore a w.
The states orthogonal to |0 N −1 1 1 , [N]; 0, 0 with the sameμ can be written as which have associated Young operators of the form Therefore, the corresponding SABEs are The functions |0 N −1 v 1 , [N − 1, 1]; k, v 1 for v > 1 have an equivalent form. Furthermore, the action of the collective excitationĴ Following the same reasoning, the remaining states withμ = 0 N −2 1 2 that are orthogonal to all the states above are The operators that act on the ground-state to generate all the SABEs relevant up to the triply excited manifold are tabulated in Appendix A.

Collective couplings.
After determining the SABEs, it is possible to evaluate the matrix elements of the interaction Hamiltonian. Let's consider the collective transition of a quantum from level v * to v * + s, i.e., the action of the operatorĴ v * +s,v * + on a SABE with spectral configurationμ. As a result of the transition, the new SABE belongs to the spectral configurationμ ′ , which relates to the initial spectral configuration according to Additionally, the occupation of levels in the two spectral configurations fulfill Given the definition ofĤ int (N) in (4), the only possible transitions are those between neighboring levels; consequentially, s = ±1. It can be shown that under this constraint, for a given manifold, the number of pairings {μ ′ ,μ} afforded by photonic (de)excitations is The matrix elements of the interaction Hamiltonian are where, in agreement with the discussion above, while the coefficient B accounts for the (de)excitation of the EM mode.
To describe the redistribution of quanta among the emitters, let's define the operatorsĴ where r ∈ {−, 0, +}, that generate the su(2) algebra with Casimir operator These operators allow the definition of the quantities and J (λ) where the explicit labels in bras and kets are the only relevant ones in the calculation of the indicated matrix elements as pointed out in the study of two-level systems, where the SABEs are known as Dicke states, and J (λ) µ ′ ,μ as the Dicke cooperation number [71]. To be specific and J (λ) According to the prescription of quantum angular momentum, the contribution of collective (de)excitations to the coupling between states is given by Couplings are said to be superradiant if L µ ′ ,μ < 1 [71]. Appendix B shows a couple of neat ways to calculate these coefficients.
The last component of (54), the coefficient C two strings of the same raising operators applied in different order to a SABE will not, in general, produce the same state. This situation creates an ambiguity whenĤ (λ,y) nexc requires more than one function with the sameμ, i,e, when K λµ > 1; otherwise, C and apply to it consecutive excitations to get to a state withμ = 0 N −2 1 1 2 1 . There are two possible paths: and While the identifications basis spanned by the previous states, then the coefficients can be extracted from the identity L (λ) For instance, let's set and and compute their orthogonal complements: and whereÎ is the identity operator. The values of the coefficients calculated with both working basis are shown in table 5. Since the bare energies, ε nexc , are independent of the index w, and

Algorithm for Hamiltonian separation
In this section we lay out the specific steps to write the block diagonalized Hamiltonian of a given excitation manifold.
(i) Define n exc , the manifold.
(ii) List all relevantμ, the spectral configurations.
(v) Determine the composition of eachμ in terms of λ according to Young's rule in (20).
(a) Find out the dominance relations among partitions according to (21). (b) Compute the pertinent Kostka numbers (K λµ ).
(vi) Write the diagonal ofĤ register the blocks of the Hamiltonian.
(a) Among the available spectral configurations, determine all pairings {μ ′ ,μ} such that The number of pairings inside a given manifold must be n int (n exc ) as calculated with (53) 3. Calculate the coupling matrix element according to (54).
(viii) Collect the calculated couplings in their respective blocks.

The triply excited manifold
In this section we illustrate the implementation of the algorithm when dealing with the system when it holds three quanta, i.e., n exc = 3. Table 6 compiles the quantities computed according to steps 1 − 4 of the algorithm in the previous section, as well as the dimensions of the irreps in step 6.b.
In table 7, we show the Kostka numbers that indicate the composition of the spaces [N ] spanned by wavefunctions with the sameμ in terms of symmetrized subspaces labeled by λ. Notice that K λµ = 0 implies that λ µ.
The factors required to calculate the off-diagonal elements, as prescribed by step 7, up until 7.c.i, are on display in table 8. Step 7.c.ii was worked in detail in the discussion leading to table 5. What is left to this section are the explicit forms of the overlaps: With all these elements, the Hopfield-Bogoliubov forms of the block-diagonalized Hamiltonian are and where 1 n is the n × n identity matrix. The matrices above are displayed in a diagrammatic form such that the horizontal lines represent energy levels, and therefore their labels correspond to diagonal elements, while the double-headed arrows indicate the couplings with their labels corresponding to off-diagonal matrix elements. In constructing (75b), the couplings involving the states withμ = 0 N −2 1 1 2 1 were worked out with the basis {|A , |B }.
In the case of the totally symmetric subspace, we can recognize a ladder of four superradiantly coupled levels, corresponding to the Dicke states. These are weakly connected to a two-level superradiant interaction between the SABEs with one excitation in v = 2. In turn, this interaction connects weakly with the remaining SABE in which the third excited state is populated. For the N − 1-degenerate subspace with λ = [N − 1, 1], the superradiantly connected Dicke states form a three-level system in Ξ configuration. They weakly couple to a Λ three-level system produced by the degeneracy introduced by the two-fold contribution ofμ = 0 N −2 1 1 2 1 to this symmetry. As in the case above, the remaining SABE interacts only weakly with this array.

Matrices for lower manifolds
Finally, we present the Hopfield-Bogoliubov form of the Hamiltonian operators for 0 ≤ n exc ≤ 2. Their properties have been exhaustively discussed elsewhere. (78c)

Properties of eigenstates.
To showcase the usefulness of the formalism, in this section, we analyze the behavior of the energy spectrum and photon content of the eigenstates as a function of parameters in the Hamiltonian, such as anharmonicity, intensity of coupling and detuning. If each of the emitters is a harmonic oscillator (HO), we can definê and the eigenstates can be represented as excitations of the polariton modes: The quantities ∆ = ω − (ε 1 − ε 0 )/ , and Ω 10 = [∆ 2 + 4(g 10 / ) 2 N] 1/2 are the detuning and Rabi frequency, respectively, and the labels UP, LP and D stand for upper, lower and dark polaritons, respectively. When anharmonicity is considered, the excitations in the polaritonic modes are no longer good quantum numbers; however, the actual eigenstates are similar enough to those in the harmonic case that the labels introduced in (80) can still be consistent with some of the features displayed by these states. We introduce anharmonicities by considering each emitter as a Morse oscillator (MO) [72], i.e., the single-emitter Hamiltonians include a potential energy function of the form where R is the mass-scaled length of the oscillator with value at equilibrium R e , is the dissociation energy, and The corresponding eigenenergies are given by with the implication that the number of bound-states is ⌊r MO ⌋, where therefore, the potential becomes harmonic as D e → ∞ (a → 0) [73]. From (83), it follows that, in the neighborhood of the equilibrium length, the emitters oscillate with frequency On the other hand, the difference in energies of contiguous levels fulfill therefore, this model introduces a mechanical anharmonicity characterized by −a 2 . There is also an electric anharmonicity, which stems from the fact that [74] It is pertinent to remark that the dipole moment operator induces transitions between all the levels in a MO; however, the model in this work focuses only on the v → v ± 1 processes. Figure 2 compares the eigenenergies as functions of the Rabi splitting, Ω 10 = 2g 10 √ N , of systems with different anharmonicities, and in which the frequency of the EM mode is resonant with the transition 0 → 1 of the emitters. As expected, the slope of the energy as a function of the coupling intensity is proportional to the number of quanta in the non-dark modes, with the sign of the contribution being positive for UP, and negative for LP. The introduction of anharmonicity lifts the degeneracy of states with multiple excitations in the dark modes, which has some remarkable consequences.
In the harmonic case, the quantum numbers facilitate the identification of the symmetry of a state just by inspection of its multiplicity. For instance, the state |3(UP) is unique in regards of its spectral configuration; therefore, it belongs to    In figure 3, the energy of the eigenstates is plotted against the number of coupled molecules with the thickness of the lines illustrating the degeneracy of each energy level. In each case, the thickness corresponds to log[dim(S λ )] and should not be mistaken by any broadening mechanism such as dissipation or disorder. It can be seen that both, degeneracy and energy separation, increase with the number of molecules. Furthermore, the most degenerate levels are those with constant energy. Among such states there are dark, with no component whatsoever from the EM mode, and states where v(UP) = v(LP) in which the contributions from both polaritonic modes balance each other energetically. This observation informs that light-matter coupling might not impact processes that engage the matter component, such as chemical reactions, even when anharmonicities are taken into account. Moreover, the resolution in energies also suggests that naturally occurring broadenings shall smear bands of nearby energy levels making them effectively indistinguishable.
The formalism presented in this work also allows to calculate observables associated with operators for which the permutational symmetry holds. One clear example is the photon number operator,v 0 =â † 0â 0 , which measures the photon contents of a given state. In the symmetrized basis, this operator is diagonal and has the form |n exc ,μ, λ; y i , w j n exc ,μ, λ; y i , w j |.
In figure 4 the photon content of the eigenstates is plotted as a function of the detuning for a fixed value of the collective coupling. The inclusion of anharmonicity should be inconsequential except for the change in the labeling of eigenstates.

Conclusions.
In this paper, we address the problem of the coupling between N identical dipoles, each with an arbitrary spectral structure, and a harmonic electromagnetic mode confined in a cavity. We have introduced tools from Group Theory that capitalize from the permutational symmetry of the system to simplify the time-independent Schrödinger equation. In the symmetry-adapted basis, the Hamiltonian breaks down into manageable matrices whose dimension is independent of the number of emitters and grows subexponentially with the number of excitations. While the total number of these matrices does depend polynomially on the number of emitters, they encode highly degenerate subspaces; therefore, the effective number of matrices to diagonalize is actually small. Since the method here presented does not rely on the explicit form of the basis transformation, the off-diagonal matrix elements are constructed by taking advantage of the fact that the structure of the Hilbert space also displays the symmetry of the special unitary group, and therefore can be described with the tools from angular momentum theory. This procedure exhibits the (super)radiant character of the transitions. We have also explored the immediate implications of including anharmonicities on the energetic and combinatorial characterization of the eigenstates. Finally, we have calculated the photon content of the eigenstates to exemplify the utility of the method in the extraction of observables associated with symmetry-preserving operators.    Table B1: Contribution of transitions in emitter space to the coupling coefficients betweenμ = 0 N −2 1 2 andμ ′ = 0 N −3 1 3 . Notice that nμ(0) = N − 2 and nμ′(1) = 3.
In this section, we present two approaches, one analytical and other computational, to the calculation of the contribution of the collective transitions in the Hilbert space of the emitters to the couplings between states in (54). First, we introduce a by-hand method to compute individual instances of (63). Taking advantage of the fact that nμ(v * + s) + nμ(v * ) = nμ′(v * + s) + nμ′(v * ), we can write L (λ) µ ′μ . When ρ J = 0, which is a typical situation among low excitation manifolds, the term inside the radical has an intuitive interpretation since it can be read as the product of the number of emitters inμ available for excitation with the number of emitters inμ ′ available for deexcitation. Furthermore, ρ J increases as λ becomes less dominant, i.e., strays away from [N]; this provides with an automatic way to determine the symmetry labels for which a given coupling is worth to calculate. table B1 illustrates this calculation for the couplings betweenμ = 0 N −2 1 2 andμ ′ = 0 N −3 1 3 .
Second, we discuss a computational method to get not only multiple values of the sought coefficients, but also the SABEs connected by a particular raising operator. Given thatĴ (μ ′ ,μ) + |μ, λ; y, w = L (λ) it is possible to writê µ ′ ,μ (w ′ , w)|μ ′ , λ; y i , w ′ j ′ μ, λ; y i , w j |, The latter implies that, if the raising operator is written in the basis of the PBEs, the coefficients L (λ) µ ′ ,μ (w ′ , w) are its singular values, while the unitary operatorsÛμ′ and Uμ correspond to the matrices of singular vectors, which provide with the changes of basis between PBEs and SABEs. In short, singular value decomposition ofĴ (μ ′ ,μ) + in an arbitrary basis yields the coupling coefficients as well as the symmetrized states in the subspaces of functions with the spectral configurations involved in the transition. The apparent contradiction in (B.3) that there are symmetry labels such that λ µ ′ yet λ µ is solved by recognizing the set of states {|μ ′ , λ; y, t ′ : λ µ} as the null-space of the raising operator, i.e., L (λ µ) µ ′ ,μ = 0. This feature is consistent with the fact that dim(M µ ′ ) > dim(M µ ).
To conclude, we work the example of the coupling, in a system with N = 4, betweeñ µ = 0 3 1 1 andμ ′ = 0 2 1 2 . The subspaces and couplings are given by In this basis, the raising operator is (B.6) Its singular value decomposition yields