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

Optical microcavities and similar devices enable the dipolar interaction between the electromagnetic (EM) field they confine and the degrees of freedom of quantum emitters in the same energy range [1–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 research fields such as quantum optics [6], excitonic and two-dimensional materials science [7–9], and chemistry [10, 11]. In many systems of interest, observable effects emerge only when a large number of atoms or molecules cooperatively interact with a single photon mode through their dipole moments [12, 13].

Polaritons have been produced using diverse setups and materials whose description requires various theoretical frameworks. For instance, real spins, nitrogen-vacancy centers [14], qubits [15], superconducting quantum interference devices [16], quantum dots [17, 18], and electronic transitions in organic molecules [19] are approached 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–27], studied from theory [28, 29], and motivated innovative theoretical ideas for signal enhancement [30–32].

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) [33–35]. From these studies, it is known that analytical solutions exist for the coupling of r-level systems to r − 1 EM modes [36]. In the context of the problem of N molecules coupled to a single EM mode, the RWA allows separating the Hamiltonian of the system enabling the computation of exact solutions. However, in this kind of system 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 [37, 38]. However, the literature on the stationary features of such systems, which are discussed in depth in the present work, is quite limited [39].

Permutational symmetry has been exploited in the few body limit to study systems such as nuclear structure [40], quantum circuits [41], magnons [42], ultracold atoms [43], nuclear spins within molecules [44], and the one-dimensional Hubbard model [45]. In the general case, the properties of the symmetrized quantum states have also been extensively discussed [39, 46]. This algebraic approach has been applied to the problem of light–matter coupling in the collective regime for two-level systems [46, 47]; however, to the best of our knowledge, there has not been an effort to present concrete results beyond this approximation.

There has been also a recently increasing number of efforts devoted to exactly solvable problems regarding light–matter coupling [48, 49], as well as their generalizations and extensions [50–53]. Among these, the many-body systems are particularly attractive [54–56].

In the present work, we provide a sophisticated yet insightful 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 symmetry, give rise to cooperative couplings that depart from the well-known factor of $\sqrt{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 block-diagonalization. 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.

The Hamiltonian

$\begin{equation}{\hat{H}}_{\text{bare}}(N)=\hslash \omega \left({\hat{a}}_{0}^{{\dagger}}{\hat{a}}_{0}+\frac{1}{2}\right)+\sum\limits _{i=1}^{N}\sum\limits _{v=0}^{r}{\varepsilon }_{v}{\hat{\sigma }}_{v,v}^{(i)},\end{equation} \tag{ 1 }$

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, ${\hat{a}}_{0}^{({\dagger})}$ is the annihilation (creation) operator of the EM mode, ɛ_v is the energy of the vth state of the emitter, and ${\hat{\sigma }}_{v,u}^{(i)}$ is the transition operator between levels u and v in the ith emitter, i.e.

$\begin{equation}{\hat{\sigma }}_{\alpha ,\beta }^{(i)}\left\vert {v}_{0}{v}_{1}\dots {v}_{i}\dots {v}_{N}\right\rangle ={\delta }_{{v}_{i}\beta }\left\vert {v}_{0}{v}_{1}\dots {\alpha }_{i}\dots {v}_{N}\right\rangle .\end{equation} \tag{ 2 }$

In (2), the eigenstates of ${\hat{H}}_{\text{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 is

$\begin{equation}\hat{H}(N)={\hat{H}}_{\text{bare}}(N)+{\hat{H}}_{\text{int}}(N),\end{equation} \tag{ 3 }$

where

$\begin{equation}{\hat{H}}_{\text{int}}(N)=\sum\limits _{v=0}^{r-1}\left({g}_{v,v+1}{\hat{a}}_{0}^{{\dagger}}{\hat{J}}_{-}^{(v+1,v)}+{g}_{v+1,v}{\hat{J}}_{+}^{(v+1,v)}{\hat{a}}_{0}\right),\end{equation} \tag{ 4 }$

with collective ladder operators for the emitters,

$\begin{equation}{\hat{J}}_{{\pm}}^{\left(v+\frac{1{\pm}1}{2}u,v-\frac{1\mp 1}{2}u\right)}=\sum\limits _{i=1}^{N}{\hat{\sigma }}_{v{\pm}u,v}^{(i)}\quad u{ >}0.\end{equation} \tag{ 5 }$

Here, the coupling constant ${g}_{v{\pm}u,v}=\sqrt{\frac{\hslash \omega }{2{{\epsilon}}_{0}\mathcal{V}}}\left\langle v{\pm}u\right\vert e\hat{x}\left\vert v\right\rangle$ is the product of the amplitude of the single-photon electric field, confined to a mode volume $\mathcal{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, which is a useful approximation at low energies.

The operator ${\hat{n}}_{\text{exc}}={\hat{a}}_{0}^{{\dagger}}{\hat{a}}_{0}+{\sum }_{i=1}^{N}{\sum }_{v=0}^{r}v{\hat{\sigma }}_{v,v}^{(i)}$ acts on the PBEs according to

$\begin{equation}{\hat{n}}_{\text{exc}}\left\vert {v}_{0}{v}_{1}\dots {v}_{N}\right\rangle =\sum\limits _{j=0}^{N}{v}_{j}\left\vert {v}_{0}{v}_{1}\dots {v}_{N}\right\rangle ,\end{equation} \tag{ 6 }$

thus indicating the total number of excitations of a given state; moreover, since $\left[\hat{H}(N),{\hat{n}}_{\text{exc}}\right]=0$, this operator corresponds to a constant of motion. Therefore, the Hamiltonian can be recast in the form

$\begin{equation}\hat{H}(N)=\sum\limits _{{n}_{\text{exc}}}{\hat{H}}_{{n}_{\text{exc}}}(N),\end{equation} \tag{ 7 }$

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. Appendix A displays their explicit form for 0 ⩽ n_exc ⩽ 3.

While separating (3) in excitation manifolds significantly simplifies the diagonalization of this Hamiltonian, it is important to remark that the matrices generated with this approach can quickly become intractable as the number of molecules increases. For instance, ${\hat{H}}_{{n}_{\text{exc}}}(N)$ is an operator in $\left(\genfrac{}{}{0.0pt}{}{N+{n}_{\text{exc}}}{{n}_{\text{exc}}}\right)$ 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 [57–59].

Further simplification of (3) requires to identify the additional symmetries of the Hamiltonian. An inspection of the explicit form of ${\hat{H}}_{{n}_{\text{exc}}}(N)$ (see appendix A) reveals that PBEs with the same distribution of quanta are degenerate eigenstates of ${\hat{H}}_{\text{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 operators

$\begin{equation}{\hat{J}}_{0}^{(u,v)}=\frac{1}{2}\left[{\hat{J}}_{+}^{(u,v)},{\hat{J}}_{-}^{(u,v)}\right]\quad 0{\leqslant}v{< }u{\leqslant}r,\end{equation} \tag{ 8 }$

which commute with ${\hat{H}}_{\text{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:

$\begin{equation}\tilde {\mu }=\prod\limits _{v=0}^{r}{v}^{{n}_{\tilde {\mu }}(v)},\end{equation} \tag{ 9 }$

a notation device that indicates the number of emitters ${n}_{\tilde {\mu }}(v)$ in the vth excited state. It must be noted that the product operator in (9) does not represent an arithmetic operation acting on integers, but rather a symbolic one acting on labels. See table B1 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; however, the spectral configuration itself does not include any information about the EM mode. The BEs can thus be partially characterized with three labels: the excitation manifold (n_exc), the spectral configuration ($\tilde {\mu }$), and a degeneracy index (y) whose meaning depends on the basis. These states fulfill

$\begin{equation}{\hat{H}}_{\text{bare}}(N)\left\vert {n}_{\text{exc}},\tilde {\mu };y\right\rangle ={\varepsilon }_{{n}_{\text{exc}}}^{(\tilde {\mu })}\left\vert {n}_{\text{exc}},\tilde {\mu };y\right\rangle ,\end{equation} \tag{ 10 }$

where

$\begin{equation}{\varepsilon }_{{n}_{\text{exc}}}^{(\tilde {\mu })}=\hslash \omega \left({v}_{0}^{({n}_{\text{exc}},\tilde {\mu })}+\frac{1}{2}\right)+\sum\limits _{v=0}^{{n}_{\text{exc}}}{n}_{\tilde {\mu }}(v){\varepsilon }_{v}\end{equation} \tag{ 11 }$

is the characteristic bare energy for states with those labels, and ${v}_{0}^{({n}_{\text{exc}},\tilde {\mu })}={n}_{\text{exc}}-{\sum }_{v=1}^{{n}_{\text{exc}}}v{n}_{\tilde {\mu }}(v)$ is the corresponding number of quanta in the EM mode.

Since ${\sum }_{v=0}^{r}{n}_{\tilde {\mu }}(v)=N$ for a given $\tilde {\mu }$, 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.

$\begin{equation}\boldsymbol{\mu }=[{\mu }_{1},{\mu }_{2},\dots ,{\mu }_{r+1}],\end{equation} \tag{ 12 }$

where ${\mu }_{i}={n}_{\tilde {\mu }}({v}_{i-1})$ under the constraint that μ_i ⩽ μ_i 'for i > i' [60]. In what follows, all the empty elements, μ_i = 0, will be omitted for brevity, as illustrated in table B1. Furthermore, the spectral configurations can be identified as partitions of the number of excitations distributed among the emitters, i.e. for every $\tilde {\mu }$ there exists a regular partition $\boldsymbol{\nu }(\tilde {\mu })$ such that

$\begin{equation}\boldsymbol{\nu }(\tilde {\mu })=\left\{{v}^{(i)}:1{\leqslant}v{\leqslant}r,0{\leqslant}i{\leqslant}{n}_{\tilde {\mu }}(v)\right\},\end{equation} \tag{ 13 }$

as illustrated in table B1. The above considerations lead to the result that the number of possible spectral configurations within a given manifold is

$\begin{equation}\vert {n}_{\text{exc}}\vert =\sum\limits _{n=0}^{{n}_{\text{exc}}}p(n),\end{equation} \tag{ 14 }$

where p(n) is the number of partitions of the integer n. These numbers can be extracted from the expansion [61]

$\begin{equation}\prod\limits _{k=1}^{m}{\left(1-{x}^{k}\right)}^{-1}=\sum\limits _{n=0}^{m}p(n){x}^{n}+O\left({x}^{m+1}\right).\end{equation} \tag{ 15 }$

Since the emitters are considered as identical, $\hat{H}(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 $\hat{\pi }\in {S}_{N}$, $\left[\hat{H}(N),\hat{\pi }\right]=0$.

The permutation operators act on the BEs according to

$\begin{equation}\hat{\pi }\left\vert {n}_{\text{exc}},\tilde {\mu };y\right\rangle =\left\vert {n}_{\text{exc}},\tilde {\mu };{y}^{\prime }\right\rangle ,\end{equation} \tag{ 16 }$

which means that BEs with common n_exc and $\tilde {\mu }$ carry a representation of S_N , as illustrated in appendix B.

A partition μ identifies representations of S_N that are, in general, reducible; they receive the name of permutation modules and are denoted by M ^μ [62]. The dimension of this representation is the number of BEs with the same μ , i.e.

$\begin{equation}\enspace \mathrm{dim}({M}^{\boldsymbol{\mu }})=\frac{N!}{\prod\limits _{i=1}^{r+1}{\mu }_{i}!}.\end{equation} \tag{ 17 }$

The permutation modules can, in turn, be decomposed in irreducible representations (irreps) according to Young's rule [63]:

$\begin{equation}{M}^{\boldsymbol{\mu }}=\underset{\boldsymbol{\lambda }{\trianglerighteq}\boldsymbol{\mu }}{\bigoplus}{K}_{\boldsymbol{\lambda }\boldsymbol{\mu }}{S}^{\boldsymbol{\lambda }},\end{equation} \tag{ 18 }$

where S ^λ symbolizes the irreps of the permutational group, also known as Specht modules [64].

The direct sum in (18) runs over all partitions λ of N that dominate (denoted by the symbol ⊵) over μ , i.e. they fulfill [63]

$\begin{equation}\sum\limits _{i=1}^{j}{\lambda }_{i}{\geqslant}\sum\limits _{i=1}^{j}{\mu }_{i}\quad \forall j:1{\leqslant}j{\leqslant}r.\end{equation} \tag{ 19 }$

For instance, consider the spectral configurations in the triply excited manifold: 0^N , 0^N−11¹, 0^N−21², 0^N−12¹, 0^N−31³, 0^N−21¹2¹, and 0^N−13¹. Since they correspond to the same regular partition, the spectral configura-tions 0^N−11¹, 0^N−12¹ and 0^N−13¹ have identical permutational properties. Let us identify how each partition relates to [N − 2, 1, 1] in terms of dominance. Following (19), it is possible to conclude that

$\begin{equation}\begin{aligned}& [N]{\trianglerighteq}[N-2,1,1],\\ & [N-1,1]{\trianglerighteq}[N-2,1,1],\\ & [N-2,2]{\trianglerighteq}[N-2,1,1],\\ & [N-3,3]{\ntrianglerighteq}[N-2,1,1],\\ & [N-2,1,1]{\trianglerighteq}[N-2,1,1].\end{aligned}\end{equation} \tag{ 20 }$

The coefficients K _λμ in (18) are known as Kostka numbers [65]; 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 [66, 67], these coefficients can be found through their combinatorial interpretation: the number of semi-standard Young tableaux (SSYT) of shape λ and content μ [63]. 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, i.e. with decreasing row-length from top to bottom, which is consistent with the regular form of the partitions (see table B1 for selected examples). An 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 [63]. For instance, the SSYT of shape [N − 1, 1] and content [N − 2, 1, 1] are

therefore, K_{[N−1,1],[N−2,1,1]} = 2.

With the above definitions, it can be verified that the representation of S_N in the basis of all BEs with $\tilde {\mu }={0}^{N-2}{1}^{1}{2}^{1}$ fulfills

$\begin{equation}{M}^{[N-2,1,1]}={S}^{[N]}\oplus 2{S}^{[N-1,1]}\oplus {S}^{[N-2,2]}\oplus {S}^{[N-2,1,1]}.\end{equation} \tag{ 22 }$

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 [63]. For instance, the SYTs for λ = [3, 2] are

This quantity is given by the hook-length formula [68]:

$\begin{equation}\enspace \mathrm{dim}({S}^{\boldsymbol{\lambda }})=\frac{N!}{{\prod }_{i=1}^{N}{h}_{\lambda }(i)},\end{equation} \tag{ 24 }$

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. An example of this computation is illustrated in appendix B.

A direct result from Young's rule is that

$\begin{equation}\enspace \mathrm{dim}({M}^{\boldsymbol{\mu }})=\sum\limits _{\boldsymbol{\lambda }{\trianglerighteq}\boldsymbol{\mu }}{K}_{\boldsymbol{\lambda }\boldsymbol{\mu }}\enspace \mathrm{dim}({S}^{\boldsymbol{\lambda }}).\end{equation} \tag{ 25 }$

The dimensions of the permutation and Specht modules for the spectral configurations found up to n_exc = 4 are displayed in table 1.

Table 1. Spectral configurations, EM mode excitations and degeneracies for the first five excitation manifolds.

	nexc
0	1	2	3	4
$\tilde {\mu }$	${v}_{0}^{({n}_{\text{exc}},\tilde {\mu })}$					dim(M μ )	$\mathrm{dim}({S}^{\boldsymbol{\lambda }})$ a
0N	0	1	2	3	4	1	1
0N−111		0	1	2	3	N	N − 1
0N−212			0	1	2	N(N − 1)/2	N(N − 3)/2
0N−121			0	1	2	N
0N−313				0	1	N(N − 1)(N − 2)/6	N(N − 1)(N − 5)/6
0N−21121				0	1	N(N − 1)	(N − 1)(N − 2)/2
0N−131				0	1	N
0N−414					0	N(N − 1)(N − 2)(N − 3)/24	N(N − 1)(N − 2)(N − 7)/24
0N−31221					0	N(N − 1)(N − 2)/2	N(N − 2)(N − 4)/3
0N−222					0	N(N − 1)/2
0N−21131					0	N(N − 1)
0N−141					0	N
|nexc|	1	2	4	7	12
$\vert {\lambda }_{{n}_{\text{exc}}}\vert$	1	2	3	5	7

^a λ = μ when applicable.

The irreps identify the smallest possible subspaces that remain mutually excluded under permutations. In other words, the space of BEs with the same $\tilde {\mu }$ can be split into sets of symmetry adapted linear combinations of BEs (SABEs) labeled by λ for which

$\begin{equation}\hat{\pi }\left\vert {n}_{\text{exc}},\tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\rangle =\left\vert {n}_{\text{exc}},\tilde {\mu },\boldsymbol{\lambda };{\mathbf{y}}^{\prime },\mathbf{w}\right\rangle ,\end{equation} \tag{ 26 }$

for all $\hat{\pi }\in {S}_{n}$, where y identifies a particular SYT of shape λ , and 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

$\begin{equation}{\left({\mathbb{C}}^{r+1}\right)}^{\otimes N}\cong \underset{\boldsymbol{\lambda }{\vdash}N}{\bigoplus}\left({S}^{\boldsymbol{\lambda }}\otimes {V}_{r+1}^{\boldsymbol{\lambda }}\right),\end{equation} \tag{ 27 }$

where ${\mathbb{C}}^{r+1}$ is the vector space spanned by the energy eigenstates of each molecule, and therefore ${\left({\mathbb{C}}^{r+1}\right)}^{\otimes N}$ is the vector space spanned by the BEs. The expression λ ⊢ N indicates that λ is a partition of N, and ${V}_{r+1}^{\boldsymbol{\lambda }}$ 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 [69, 70], 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 acting on the internal degrees of freedom of each molecule. Furthermore, the decomposition in (27) provides with the meaning of all the indices in (26).

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 an 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 B2.

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}^{\boldsymbol{\lambda }}$ can be enumerated with SSYT, and thus uniquely identified with the indices $\tilde {\mu }$, λ , 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 B3.

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 $\tilde {\mu }={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^[N] = S^[N]. These facts are consistent with the uniqueness of the state within each manifold where all the excitations reside in the EM mode, which is denoted by $\left\vert {n}_{\text{exc}},{0}^{N},[N];0,0\right\rangle$.

Since $\left[\hat{\pi },{\hat{J}}_{+}^{(u,v)}\right]=0$, the states $\left\vert {n}_{\text{exc}},\tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\rangle$ behave identically to ${\hat{a}}_{0}{\hat{J}}_{+}^{u,v}\left\vert {n}_{\text{exc}},\tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\rangle$ under permutations. This means that, even though the operator ${\hat{a}}_{0}{\hat{J}}_{+}^{(u,v)}$ modifies the spectral configuration, the states it couples carry the same irrep. Consequently

$\begin{equation}{\hat{a}}_{0}{\hat{J}}_{+}^{(u,0)}\left\vert {n}_{\text{exc}},{0}^{N},[N];0,0\right\rangle \propto \left\vert {n}_{\text{exc}},{0}^{N-1}{v}^{1},[N];0,0\right\rangle .\end{equation} \tag{ 28 }$

However, M^[N−1,1] = S^[N] ⊕ S^[N−1,1], which means that the remaining SABEs with $\tilde {\mu }={0}^{N-1}{u}^{1}$ carry the irrep with λ = [N − 1, 1]. These states, $\left\vert {n}_{\text{exc}},{0}^{N-1}{v}^{1},[N-1,1];k,{1}^{1}\right\rangle$ with 1 ⩽ k ⩽ N − 1, can be obtained through Gram–Schmidt orthogonalization over the basis of PBEs but with $\left\vert {n}_{\text{exc}},{0}^{N-1}{v}^{1},[N];0,0\right\rangle$ replacing one of the states and remaining fixed as seed of the procedure. The structure of the doubly excited manifold can be understood in terms of this procedure as illustrated in figure 1. Appendix B includes the example of reduced representations of the group S₃, and appendix C discusses the generalities of the expression of the SABEs in terms of the PBEs.

Zoom In Zoom Out Reset image size

The spectral configurations, $\tilde {\mu }$, giving rise to partitions that are valid labels, λ , for the irreps fulfill

$\begin{equation}{n}_{\tilde {\mu }}(v){\geqslant}{n}_{\tilde {\mu }}({v}^{\prime })\quad \forall \enspace v,{v}^{\prime }:{v}^{\prime }{ >}v{ >}0;\end{equation} \tag{ 29 }$

this constraint produces the so-called partitions with weakly decreasing multiplicities [71]. The enumeration of these partitions imply that the total number of irreps in the n_excth manifold is

$\begin{equation}\vert {\lambda }_{{n}_{\text{exc}}}\vert =\sum\limits _{n=0}^{{n}_{\text{exc}}}q(n),\end{equation} \tag{ 30 }$

where q(n) is the number of partitions of n with weakly decreasing multiplicities as illustrated in table 1. These numbers can be extracted from the expansion [60, 61]

$\begin{equation}\prod\limits _{k=2}^{m}{\left(1-{x}^{\left(\genfrac{}{}{0.0pt}{}{k}{2}\right)}\right)}^{-1}=\sum\limits _{n=0}^{m}q(n){x}^{n}+O\left({x}^{m+1}\right).\end{equation} \tag{ 31 }$

The main result of this section is that the Hamiltonians of the excitation manifolds can be split according to

$\begin{equation}{\hat{H}}_{{n}_{\text{exc}}}(N)=\sum\limits _{i=1}^{\left\vert {\lambda }_{{n}_{\text{exc}}}\right\vert }{\hat{H}}_{{n}_{\text{exc}}}^{({\boldsymbol{\lambda }}_{i})},\end{equation} \tag{ 32 }$

where

$\begin{equation}{\hat{H}}_{{n}_{\text{exc}}}^{(\boldsymbol{\lambda })}=\sum\limits _{j=1}^{\mathrm{dim}\left({S}^{\boldsymbol{\lambda }}\right)}{\hat{H}}_{{n}_{\text{exc}}}^{\left(\boldsymbol{\lambda },{\mathbf{y}}_{j}\right)},\end{equation} \tag{ 33 }$

with each ${\hat{H}}_{{n}_{\text{exc}}}^{(\boldsymbol{\lambda },\mathbf{y})}$ encoding the energies and interactions of spectral configurations with the same symmetry. Equation (33) implies that the states obtained from diagonalizing ${\hat{H}}_{{n}_{\text{exc}}}^{(\boldsymbol{\lambda })}$ are dim(S ^λ )-fold degenerate. In particular, for the diagonal part of the Hamiltonian, we have

$\begin{equation}{\hat{H}}_{\text{bare},{n}_{\text{exc}}}^{(\boldsymbol{\lambda })}=\sum\limits _{\tilde {\mu }:\boldsymbol{\mu }{\trianglelefteq}\boldsymbol{\lambda }}\sum\limits _{i=1}^{\mathrm{dim}({S}^{\boldsymbol{\lambda }})}\sum\limits _{j=1}^{{K}_{\boldsymbol{\lambda }\boldsymbol{\mu }}}{\varepsilon }_{{n}_{\text{exc}}}^{(\tilde {\mu })}\left\vert {n}_{\text{exc}},\tilde {\mu },\boldsymbol{\lambda };{\mathbf{y}}_{i},{\mathbf{w}}_{j}\right\rangle \left\langle {n}_{\text{exc}},\tilde {\mu },\boldsymbol{\lambda };{\mathbf{y}}_{i},{\mathbf{w}}_{j}\right\vert ,\end{equation} \tag{ 34 }$

Notice that

$\begin{equation}\enspace \mathrm{dim}\left({\hat{H}}_{{n}_{\text{exc}}}^{(\boldsymbol{\lambda },\mathbf{y})}\right)=\sum\limits _{\boldsymbol{\mu }{\trianglelefteq}\boldsymbol{\lambda }}{K}_{\boldsymbol{\lambda }\boldsymbol{\mu }}.\end{equation} \tag{ 35 }$

The simplification achieved with this strategy is significant since, for each excitation manifold, the problem has been reduced to the diagonalization of $\vert {\lambda }_{{n}_{\text{exc}}}\vert \sim \enspace \mathrm{exp}(A{n}_{\text{exc}}^{1/3})$ matrices with dimensions in the neighborhood of $\vert {n}_{\text{exc}}\vert \sim \mathrm{exp}(\sqrt{C{n}_{\text{exc}}})$ [72], both independent of N, which is a substantial gain over the original $\left(\genfrac{}{}{0.0pt}{}{N+{n}_{\text{exc}}}{{n}_{\text{exc}}}\right)$-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.

After determining the SABEs, it is possible to evaluate the matrix elements of the interaction Hamiltonian. Let us consider the collective transition of a quantum from level v^* to v^* +  s, i.e. the action of the operator ${\hat{J}}_{+}^{{v}^{{\ast}}+s,{v}^{{\ast}}}$ on a SABE with spectral configuration $\tilde {\mu }$. As a result of the transition, the new SABE belongs to the spectral configuration ${\tilde {\mu }}^{\prime }$, which relates to the initial spectral configuration according to

$\begin{equation}s=\sum\limits _{v=1}^{{n}_{\text{exc}}}v[{n}_{{\tilde {\mu }}^{\prime }}(v)-{n}_{\tilde {\mu }}(v)].\end{equation} \tag{ 36 }$

Additionally, the occupation of levels in the two spectral configurations fulfill

$\begin{equation}{n}_{\tilde {\mu }}({v}^{{\ast}})-{n}_{{\tilde {\mu }}^{\prime }}({v}^{{\ast}})={n}_{{\tilde {\mu }}^{\prime }}({v}^{{\ast}}+s)-{n}_{\tilde {\mu }}({v}^{{\ast}}+s)=1.\end{equation} \tag{ 37 }$

Given the definition of ${\hat{H}}_{\text{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 $\left\{{\tilde {\mu }}^{\prime },\tilde {\mu }\right\}$ afforded by photonic (de)excitations is

$\begin{equation}{n}_{\text{int}}({n}_{\text{exc}})=\sum\limits _{n=0}^{{n}_{\text{exc}}-1}\sum\limits _{m=0}^{n}p(m).\end{equation} \tag{ 38 }$

The matrix elements of the interaction Hamiltonian are

$\begin{equation}\left\langle {n}_{\text{exc}}^{\prime },{\tilde {\mu }}^{\prime },{\boldsymbol{\lambda }}^{\prime };{\mathbf{y}}^{\prime },{\mathbf{w}}^{\prime }\right\vert {\hat{H}}_{\text{int}}(N)\left\vert {n}_{\text{exc}},\tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\rangle ={g}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}{B}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{({n}_{\text{exc}})}{L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}{C}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}({\mathbf{w}}^{\prime },\mathbf{w}){\delta }_{{n}_{\text{exc}}^{\prime }{n}_{\text{exc}}}{\delta }_{{\boldsymbol{\lambda }}^{\prime }\boldsymbol{\lambda }}{\delta }_{{\mathbf{y}}^{\prime }\mathbf{y}},\end{equation} \tag{ 39 }$

where, in agreement with the discussion above,

$\begin{equation}{g}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}={g}_{{v}^{{\ast}}+s,{v}^{{\ast}}}{\delta }_{\vert s\vert ,1}{\delta }_{{n}_{{\tilde {\mu }}^{\prime }}({v}^{{\ast}}),{n}_{\tilde {\mu }}({v}^{{\ast}})-1}{\delta }_{{n}_{{\tilde {\mu }}^{\prime }}({v}^{{\ast}}+s),{n}_{\tilde {\mu }}({v}^{{\ast}}+s)+1},\end{equation} \tag{ 40 }$

while the coefficient

$\begin{equation}{B}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{({n}_{\text{exc}})}={\left({v}_{0}^{({n}_{\text{exc}},\tilde {\mu })}+{\delta }_{s,-1}\right)}^{1/2}\end{equation} \tag{ 41 }$

accounts for the (de)excitation of the EM mode.

To describe the redistribution of quanta among the emitters, let us define the operators

$\begin{equation}{\hat{J}}_{r}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}={\hat{J}}_{r}^{({v}^{{\ast}}+s,{v}^{{\ast}})},\end{equation} \tag{ 42 }$

where r ∈ {−, 0, +}, that generate the $\mathfrak{s}\mathfrak{u}(2)$ algebra with Casimir operator

$\begin{equation}{\hat{\mathbf{J}}}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{2}={\left({\hat{J}}_{0}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}\right)}^{2}+\frac{1}{2}\left({\hat{J}}_{+}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}{\hat{J}}_{-}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}+{\hat{J}}_{-}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}{\hat{J}}_{+}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}\right).\end{equation} \tag{ 43 }$

These operators allow the definition of the quantities

$\begin{equation}{M}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}=\left\langle \tilde {\mu }\right\vert {\hat{J}}_{0}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}\left\vert \tilde {\mu }\right\rangle \end{equation} \tag{ 44 }$

and

$\begin{equation}{J}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}\left({J}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}+1\right)=\left\langle \tilde {\mu },\boldsymbol{\lambda }\right\vert {\hat{\mathbf{J}}}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{2}\left\vert \tilde {\mu },\boldsymbol{\lambda }\right\rangle ,\end{equation} \tag{ 45 }$

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}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}$ as the Dicke cooperation number [73]. To be specific

$\begin{equation}{M}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}=\frac{{n}_{\tilde {\mu }}({v}^{{\ast}}+s)-{n}_{\tilde {\mu }}({v}^{{\ast}})}{2},\end{equation} \tag{ 46 }$

and

$\begin{equation}{J}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}=\underset{{\tilde {\mu }}^{{\prime\prime}}}{\mathrm{max}}\left(\left\langle {\tilde {\mu }}^{{\prime\prime}},\boldsymbol{\lambda }\right\vert {\hat{J}}_{0}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}\left\vert {\tilde {\mu }}^{{\prime\prime}},\boldsymbol{\lambda }\right\rangle \right).\end{equation} \tag{ 47 }$

According to the prescription of quantum angular momentum, the contribution of collective (de)excitations to the coupling between states is given by

$\begin{equation}{L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}={\left[\left({J}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}+{M}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}\right)\left({J}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}-{M}_{\tilde {\mu },{\tilde {\mu }}^{\prime }}\right)\right]}^{1/2}.\end{equation} \tag{ 48 }$

Couplings are said to be superradiant if ${L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}{ >}1$, and subradiant if ${L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}{< }1$ [73]. Appendix D shows a couple of neat ways to calculate these coefficients.

The last component of (39), the coefficient ${C}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}({\mathbf{w}}^{\prime },\mathbf{w})$, acknowledges the fact that, since $\left[{\hat{J}}_{+}^{({v}^{\prime },v)},{\hat{J}}_{+}^{({u}^{\prime },u)}\right]\propto 1-{\delta }_{uv}{\delta }_{{u}^{\prime }{v}^{\prime }}$, 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 ${\hat{H}}_{{n}_{\text{exc}}}^{(\boldsymbol{\lambda },\mathbf{y})}$ requires more than one function with the same $\tilde {\mu }$, i,e, when K _λμ > 1; otherwise, ${C}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}({\mathbf{w}}^{\prime },\mathbf{w})=1$. For instance, consider the state $\left\vert 3,{0}^{N-1}{1}^{1},[N-1,1];2,1\right\rangle =(\left\vert {2}_{0}{1}_{1}\right\rangle -\left\vert {2}_{0}{1}_{2}\right\rangle )/\sqrt{2}$, and apply to it consecutive excitations to get to a state with $\tilde {\mu }={0}^{N-2}{1}^{1}{2}^{1}$. There are two possible paths:

$\begin{align}\frac{\left\vert {2}_{0}{1}_{1}\right\rangle -\left\vert {2}_{0}{1}_{2}\right\rangle }{\sqrt{2}}& {\hat{a}}_{0}{\hat{J}}_{+}^{(1,0)}{\to }\sum\limits _{i=3}^{N}\frac{\left\vert {1}_{0}{1}_{1}{1}_{i}\right\rangle -\left\vert {1}_{0}{1}_{2}{1}_{i}\right\rangle }{\sqrt{2(N-2)}}\\ & {\hat{a}}_{0}{\hat{J}}_{+}^{(2,1)}{\to }\sum\limits _{i=3}^{N}\frac{\left\vert {2}_{1}{1}_{i}\right\rangle -\left\vert {2}_{2}{1}_{i}\right\rangle +\left\vert {1}_{1}{2}_{i}\right\rangle -\left\vert {1}_{2}{2}_{i}\right\rangle }{2\sqrt{N-2}},\end{align} \tag{ 49a }$

and

$\begin{align}\frac{\left\vert {2}_{0}{1}_{1}\right\rangle -\left\vert {2}_{0}{1}_{2}\right\rangle }{\sqrt{2}}& {\hat{a}}_{0}{\hat{J}}_{+}^{(2,1)}{\to }\frac{\left\vert {1}_{0}{2}_{1}\right\rangle -\left\vert {1}_{0}{2}_{2}\right\rangle }{\sqrt{2}}\\ & {\hat{a}}_{0}{\hat{J}}_{+}^{(1,0)}{\to }\frac{1}{\sqrt{2(N-1)}}\left[\left\vert {2}_{1}{1}_{2}\right\rangle -\left\vert {1}_{1}{2}_{2}\right\rangle +\sum\limits _{i=3}^{N}(\left\vert {2}_{1}{1}_{i}\right\rangle -\left\vert {2}_{2}{1}_{i}\right\rangle )\right].\end{align} \tag{ 49b }$

While the identifications

$\begin{equation}\begin{aligned}\left\vert 3,{0}^{N-2}{1}^{2},[N-1,1];2,1\right\rangle & =\sum\limits _{i=3}^{N}\frac{\left\vert {1}_{0}{1}_{1}{1}_{i}\right\rangle -\left\vert {1}_{0}{1}_{2}{1}_{i}\right\rangle }{\sqrt{2(N-2)}},\\ \left\vert 3,{0}^{N-1}{2}^{1},[N-1,1];2,1\right\rangle & =\frac{\left\vert {1}_{0}{2}_{1}\right\rangle -\left\vert {1}_{0}{2}_{2}\right\rangle }{\sqrt{2}}\end{aligned}\end{equation} \tag{ 50 }$

are immediate, the assignment of $\left\vert 3,{0}^{N-2}{1}^{1}{2}^{1},[N-1,1];2,1\right\rangle$ and $\left\vert 3,{0}^{N-2}{1}^{1}{2}^{1},[N-1,1];2,2\right\rangle$ is unclear since the two obtained states are linearly independent yet not mutually orthogonal. Bypassing this ambiguity requires the definition of an orthonormal basis spanned by the previous states, then the coefficients can be extracted from the identity

$\begin{equation}{L}_{{\tilde {\mu }}^{\prime }\tilde {\mu }}^{(\boldsymbol{\lambda })}{C}_{{\tilde {\mu }}^{\prime }\tilde {\mu }}^{(\boldsymbol{\lambda })}({\mathbf{w}}^{\prime },\mathbf{w})=\left\langle \tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\vert {\hat{J}}_{-}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}\left\vert {\tilde {\mu }}^{\prime },\boldsymbol{\lambda };\mathbf{y},{\mathbf{w}}^{\prime }\right\rangle .\end{equation} \tag{ 51 }$

For instance, let us set

$\begin{equation}\left\vert A\right\rangle =\sum\limits _{i=3}^{N}\frac{\left\vert {2}_{1}{1}_{i}\right\rangle -\left\vert {2}_{2}{1}_{i}\right\rangle +\left\vert {1}_{1}{2}_{i}\right\rangle -\left\vert {1}_{2}{2}_{i}\right\rangle }{2\sqrt{N-2}},\end{equation} \tag{ 52 }$

and

$\begin{equation}\left\vert X\right\rangle =\frac{\left\vert {2}_{1}{1}_{2}\right\rangle -\left\vert {1}_{1}{2}_{2}\right\rangle +{\sum }_{i=3}^{N}(\left\vert {2}_{1}{1}_{i}\right\rangle -\left\vert {2}_{2}{1}_{i}\right\rangle )}{\sqrt{2(N-1)}},\end{equation} \tag{ 53 }$

and compute their orthogonal complements:

$\begin{equation}\left\vert B\right\rangle =\frac{\hat{I}-\left\vert A\right\rangle \left\langle A\right\vert }{{\left(1-\vert \left\langle X\vert A\right\rangle {\vert }^{2}\right)}^{1/2}}\left\vert X\right\rangle ,\end{equation} \tag{ 54 }$

and

$\begin{equation}\left\vert Y\right\rangle =\frac{\hat{I}-\left\vert X\right\rangle \left\langle X\right\vert }{{\left(1-\vert \left\langle X\vert A\right\rangle {\vert }^{2}\right)}^{1/2}}\left\vert A\right\rangle ,\end{equation} \tag{ 55 }$

where $\hat{I}$ is the identity operator. The values of the coefficients calculated with both working basis are shown in table 2. Since the bare energies, ${\varepsilon }_{{n}_{\text{exc}}}^{(\tilde {\mu })}$, are independent of the index w, and ${\sum }_{{t}^{\prime }=1}^{{K}_{{\boldsymbol{\lambda }\boldsymbol{\mu }}^{\prime }}}{\left\vert {C}_{{\tilde {\mu }}^{\prime }\tilde {\mu }}^{(\boldsymbol{\lambda })}({\mathbf{w}}^{\prime },\mathbf{w})\right\vert }^{2}=1$, the choice of basis is immaterial for the eigenvalues of ${\hat{H}}_{{n}_{\text{exc}}}^{(\boldsymbol{\lambda },\mathbf{y})}$.

Table 2. Coefficients of projected couplings involving ${\tilde {\mu }}^{\prime }={0}^{N-2}{1}^{1}{2}^{1}$ for repeated irreps with λ = [N − 1, 1].

$\left\vert 3,{\tilde {\mu }}^{\prime },\boldsymbol{\lambda };2,1\right\rangle$	$\left\vert 3,{\tilde {\mu }}^{\prime },\boldsymbol{\lambda };2,2\right\rangle$	$\tilde {\mu }$	${C}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}({1}^{1},{1}^{1})$	${C}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}({1}^{1},{2}^{1})$
$\left\vert A\right\rangle$	$\left\vert B\right\rangle$	0N−212	1	0
		0N−121	$\left\langle X\vert A\right\rangle$	$\left\langle X\vert B\right\rangle$
$\left\vert X\right\rangle$	$\left\vert Y\right\rangle$	0N−212	$\left\langle A\vert X\right\rangle$	$\left\langle A\vert Y\right\rangle$
		0N−121	1	0

In summary, the full symmetrized blocks of the interaction Hamiltonian are given by

$\begin{equation}{\hat{H}}_{\text{int},{n}_{\text{exc}}}^{(\boldsymbol{\lambda },\mathbf{y})}=\sum\limits _{\tilde {\mu }:\boldsymbol{\mu }{\trianglelefteq}\boldsymbol{\lambda }}\;\sum\limits _{{\tilde {\mu }}^{\prime }:{\boldsymbol{\mu }}^{\prime }{\trianglelefteq}\boldsymbol{\lambda }}\sum\limits _{i=1}^{{K}_{\boldsymbol{\lambda }\boldsymbol{\mu }}}\sum\limits _{{i}^{\prime }=1}^{{K}_{\boldsymbol{\lambda }{\boldsymbol{\mu }}^{\prime }}}{g}_{{\tilde {\mu }}^{\prime }\tilde {\mu }}{B}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{({n}_{\text{exc}})}{L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}{C}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}({\mathbf{w}}_{{i}^{\prime }}^{\prime },{\mathbf{w}}_{i})\left\vert {n}_{\text{exc}},{\tilde {\mu }}^{\prime },\boldsymbol{\lambda };\mathbf{y},{\mathbf{w}}_{{i}^{\prime }}^{\prime }\right\rangle \left\langle {n}_{\text{exc}},\tilde {\mu },\boldsymbol{\lambda };\mathbf{y},{\mathbf{w}}_{i}\right\vert .\end{equation} \tag{ 56 }$

In this section we lay out the specific steps to write the block diagonalized Hamiltonian of a given excitation manifold.

• (a)  
  Define n_exc, the manifold.
• (b)  
  List all relevant $\tilde {\mu }$, the spectral configurations.

• (c)  
  Calculate all the relevant bare energies, ${\varepsilon }_{{n}_{\text{exc}}}^{(\tilde {\mu })}$, with (11).
• (d)  
  List all the allowed λ , the symmetry defining partitions.

  • 1.  
    Enumerate the spectral configurations as μ , i.e. as partitions in regular form.
  • 2.  
    Discard redundancies. The remaining partitions are the possible λ . Their number must add up to $\vert {\lambda }_{{n}_{\text{exc}}}\vert$ as calculated with (30) and (31).
• (e)  
  Determine the composition of each $\tilde {\mu }$ in terms of λ according to Young's rule in (18).

  • 1.  
    Find out the dominance relations among partitions according to (19).
  • 2.  
    Compute the pertinent Kostka numbers (K _λμ ).
• (f)  
  Write the diagonal of ${\hat{H}}_{{n}_{\text{exc}}}^{(\boldsymbol{\lambda })}$.

• (g)  
  Calculate the couplings.

• (h)  
  Collect the calculated couplings in their respective blocks.

6.1. 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 3 compiles the quantities computed according to steps (a)–(d) of the algorithm in the previous section, as well as the dimensions of the irreps in step (f) 2.

Table 3. Relevant partitions, spectral configurations, excitations in the EM mode, bare energies, and module dimensions calculated in the block diagonalization of ${\hat{H}}_{3}(N)$.

n	p(n)	ν	$\tilde {\mu }$	${v}_{0}^{(3,\tilde {\mu })}$	${\varepsilon }_{3}^{(\tilde {\mu })}$	μ	q(n)	λ	dim(S λ )
0	1	[0]	0N	3	$\frac{7}{2}\hslash \omega +N{\varepsilon }_{0}$	[N]	1	[N]	1
1	1	[1]	0N−111	2	$\frac{5}{2}\hslash \omega +(N-1){\varepsilon }_{0}+{\varepsilon }_{1}$	[N − 1, 1]	1	[N − 1, 1]	N − 1
2	2	[2]	0N−121	1	$\frac{3}{2}\hslash \omega +(N-1){\varepsilon }_{0}+{\varepsilon }_{2}$	[N − 1, 1]	1	—	—
		[1, 1]	0N−212	1	$\frac{3}{2}\hslash \omega +(N-2){\varepsilon }_{0}+2{\varepsilon }_{1}$	[N − 2, 2]		[N − 2, 2]	N(N − 3)/2
3	3	[3]	0N−131	0	$\frac{\hslash \omega }{2}+(N-1){\varepsilon }_{0}+{\varepsilon }_{3}$	[N − 1, 1]	2	—	—
		[2, 1]	0N−21121	0	$\frac{\hslash \omega }{2}+(N-2){\varepsilon }_{0}+{\varepsilon }_{1}+{\varepsilon }_{2}$	[N − 2, 1, 1]		[N − 2, 1, 1]	(N − 1)(N − 2)/2
		[1, 1, 1]	0N−313	0	$\frac{\hslash \omega }{2}+(N-3){\varepsilon }_{0}+3{\varepsilon }_{1}$	[N − 3, 3]		[N − 3, 3]	N(N − 1)(N − 5)/6
	$\left\vert {\tilde {\mu }}_{3}\right\vert =7$						$\left\vert {\lambda }_{3}\right\vert =5$

In table 4, we show the Kostka numbers that indicate the composition of the spaces spanned by wavefunctions with the same $\tilde {\mu }$ in terms of symmetrized subspaces labeled by λ . Notice that K _λμ = 0 implies that $\boldsymbol{\lambda }{\ntrianglerighteq}\boldsymbol{\mu }$.

Table 4. Kostka numbers, K _λμ , relating the permutation and Specht modules that appear in the triply excited manifold.

	μ
λ	[N]	[N − 1, 1]	[N − 2, 2]	[N − 2, 1, 1]	[N − 3, 3]
[N]	1	1	1	1	1
[N − 1, 1]	0	1	1	2	1
[N − 2, 2]	0	0	1	1	1
[N − 2, 1, 1]	0	0	0	1	0
[N − 3, 3]	0	0	0	0	1

The factors required to calculate the off-diagonal elements, as prescribed by step (g), up until [(g) 3. (i)], are on display in table 5. Step (g) 3. (ii) was worked in detail in the discussion leading to table 2. What is left to this section are the explicit forms of the overlaps:

$\begin{equation*}\left\langle X\vert A\right\rangle ={\left(\frac{N-2}{2(N-1)}\right)}^{1/2},\end{equation*}$

and

$\begin{equation*}\left\langle X\vert B\right\rangle ={\left(\frac{N}{2(N-1)}\right)}^{1/2}=\left\langle A\vert Y\right\rangle .\end{equation*}$

Table 5. Energy levels involved and contributions from transitions in the EM and emitter modes to the couplings between spectral configurations in the triply excited manifold.

$\tilde {\mu }$	${\tilde {\mu }}^{\prime }$	v*	v* + s	${B}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(3)}$	${L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{[N]}$	${L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{[N-1,1]}$	${L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{[N-2,2]}$
0N	0N−111	0	1	$\sqrt{3}$	$\sqrt{N}$
0N−111	0N−212	0	1	$\sqrt{2}$	$\sqrt{2(N-1)}$	$\sqrt{N-2}$
0N−111	0N−121	1	2	$\sqrt{2}$	1	1
0N−212	0N−313	0	1	1	$\sqrt{3(N-2)}$	$\sqrt{2(N-3)}$	$\sqrt{N-4}$
0N−212	0N−21121	1	2	1	$\sqrt{2}$	$\sqrt{2}$	$\sqrt{2}$
0N−121	0N−21121	0	1	1	$\sqrt{N-1}$	$\sqrt{N-1}$
0N−121	0N−131	2	3	1	1	1

With all these elements, the Hopfield–Bogoliubov forms of the block-diagonalized Hamiltonian are

$\begin{equation}{\mathbf{H}}_{3}^{[N-3,3]}={\mathbf{1}}_{N(N-1)(N-5)/6}\otimes \left(\begin{matrix}\hfill {\varepsilon }_{3}^{({0}^{N-3}{1}^{3})}\hfill \end{matrix}\right),\end{equation} \tag{ 60d }$

and

$\begin{equation}{\mathbf{H}}_{3}^{[N-2,1,1]}={\mathbf{1}}_{(N-1)(N-2)/2}\otimes \left(\begin{matrix}\hfill {\varepsilon }_{3}^{({0}^{N-2}{1}^{1}{2}^{1})}\hfill \end{matrix}\right),\end{equation} \tag{ 60e }$

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 (60b), the couplings involving the states with $\tilde {\mu }={0}^{N-2}{1}^{1}{2}^{1}$ were worked out with the basis $\left\{\left\vert A\right\rangle ,\left\vert B\right\rangle \right\}$.

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 $\tilde {\mu }={0}^{N-2}{1}^{1}{2}^{1}$ to this symmetry. As in the case above, the remaining SABE interacts only weakly with this array. For the symmetry λ = [N − 2, 2] the subspace corresponds to a three-level system with configuration Λ and degeneracy degree of N(N − 3)/2. However, only one of the couplings is superradiant while the other is weak. The remaining subspaces with symmetries [N − 3, 3] and [N − 2, 1, 1] are dark, and therefore there are no interactions among the states with degeneracy degrees N(N − 1)(N − 5)/6 and (N − 1)(N − 2)/2, respectively.

6.2. 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 [12, 74, 75].

$\begin{equation}{\mathbf{H}}_{0}^{[N]}=\left(\begin{matrix}\hfill {\varepsilon }_{0}^{({0}^{N})}\hfill \end{matrix}\right),\end{equation} \tag{ 61 }$

$\begin{equation}{\mathbf{H}}_{1}^{[N]}=\left(\begin{matrix}\hfill {\varepsilon }_{1}^{({0}^{N})}\hfill & \hfill {g}_{10}\sqrt{N}\hfill \\ \hfill {g}_{01}\sqrt{N}\hfill & \hfill {\varepsilon }_{1}^{({0}^{N-1}{1}^{1})}\hfill \end{matrix}\right),\end{equation} \tag{ 62a }$

$\begin{equation}{\mathbf{H}}_{1}^{[N-1,1]}={\mathbf{1}}_{N-1}\otimes \left(\begin{matrix}\hfill {\varepsilon }_{1}^{({0}^{N-1}{1}^{1})}\hfill \end{matrix}\right),\end{equation} \tag{ 62b }$

$\begin{equation}{\mathbf{H}}_{2}^{[N]}=\left(\begin{matrix}\hfill {\varepsilon }_{2}^{({0}^{N})}\hfill & \hfill {g}_{10}\sqrt{2N}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {g}_{01}\sqrt{2N}\hfill & \hfill {\varepsilon }_{2}^{({0}^{N-1}{1}^{1})}\hfill & \hfill {g}_{10}\sqrt{2(N-1)}\hfill & \hfill {g}_{21}\hfill \\ \hfill 0\hfill & \hfill {g}_{01}\sqrt{2(N-1)}\hfill & \hfill {\varepsilon }_{2}^{({0}^{N-2}{1}^{2})}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {g}_{12}\hfill & \hfill 0\hfill & \hfill {\varepsilon }_{2}^{({0}^{N-1}{2}^{1})}\hfill \end{matrix}\right),\end{equation} \tag{ 63a }$

$\begin{equation}{\mathbf{H}}_{2}^{[N-1,1]}={\mathbf{1}}_{N-1}\otimes \left(\begin{matrix}\hfill {\varepsilon }_{2}^{({0}^{N-1}{1}^{1})}\hfill & \hfill {g}_{10}\sqrt{N-2}\hfill & \hfill {g}_{21}\hfill \\ \hfill {g}_{01}\sqrt{N-2}\hfill & \hfill {\varepsilon }_{2}^{({0}^{N-2}{1}^{2})}\hfill & \hfill 0\hfill \\ \hfill {g}_{12}\hfill & \hfill 0\hfill & \hfill {\varepsilon }_{2}^{({0}^{N-1}{2}^{1})}\hfill \end{matrix}\right)\end{equation} \tag{ 63b }$

$\begin{equation}{\mathbf{H}}_{2}^{[N-2,2]}={\mathbf{1}}_{N(N-3)/2}\otimes \left(\begin{matrix}\hfill {\varepsilon }_{2}^{{0}^{N-2}{1}^{2}}\hfill \end{matrix}\right).\end{equation} \tag{ 63c }$

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 define

$\begin{equation}{\hat{a}}_{i}={\hat{\sigma }}_{v-1,v}^{(i)}\quad v{ >}0,\end{equation} \tag{ 64 }$

and the eigenstates can be represented as excitations of the polariton modes:

$\begin{equation}\left\vert v(\mathrm{U}\mathrm{P})v(\mathrm{L}\mathrm{P})v({\mathrm{D}}_{1})\dots v({\mathrm{D}}_{N-1})\right\rangle ={({a}_{+}^{{\dagger}})}^{v(\mathrm{U}\mathrm{P})}{({a}_{-}^{{\dagger}})}^{v(\mathrm{L}\mathrm{P})}\prod\limits _{k=1}^{N-1}{({\hat{a}}_{\mathrm{D}(k)})}^{v({\mathrm{D}}_{k})}\left\vert 0\right\rangle ,\end{equation} \tag{ 65 }$

where

$\begin{equation}{\hat{a}}_{{\pm}}^{{\dagger}}=\frac{1}{\sqrt{2{{\Omega}}_{10}}}\left({\pm}\sqrt{{{\Omega}}_{10}{\pm}{\Delta}}{\hat{a}}_{0}^{{\dagger}}+\sqrt{\frac{{{\Omega}}_{10}\mp {\Delta}}{N}}\sum\limits _{i=1}^{N}{\hat{a}}_{i}^{{\dagger}}\right),\end{equation} \tag{ 66a }$

$\begin{equation}{\hat{a}}_{\mathrm{D}(k)}=\frac{1}{\sqrt{k}}\left(\sqrt{k-1}{\hat{a}}_{k}^{{\dagger}}-\frac{1}{\sqrt{k-1}}\sum\limits _{i=1}^{k-1}{\hat{a}}_{i}^{{\dagger}}\right).\end{equation} \tag{ 66b }$

The quantities Δ = ω − (ɛ₁ − ɛ₀)/ℏ, and ${{\Omega}}_{10}={[{{\Delta}}^{2}+4{({g}_{10}/\hslash )}^{2}N]}^{1/2}$ are the detuning between the EM mode and the 0 → 1 transition of the emitters, and Rabi frequency, respectively. 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 (65) can still be consistent with some of the features displayed by these states.

The method introduced in this work requires only of the energy spectrum of a single emitter and is therefore compatible with any approach to determine or approximate it. Since it is a good approximation for covalent interactions, we introduce anharmonicities by considering each emitter as a Morse oscillator (MO) [76], i.e. the single-emitter Hamiltonians include a potential energy function of the form

$\begin{equation}V(R)={D}_{\mathrm{e}}{\left(1-{\mathrm{e}}^{-\sqrt{8{D}_{\mathrm{e}}}{\chi }_{\mathrm{e}}(R-{R}_{\mathrm{e}})/\hslash }\right)}^{2},\end{equation} \tag{ 67 }$

where R is the mass-scaled length of the oscillator with value at equilibrium R_e, D_e = V(∞) − V(R_e) is the dissociation energy, and χ_e is the so-called anharmonicity constant. The corresponding eigenenergies are given by

$\begin{equation}{\varepsilon }_{v}=4{D}_{\mathrm{e}}{\chi }_{\mathrm{e}}\left(v+\frac{1}{2}\right)\left[1-{\chi }_{\mathrm{e}}\left(v+\frac{1}{2}\right)\right]\quad {\varepsilon }_{v}{\leqslant}{D}_{\mathrm{e}},\end{equation} \tag{ 68 }$

with the implication that the number of bound-states is ⌊r_MO⌋, where

$\begin{equation}{r}_{\text{MO}}=\frac{1-{\chi }_{\mathrm{e}}}{2{\chi }_{\mathrm{e}}};\end{equation} \tag{ 69 }$

therefore, the potential becomes harmonic as D_e → ∞ (χ_e → 0) [77].

From (67), it follows that, in the neighborhood of the equilibrium length, the emitters oscillate with frequency

$\begin{equation}{\omega }_{\text{MO}}=\frac{4{D}_{\mathrm{e}}{\chi }_{\mathrm{e}}}{\hslash }.\end{equation} \tag{ 70 }$

On the other hand, the difference in energies of contiguous levels fulfill

$\begin{equation}{\varepsilon }_{v+1}-{\varepsilon }_{v}=\hslash {\omega }_{\text{MO}}\left[1-2{\chi }_{\mathrm{e}}\left(v+1\right)\right];\end{equation} \tag{ 71 }$

therefore, χ_e characterizes the mechanical anharmonicity introduced by this model. There is also an electric anharmonicity, which stems from the fact that [78]

$\begin{equation}\frac{{g}_{v+1,v}^{\text{MO}}}{{g}_{v+1,v}^{\text{HO}}}=\frac{2}{2({r}_{\text{MO}}-v)-1}{\left[\frac{({r}_{\text{MO}}-v)({r}_{\text{MO}}-v-1)}{{\chi }_{\mathrm{e}}(2{r}_{\text{MO}}-v)}\right]}^{1/2}.\end{equation} \tag{ 72 }$

It is pertinent to remark that the dipole moment operator induces transitions between all the levels in an 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, $\hslash {{\Omega}}_{10}=2{g}_{10}\sqrt{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.

Zoom In Zoom Out Reset image size

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 $\left\vert 3(\mathrm{U}\mathrm{P})\right\rangle$ is unique in regards of its spectral configuration; therefore, it belongs to λ = [N]. On the other hand, states of the form $\left\vert 1(\mathrm{U}\mathrm{P})1(\mathrm{L}\mathrm{P})1({\mathrm{D}}_{k})\right\rangle$ have multiplicity of N − 1 and symmetry λ = [N − 1, 1]. Also, from the (N − 1)(N − 2)/2 states of the form $\left\vert 1(\mathrm{L}\mathrm{P})1({\mathrm{D}}_{k})1({\mathrm{D}}_{{k}^{\prime }})\right\rangle$, one has symmetry λ = [N] while the remaining N(N − 3)/2 are of the λ = [N − 2, 2] kind. In contrast, when anharmonicities are involved, the eigenstates can be labeled according to harmonic states of reference with the appropriate energetics. As a consequence, a state with labels $\left\vert 1(\mathrm{L}\mathrm{P})2({\mathrm{D}}_{k})\right\rangle$, which accounts for BEs with $\tilde {\mu }={0}^{N-1}{2}^{1}$, now has symmetry λ = [N] despite belonging to a multiplet of size N − 1. The lift in degeneracy, together with the diversity in slopes with which the eigenenergies depend on the coupling strength generate situations in which the ordering of levels changes for different coupling intensities. This might result in the inversion of labeling of spectroscopic signals for samples with different concentrations of emitters, which, in turn, may open the door for coherent control strategies based on reactive systems. Furthermore, conditions of near degeneracy modulated by light–matter coupling can result in enhancement of multiphoton absorption, as suggested in [32] for the doubly excited manifold.

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.

Zoom In Zoom Out Reset image size

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, ${\hat{v}}_{0}={\hat{a}}_{0}^{{\dagger}}{\hat{a}}_{0}$, which measures the photon contents of a given state. In the symmetrized basis, this operator is diagonal and has the form

$\begin{equation}{\hat{v}}_{0}=\sum\limits _{\tilde {\mu }:\boldsymbol{\mu }{\trianglelefteq}\boldsymbol{\lambda }}\sum\limits _{i=1}^{\mathrm{dim}({S}^{\boldsymbol{\lambda }})}\sum\limits _{j=1}^{{K}_{\boldsymbol{\lambda }\boldsymbol{\mu }}}{v}_{0}^{({n}_{\text{exc}},\tilde {\mu })}\left\vert {n}_{\text{exc}},\tilde {\mu },\boldsymbol{\lambda };{\mathbf{y}}_{i},{\mathbf{w}}_{j}\right\rangle \left\langle {n}_{\text{exc}},\tilde {\mu },\boldsymbol{\lambda };{\mathbf{y}}_{i},{\mathbf{w}}_{j}\right\vert .\end{equation} \tag{ 73 }$

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.

Zoom In Zoom Out Reset image size

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.

JACGA thanks Kai Schwennicke and Stephan van den Wildenberg for their useful comments and insights. JACGA, RFR, and JYZ acknowledge funding support from the Air Force Office of Scientific Research award FA9550-18-1-0289. JACGA acknowledges initial support of the UC-MEXUS-CONACYT graduate scholarship ref. 235273/472318.

No new data were created or analysed in this study.

After recognizing n_exc as a conserved quantity, the explicit forms of ${\hat{H}}_{{n}_{\text{exc}}}$ for the groundstate, as well as the first, second, and third excitation manifolds are, respectively, ${\hat{H}}_{0}(N)=\left(\frac{\hslash \omega }{2}+N{\varepsilon }_{0}\right)\left\vert 0\right\rangle \left\langle 0\right\vert$

$\begin{align}{\hat{H}}_{1}(N)& =\left(\frac{3}{2}\hslash \omega +N{\varepsilon }_{0}\right)\left\vert {1}_{0}\right\rangle \left\langle {1}_{0}\right\vert +\left(\frac{\hslash \omega }{2}+{\varepsilon }_{1}+(N-1){\varepsilon }_{0}\right)\sum\limits _{i=1}^{N}\left\vert {1}_{i}\right\rangle \left\langle {1}_{i}\right\vert \\ & \quad +\sum\limits _{i=1}^{N}\left({g}_{01}\left\vert {1}_{0}\right\rangle \left\langle {1}_{i}\right\vert +{g}_{10}\left\vert {1}_{i}\right\rangle \left\langle {1}_{0}\right\vert \right),\end{align} \tag{ 1.1a }$

$\begin{align}{\hat{H}}_{2}(N)& =\left(\frac{5}{2}\hslash \omega +N{\varepsilon }_{0}\right)\left\vert {2}_{0}\right\rangle \left\langle {2}_{0}\right\vert +\left(\frac{3}{2}\hslash \omega +{\varepsilon }_{1}+(N-1){\varepsilon }_{0}\right)\sum\limits _{i=1}^{N}\left\vert {1}_{0}{1}_{i}\right\rangle \left\langle {1}_{0}{1}_{i}\right\vert \\ & \quad +\left(\frac{\hslash \omega }{2}+2{\varepsilon }_{1}+(N-2){\varepsilon }_{0}\right)\sum\limits _{i=1}^{N-1}\sum\limits _{j=i+1}^{N}\left\vert {1}_{i}{1}_{j}\right\rangle \left\langle {1}_{i}{1}_{j}\right\vert \\ & \quad +\left(\frac{\hslash \omega }{2}+{\varepsilon }_{2}+(N-1){\varepsilon }_{0}\right)\sum\limits _{i=1}^{N}\left\vert {2}_{i}\right\rangle \left\langle {2}_{i}\right\vert \end{align} \tag{ 1.1b }$

$\begin{align}& \quad +\sum\limits _{i=1}^{N}\left[\sqrt{2}{g}_{01}\left\vert {2}_{0}\right\rangle \left\langle {1}_{0}{1}_{i}\right\vert +\left\vert {1}_{0}{1}_{i}\right\rangle \left({g}_{12}\left\langle {2}_{i}\right\vert +{g}_{01}\sum\limits _{j\ne i}\left\langle {1}_{i}{1}_{j}\right\vert \right)+\mathrm{H}.\mathrm{c}.\right],\\ {\hat{H}}_{3}(N)& =\left(\frac{7}{2}\hslash \omega +N{\varepsilon }_{0}\right)\left\vert {3}_{0}\right\rangle \left\langle {3}_{0}\right\vert +\left(\frac{5}{2}\hslash \omega +{\varepsilon }_{1}+(N-1){\varepsilon }_{0}\right)\sum\limits _{i=1}^{N}\left\vert {2}_{0}{1}_{i}\right\rangle \left\langle {2}_{0}{1}_{i}\right\vert \\ & \quad +\left(\frac{3}{2}\hslash \omega +{\varepsilon }_{2}+(N-1){\varepsilon }_{0}\right)\sum\limits _{i=1}^{N}\left\vert {1}_{0}{2}_{i}\right\rangle \left\langle {1}_{0}{2}_{i}\right\vert \\ & \quad +\left(\frac{3}{2}\hslash \omega +2{\varepsilon }_{1}+(N-2){\varepsilon }_{0}\right)\sum\limits _{i=1}^{N-1}\sum\limits _{j=i+1}^{N}\left\vert {1}_{0}{1}_{i}{1}_{j}\right\rangle \left\langle {1}_{0}{1}_{i}{1}_{j}\right\vert \\ & \quad +\left(\frac{\hslash \omega }{2}+{\varepsilon }_{3}+(N-1){\varepsilon }_{0}\right)\sum\limits _{i=1}^{N}\left\vert {3}_{i}\right\rangle \left\langle {3}_{i}\right\vert \\ & \quad +\left(\frac{\hslash \omega }{2}+{\varepsilon }_{2}+{\varepsilon }_{1}+(N-2){\varepsilon }_{0}\right)\sum\limits _{i=1}^{N}\sum\limits _{j\ne i}\left\vert {2}_{i}{1}_{j}\right\rangle \left\langle {2}_{i}{1}_{j}\right\vert \\ & \quad +\left(\frac{\hslash \omega }{2}+3{\varepsilon }_{1}+(N-3){\varepsilon }_{0}\right)\sum\limits _{i=1}^{N-2}\sum\limits _{j=i+1}^{N-1}\sum\limits _{k=j+1}^{N}\left\vert {1}_{i}{1}_{j}{1}_{k}\right\rangle \left\langle {1}_{i}{1}_{j}{1}_{k}\right\vert \\ & \quad +\sum\limits _{i=1}^{N}\left[\sqrt{3}{g}_{01}\left\vert {3}_{0}\right\rangle \left\langle {2}_{0}{1}_{i}\right\vert +\sqrt{2}\left\vert {2}_{0}{1}_{i}\right\rangle \left({g}_{12}\left\langle {1}_{0}{2}_{i}\right\vert +{g}_{01}\sum\limits _{j\ne i}\left\langle {1}_{0}{1}_{i}{1}_{j}\right\vert \right)+\mathrm{H}.\mathrm{c}.\right]\end{align} \tag{ 1.1c }$

$\begin{align}& \quad +\sum\limits _{i=1}^{N-1}\sum\limits _{j=i+1}^{N}\left\{\left\vert {1}_{0}{1}_{i}{1}_{j}\right\rangle \left[{g}_{12}\left(\left\langle {2}_{i}{1}_{j}\right\vert +\left\langle {2}_{j}{1}_{i}\right\vert \right)+{g}_{01}\sum\limits _{k\ne i,j}\left\langle {1}_{i}{1}_{j}{1}_{k}\right\vert \right]+\mathrm{H}.\mathrm{c}.\right\}\\ & \quad +\sum\limits _{i=1}^{N}\left[\left\vert {1}_{0}{2}_{i}\right\rangle \left({g}_{23}\left\langle {3}_{i}\right\vert +{g}_{01}\sum\limits _{j\ne i}\left\langle {2}_{i}{1}_{j}\right\vert \right)+\mathrm{H}.\mathrm{c}.\right],\end{align} \tag{ 1.1d }$

where $\left\vert 0\right\rangle =\left\vert {0}_{0}{0}_{1}\dots {0}_{N}\right\rangle$ is the global ground state and H.c. stands for Hermitian conjugate.

In general, these can be written in terms of the spectral configuration as

$\begin{align}{\hat{H}}_{{n}_{\text{exc}}}(N)& =\sum\limits _{\tilde {\mu }}\sum\limits _{y=1}^{\vert \tilde {\mu }\vert }{\varepsilon }_{\tilde {\mu }}\left\vert {n}_{\text{exc}},\tilde {\mu };y\right\rangle \left\langle {n}_{\text{exc}},\tilde {\mu };y\right\vert \\ & \quad +\sum\limits _{\tilde {\mu }}\left({g}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}{\hat{J}}_{+}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}{\hat{a}}_{0}\left\vert {n}_{\text{exc}},\tilde {\mu };y\right\rangle \left\langle {n}_{\text{exc}},\tilde {\mu };y\right\vert +\mathrm{H}.\mathrm{c}.\right),\end{align} \tag{ 1.2 }$

where the first sum runs over the possible |n_exc| spectral configurations within the respective manifold, and

$\begin{equation}\vert \tilde {\mu }\vert =\frac{N!}{{\prod }_{i=1}^{r}{\mu }_{i}!}.\end{equation} \tag{ 1.3 }$

This section includes useful illustrations for some of the tools and concepts introduced throughout the main text.

B.1. Spectral configurations and Young tableaux

Tables B1–B3 display examples of the relations between the spectral configuration and the main partitions to consider, as well as the corresponding Young diagrams, standard Young tableaux, and semistandard Young tableaux.

Table B1. Examples of spectral configurations, partitions, and Young diagrams corresponding to selected PBEs with ${n}_{\text{exc}}=0$.

PBE	$\tilde {\mu }$	μ	ν	Young diagram a
$\left\vert 0\right\rangle$	0N	[N]	[0]
$\left\vert {1}_{i}\right\rangle$	0N−111	[N − 1, 1]	[1]
$\left\vert {1}_{i}{1}_{j}\right\rangle$	0N−212	[N − 2, 2]	[1, 1]
$\left\vert {2}_{i}\right\rangle$	0N−121	[N − 1, 1]	[2]
$\left\vert {2}_{i}{4}_{j}{5}_{k}\right\rangle$	0N−3214151	[N − 3, 1, 1, 1]	[5, 4, 2]
$\left\vert {1}_{i}{1}_{j}{2}_{k}{2}_{\ell }{2}_{m}{4}_{n}\right\rangle$	122341	[3, 2, 1]	[4, 2, 2, 2, 1, 1]

^a N = 6.

Table B2. Examples of standard Young tableaux and their associated indices for N = 6.

SYT	λ	y
	[6]	0
	[5, 1]	6
	[4, 2]	(3, 5)
	[3, 3]	(3, 5, 6)
	[4, 1, 1]	(3; 4)
	[3, 2, 1]	(2, 5; 6)

Table B3. Examples of semistandard Young tableaux and their associated indices for N = 6 and $\tilde {\mu }={0}^{3}{1}^{2}{2}^{1}$.

SSYT	λ	w
	[6]	0
	[5, 1]	11
	[5, 1]	21
	[4, 2]	1121
	[3, 3]	1221
	[4, 1, 1]	(11, 21)
	[3, 2, 1]	(12, 21)

B.2. Determination of hook-lengths

Let us consider the Young diagram for λ = [6, 4, 3, 2, 1], the hook of the cell with coordinates (2, 2) is

and the hook-lengths for each cell are

B.3. Reduction of representations of S₃

The six elements of S₃ in the basis of PBEs with $\tilde {\mu }={1}^{1}$ are

$\begin{equation}\begin{aligned}\hat{\pi }()=\left\vert {1}_{1}\right\rangle \left\langle {1}_{1}\right\vert +\left\vert {1}_{2}\right\rangle \left\langle {1}_{2}\right\vert +\left\vert {1}_{3}\right\rangle \left\langle {1}_{3}\right\vert ,\\ \hat{\pi }(12)=\left\vert {1}_{1}\right\rangle \left\langle {1}_{2}\right\vert +\left\vert {1}_{2}\right\rangle \left\langle {1}_{1}\right\vert +\left\vert {1}_{3}\right\rangle \left\langle {1}_{3}\right\vert ,\\ \hat{\pi }(23)=\left\vert {1}_{1}\right\rangle \left\langle {1}_{1}\right\vert +\left\vert {1}_{2}\right\rangle \left\langle {1}_{3}\right\vert +\left\vert {1}_{3}\right\rangle \left\langle {1}_{2}\right\vert ,\\ \hat{\pi }(31)=\left\vert {1}_{1}\right\rangle \left\langle {1}_{3}\right\vert +\left\vert {1}_{2}\right\rangle \left\langle {1}_{2}\right\vert +\left\vert {1}_{3}\right\rangle \left\langle {1}_{1}\right\vert ,\\ \hat{\pi }(123)=\left\vert {1}_{1}\right\rangle \left\langle {1}_{3}\right\vert +\left\vert {1}_{2}\right\rangle \left\langle {1}_{1}\right\vert +\left\vert {1}_{3}\right\rangle \left\langle {1}_{2}\right\vert ,\\ \hat{\pi }(132)=\left\vert {1}_{1}\right\rangle \left\langle {1}_{2}\right\vert +\left\vert {1}_{2}\right\rangle \left\langle {1}_{3}\right\vert +\left\vert {1}_{3}\right\rangle \left\langle {1}_{1}\right\vert .\end{aligned}\end{equation} \tag{ 2.1 }$

Let us investigate the consequence of changing the basis to that of the SABEs.

Starting from the state $\left\vert 1,{0}^{3},[3];0,0\right\rangle$, the action of the operator ${\hat{a}}_{0}{\hat{J}}_{+}^{(1,0)}$ yields

$\begin{equation}\left\vert 1,{0}^{2}{1}^{1},[3];0,0\right\rangle =\frac{1}{\sqrt{3}}\left(\left\vert {1}_{1}\right\rangle +\left\vert {1}_{2}\right\rangle +\left\vert {1}_{3}\right\rangle \right).\end{equation} \tag{ 2.2 }$

The remaining two SABEs with $\tilde {\mu }={0}^{2}{1}^{1}$ are

$\begin{equation}\left\vert 1,{0}^{2}{1}^{1},[2,1];2,{1}^{1}\right\rangle =\frac{\left\vert {1}_{1}\right\rangle -\left\vert {1}_{2}\right\rangle }{\sqrt{2}},\end{equation} \tag{ 2.3a }$

and

$\begin{equation}\left\vert 1,{0}^{2}{1}^{1},[2,1];3,{1}^{1}\right\rangle =\frac{-\left\vert {1}_{1}\right\rangle -\left\vert {1}_{2}\right\rangle +2\left\vert {1}_{3}\right\rangle }{\sqrt{6}}.\end{equation} \tag{ 2.3b }$

The permutation operators thus become

$\begin{equation}\begin{aligned}\hat{\pi }()& =\left\vert [3];0\right\rangle \left\langle [3];0\right\vert +\left\vert [2,1];2\right\rangle \left\langle [2,1];2\right\vert +\left\vert [2,1];3\right\rangle \left\langle [2,1];3\right\vert ,\\ \hat{\pi }(12)& =\left\vert [3];0\right\rangle \left\langle [3];0\right\vert -\left\vert [2,1];2\right\rangle \left\langle [2,1];2\right\vert +\left\vert [2,1];3\right\rangle \left\langle [2,1];3\right\vert ,\\ \hat{\pi }(23)& =\left\vert [3];0\right\rangle \left\langle [3];0\right\vert +\frac{1}{2}\left(\left\vert [2,1];2\right\rangle \left\langle [2,1];2\right\vert -\left\vert [2,1];3\right\rangle \left\langle [2,1];3\right\vert \right)\\ & \quad +\frac{\sqrt{3}}{2}\left(\left\vert [2,1];2\right\rangle \left\langle [2,1];3\right\vert +\left\vert [2,1];3\right\rangle \left\langle [2,1];2\right\vert \right),\\ \hat{\pi }(31)& =\left\vert [3];0\right\rangle \left\langle [3];0\right\vert +\frac{1}{2}\left(\left\vert [2,1];2\right\rangle \left\langle [2,1];2\right\vert -\left\vert [2,1];3\right\rangle \left\langle [2,1];3\right\vert \right)\\ & \quad -\frac{\sqrt{3}}{2}\left(\left\vert [2,1];2\right\rangle \left\langle [2,1];3\right\vert +\left\vert [2,1];3\right\rangle \left\langle [2,1];2\right\vert \right),\\ \hat{\pi }(123)& =\left\vert [3];0\right\rangle \left\langle [3];0\right\vert -\frac{1}{2}\left(\left\vert [2,1];2\right\rangle \left\langle [2,1];2\right\vert +\left\vert [2,1];3\right\rangle \left\langle [2,1];3\right\vert \right)\\ & \quad -\frac{\sqrt{3}}{2}\left(\left\vert [2,1];2\right\rangle \left\langle [2,1];3\right\vert -\left\vert [2,1];3\right\rangle \left\langle [2,1];2\right\vert \right),\\ \hat{\pi }(132)& =\left\vert [3];0\right\rangle \left\langle [3];0\right\vert -\frac{1}{2}\left(\left\vert [2,1];2\right\rangle \left\langle [2,1];2\right\vert +\left\vert [2,1];3\right\rangle \left\langle [2,1];3\right\vert \right)\\ & \quad +\frac{\sqrt{3}}{2}\left(\left\vert [2,1];2\right\rangle \left\langle [2,1];3\right\vert -\left\vert [2,1];3\right\rangle \left\langle [2,1];2\right\vert \right),\end{aligned}\end{equation} \tag{ 2.4 }$

where the notation was simplified by making implicit the common indices n_exc = 1 and $\tilde {\mu }={0}^{2}{1}^{1}$, as well as by omitting the label w. As it can be seen, there are two clear independent subspaces labeled by [1–3], respectively.

The SABEs obtained through application of the collective excitations, ${\hat{J}}_{+}^{(u,v)}$, and Gram–Schmidt orthogonalization have been discussed in the literature under the name of Schur–Weyl states [42, 79–81], or Gelfand–Tsetlin states [82–85]. 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. The following formulae is proven to work for SABEs up to the triply excited manifold, but we have reason to believe that it prepares the ground for a method that also works with upper excitations. 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, ${\hat{a}}_{0}^{({\dagger})}$, are not included in the following discussion, i.e. the SABEs are denoted by $\left\vert \tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\rangle$.

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 ground-state and some of its collective excitations are

$\begin{equation}\left\vert {0}^{N},[N];0,0\right\rangle =\left\vert 0\right\rangle ,\end{equation} \tag{ 3.1a }$

$\begin{equation}\left\vert {0}^{N-1}{v}^{1},[N];0,0\right\rangle =\frac{1}{\sqrt{N}}\sum\limits _{i=1}^{N}\left\vert {v}_{i}\right\rangle ,\end{equation} \tag{ 3.1b }$

$\begin{equation}\left\vert {0}^{N-2}{v}^{2},[N];0,0\right\rangle ={\left(\genfrac{}{}{0.0pt}{}{N}{2}\right)}^{-1/2}\sum\limits _{i=1}^{N-1}\sum\limits _{j=i+1}^{N}\left\vert {v}_{i}{v}_{j}\right\rangle ,\end{equation} \tag{ 3.1c }$

$\begin{equation}\left\vert {0}^{N-3}{v}^{3},[N];0,0\right\rangle ={\left(\genfrac{}{}{0.0pt}{}{N}{3}\right)}^{-1/2}\sum\limits _{i=1}^{N-2}\sum\limits _{j=i+1}^{N-1}\sum\limits _{k=j+1}^{N}\left\vert {v}_{i}{v}_{j}{v}_{k}\right\rangle ,\end{equation} \tag{ 3.1d }$

$\begin{equation}\left\vert {0}^{N-2}{u}^{1}{v}^{1},[N];0,0\right\rangle =\frac{1}{\sqrt{N(N-1)}}\sum\limits _{i=1}^{N-1}\sum\limits _{j=i+1}^{N}\left(\left\vert {u}_{i}{v}_{j}\right\rangle +\left\vert {v}_{i}{u}_{j}\right\rangle \right),\end{equation} \tag{ 3.1e }$

where λ = [N] and y = 0 indicate the SYT 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₁, i₂, ..., i_n ), let us define the vandermonde-like matrix $\hat{\mathbf{V}}({i}_{1},{i}_{2},\dots ,{i}_{n})$ with elements

$\begin{equation}{\left[\hat{\mathbf{V}}({i}_{1},{i}_{2},\dots ,{i}_{n})\right]}_{\alpha ,\beta }={\hat{\sigma }}_{\alpha -1,0}^{({i}_{\beta })}.\end{equation} \tag{ 3.2 }$

A given SYT can be associated to Young operators of the form

$\begin{equation}\hat{Y}(\boldsymbol{\lambda },\mathbf{y})=\sum\limits _{\psi }\prod\limits _{j=1}^{{\lambda }_{1}}\enspace \mathrm{det}\left(\hat{\mathbf{V}}({\mathrm{c}\mathrm{o}\mathrm{l}}_{j}[{\mathbf{y}}_{\psi }^{\prime }])\right),\end{equation} \tag{ 3.3 }$

where col_j (A) extracts the jth column of array A, and ${\mathbf{y}}_{\psi }^{\prime }$ 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 is

$\begin{align}\hat{Y}\left([3,2,1],(2,5;6)\right)& =\mathrm{det}\left(\hat{\mathbf{V}}(1,2,6)\right)\mathrm{det}\left(\hat{\mathbf{V}}(3,5)\right)\mathrm{det}\left(\hat{\mathbf{V}}(4)\right)\\ & \quad +\mathrm{det}\left(\hat{\mathbf{V}}(1,2,6)\right)\mathrm{det}\left(\hat{\mathbf{V}}(3)\right)\mathrm{det}\left(\hat{\mathbf{V}}(4,5)\right)\\ & \quad +\mathrm{det}\left(\hat{\mathbf{V}}(1,2)\right)\mathrm{det}\left(\hat{\mathbf{V}}(3,5,6)\right)\mathrm{det}\left(\hat{\mathbf{V}}(4)\right)\\ & \quad +\mathrm{det}\left(\hat{\mathbf{V}}(1,2)\right)\mathrm{det}\left(\hat{\mathbf{V}}(3)\right)\mathrm{det}\left(\hat{\mathbf{V}}(4,5,6)\right).\end{align} \tag{ 3.6 }$

If μ = λ , the wavefunction $\left\vert \tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\rangle$ is proportional to $\hat{Y}(\boldsymbol{\lambda },\mathbf{y})\left\vert {0}^{N},[N];0,0\right\rangle$. Notice that the constraint μ = λ uniquely defines an SSYT, and therefore a w.

The states orthogonal to $\left\vert {0}^{N-1}{1}^{1},[N];0,0\right\rangle$ with the same $\tilde {\mu }$ can be written as $\left\vert {0}^{N-1}{1}^{1},[N-1,1];k,{1}^{1}\right\rangle$, with 2 ⩽ k ⩽ N. The corresponding SYT are

which have associated Young operators of the form

$\begin{equation}\hat{Y}([N-1,1],k)=\sum\limits _{i=1}^{k-1}\prod\limits _{j=1}^{i-1}{\hat{\sigma }}_{0,0}^{(j)}\cdot \left\vert \begin{matrix}\hfill {\hat{\sigma }}_{0,0}^{(i)}\hfill & \hfill {\hat{\sigma }}_{0,0}^{(k)}\hfill \\ \hfill {\hat{\sigma }}_{1,0}^{(i)}\hfill & \hfill {\hat{\sigma }}_{1,0}^{(k)}\hfill \end{matrix}\right\vert \cdot \prod\limits _{{j}^{\prime }=i+1}^{N}{\hat{\sigma }}_{0,0}^{({j}^{\prime })}.\end{equation} \tag{ 3.7 }$

Therefore, the corresponding SABEs are

$\begin{equation}\left\vert {0}^{N-1}{1}^{1},[N-1,1];k,{1}^{1}\right\rangle =\frac{1}{\sqrt{k(k-1)}}\left(-\sum\limits _{i=1}^{k-1}\left\vert {1}_{i}\right\rangle +(k-1)\left\vert {1}_{k}\right\rangle \right).\end{equation} \tag{ 3.8 }$

The functions $\left\vert {0}^{N-1}{v}^{1},[N-1,1];k,{v}^{1}\right\rangle$ for v > 1 have an equivalent form. Furthermore, the action of the collective excitation ${\hat{J}}_{+}^{(1,0)}$ yields

$\begin{equation}\left\vert {0}^{N-2}{1}^{2},[N-1,1];k,{1}^{1}\right\rangle ={\left[\left(N-2\right)k\left(k-1\right)\right]}^{-1/2}\sum\limits _{i=1}^{k}\sum\limits _{j=i+1}^{N}{c}_{ij}^{(k)}\left\vert {1}_{i}{1}_{j}\right\rangle ,\end{equation} \tag{ 3.9a }$

where

$\begin{equation}{c}_{ij}^{(k)}=\begin{cases}-2\quad & \quad \text{for}\enspace 1{\leqslant}i{< }k-1\;\text{and}\;1{< }j{\leqslant}k-1,\\ k-2\quad & \quad \text{for}\enspace 1{\leqslant}i{\leqslant}k-1\;\text{and}\;j=k,\\ -1\quad & \quad \text{for}\enspace 1{\leqslant}i{\leqslant}k-1\;\text{and}\;k{< }j{\leqslant}N,\\ k-1\quad & \quad \text{for}\enspace i=k\;\text{and}\;k{< }j{\leqslant}N\end{cases}.\end{equation} \tag{ 3.9b }$

Following the same reasoning, the remaining states with $\tilde {\mu }={0}^{N-2}{1}^{2}$ that are orthogonal to all the states above are

$\begin{equation}\left\vert {0}^{N-2}{1}^{2},[N-2,2];(k,\ell ),{1}^{2}\right\rangle ={\left[k(k-1)(\ell -2)(\ell -3)\right]}^{-1/2}\sum\limits _{i=1}^{k}\sum\limits _{j=i+1}^{\ell }{c}_{ij}^{(k,\ell )}\left\vert {1}_{i}{1}_{j}\right\rangle ,\end{equation} \tag{ 3.10a }$

where

$\begin{equation}{c}_{ij}^{(k,\ell )}=\begin{cases}2\quad & \quad \text{for}\enspace 1{\leqslant}i{< }k-1\;\text{and}\;1{< }j{\leqslant}k-1,\\ -(k-2)\quad & \quad \text{for}\enspace 1{\leqslant}i{\leqslant}k-1\;\text{and}\;j=k,\\ 1\quad & \quad \text{for}\enspace 1{\leqslant}i{\leqslant}k-1\;\text{and}\;k{< }j{< }\ell ,\\ -(\ell -3)\quad & \quad \text{for}\enspace 1{\leqslant}i{\leqslant}k-1\;\text{and}\;j=\ell ,\\ -(k-1)\quad & \quad \text{for}\enspace i=k\;\text{and}\;k{< }j{< }\ell ,\\ (k-1)(\ell -3)\quad & \quad \text{for}\enspace i=k\;\text{and}\;j=\ell ,\\ 0\quad & \quad \text{otherwise}\end{cases}.\end{equation} \tag{ 3.10b }$

Tables C1 and C2 show, up to the triply excited manifold, the excitation operators, $\hat{O}\left\{\tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\}$, that produce the states with the same labels when applied to the global groundstate of the emitters, i.e.

$\begin{equation}\hat{O}\left\{\tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\}\left\vert 0\right\rangle ={\mathcal{N}}^{-2}\left\vert \tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\rangle .\end{equation} \tag{ 3.11 }$

The nested structure of the operators responds to the fact that the higher states are obtained through either collective excitations, or Gram–Schmidt orthogonalization of the lower ones.

Table C1. SABEs in terms of PBEs for spectral configurations with at most three excitations in the emitters.

$\tilde {\mu }$	λ	y	w	$\hat{O}\left\{\tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\}$	${\mathcal{N}}^{-2}$	Notes
0N	[N]	0	0	1	1
0N−1 v1	[N]	0	0	$\sum _{i=1}^{N}{\hat{\sigma }}_{v,0}^{(i)}$	N	v ∈{1, 2, 3}
	[N − 1, 1]	k	v1	$(k-1){\hat{\sigma }}_{v,0}^{(k)}-\hat{O}\left\{{0}^{k-2}{1}^{1},[k-1];0,0\right\}$	k(k − 1)	2⩽k⩽N
0N−212	[N]	0	0	$2\sum _{i=1}^{N-1}\sum _{j=i+1}^{N}{\hat{\sigma }}_{1,0}^{(i)}{\hat{\sigma }}_{1,0}^{(j)}$	2N(N − 1)
	[N − 1, 1]	k	11	$\begin{aligned}(k-2)& {\hat{\sigma }}_{1,0}^{(k)}\hat{O}\left\{{0}^{k-2}{1}^{1},[k-1];0,0\right\}-\hat{O}\left\{{0}^{k-3}{1}^{2},[k-1];0,0\right\}\\ +& \sum _{j=k+1}^{N}{\hat{\sigma }}_{1,0}^{(j)}\hat{O}\left\{{0}^{k-1}{1}^{1},[k-1,1];k,{1}^{1}\right\}\end{aligned}$	k(k − 1)(N − 2)	2⩽k⩽N
	[N − 2, 2]	(k, ℓ)	12	$(\ell -3){\hat{\sigma }}_{1,0}^{(\ell )}\hat{O}\left\{{0}^{k-1}{1}^{1},[k-1,1];k,{1}^{1}\right\}-\hat{O}\left\{{0}^{\ell -3}{1}^{2},[\ell -2,1];k,{1}^{1}\right\}$	k(k − 1)(ℓ − 2)(ℓ − 3)	$\begin{matrix}\hfill 4{\leqslant}\ell {\leqslant}N\hfill \\ \hfill 2{\leqslant}k{< }\ell \hfill \end{matrix}$
0N−313	[N]	0	0	$6\sum _{h=1}^{N-2}\sum _{i=h+1}^{N-1}\sum _{j=i+1}^{N}{\hat{\sigma }}_{1,0}^{(h)}{\hat{\sigma }}_{1,0}^{(i)}{\hat{\sigma }}_{1,0}^{(j)}$	6N(N − 1)(N − 2)
	[N − 1, 1]	k	11	$\begin{aligned}(k-3)& {\hat{\sigma }}_{1,0}^{(k)}\hat{O}\left\{{0}^{k-3}{1}^{2},[k-1];0,0\right\}-\hat{O}\left\{{0}^{k-4}{1}^{3},[k-1];0,0\right\}\\ +& 2\sum _{h=k+1}^{N}{\hat{\sigma }}_{1,0}^{(h)}\left[(k-2){\hat{\sigma }}_{1,0}^{(k)}\hat{O}\left\{{0}^{k-2}{1}^{1},[k-1];0,0\right\}-\hat{O}\left\{{0}^{k-3}{1}^{2},[k-1];0,0\right\}\right]\\ & \quad +2\sum _{j=k+1}^{N-1}\sum _{h=j+1}^{N}{\hat{\sigma }}_{1,0}^{(j)}{\hat{\sigma }}_{1,0}^{(h)}\hat{O}\left\{{0}^{k-1}{1}^{1},[k-1,1];k,{1}^{1}\right\}\end{aligned}$	2k(k − 1)N(N − 1)	2⩽k⩽N
	[N − 2, 2]	(k, ℓ)	12	$\begin{aligned}(\ell -4)& {\hat{\sigma }}_{1,0}^{(\ell )}\hat{O}\left\{{0}^{\ell -3}{1}^{2},[\ell -2,1];k,1\right\}+\sum _{h=\ell +1}^{N}{\hat{\sigma }}_{1,0}^{(h)}\hat{O}\left\{{0}^{\ell -3}{1}^{2},[\ell -3,2];(k,\ell ),{1}^{2}\right\}\\ -& \hat{O}\left\{{0}^{\ell -4}{1}^{3},[\ell -2,1];k,{1}^{1}\right\}\end{aligned}$	k(k − 1)(ℓ − 2)(ℓ − 3)(N − 4)	$\begin{matrix}\hfill 4{\leqslant}\ell {\leqslant}N\hfill \\ \hfill 2{\leqslant}k{< }\ell \hfill \end{matrix}$
	[N − 3, 3]	(k, ℓ, m)	13	$(m-5){\hat{\sigma }}_{1,0}^{(m)}\hat{O}\left\{{0}^{\ell -3}{1}^{2},[\ell -3,2];(k,\ell ),{1}^{2}\right\}+\hat{O}\left\{{0}^{m-4}{1}^{3},[m-3,2];(k,\ell ),{1}^{2}\right\}$	k(k − 1)(ℓ − 2)(ℓ − 3)(m − 4)(m − 5)	$\begin{matrix}\hfill 6{\leqslant}m{\leqslant}N\hfill \\ \hfill 4{\leqslant}\ell {< }m\hfill \\ \hfill 2{\leqslant}k{< }\ell \hfill \end{matrix}$

Table C2. SABEs in terms of PBEs with $\tilde {\mu }={0}^{N-2}{1}^{1}{2}^{1}$.

$\tilde {\mu }$	λ	y	w	$\hat{O}\left\{\tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\}$	${\mathcal{N}}^{-2}$	Notes
0N−21121	[N]	0	0	$2\sum _{i=1}^{N-1}\sum _{j=i+1}^{N}\left({\hat{\sigma }}_{2,0}^{(i)}{\hat{\sigma }}_{1,0}^{(j)}+{\hat{\sigma }}_{2,0}^{(j)}{\hat{\sigma }}_{1,0}^{(i)}\right)$	4N(N − 1)
	[N − 1, 1]	k	11	$\begin{aligned}(k-2)& \left({\hat{\sigma }}_{2,0}^{(k)}\hat{O}\left\{{0}^{k-2}{1}^{1},[k-1];0,0\right\}+{\hat{\sigma }}_{1,0}^{(k)}\hat{O}\left\{{0}^{k-2}{2}^{1},[k-1];0,0\right\}\right)\\ -& \sum _{j=k+1}^{N}\left({\hat{\sigma }}_{1,0}^{(j)}\hat{O}\left\{{0}^{k-2}{2}^{1},[k-1];0,0\right\}+{\hat{\sigma }}_{2,0}^{(j)}\hat{O}\left\{{0}^{k-2}{1}^{1},[k-1];0,0\right\}\right)\\ & \quad +(k-1)\sum _{j=k+1}^{N}\left({\hat{\sigma }}_{2,0}^{(k)}{\hat{\sigma }}_{1,0}^{(j)}+{\hat{\sigma }}_{1,0}^{(k)}{\hat{\sigma }}_{2,0}^{(j)}\right)-\hat{O}\left\{{0}^{k-3}{1}^{1}{2}^{1},[k-1],0,0\right\}\end{aligned}$	2k(k − 1)(N − 2)	2⩽k⩽N
			21	$\begin{aligned}k& \left({\hat{\sigma }}_{2,0}^{(k)}\hat{O}\left\{{0}^{k-2}{1}^{1},[k-1];0,0\right\}-{\hat{\sigma }}_{1,0}^{(k)}\hat{O}\left\{{0}^{k-2}{2}^{1},[k-1];0,0\right\}\right)\\ & +\sum _{j=k+1}^{N}\left({\hat{\sigma }}_{2,0}^{(j)}\hat{O}\left\{{0}^{k-2}{1}^{1},[k-1];0,0\right\}-{\hat{\sigma }}_{1,0}^{(j)}\hat{O}\left\{{0}^{k-2}{2}^{1},[k-1];0,0\right\}\right)\\ & \quad +(k-1)\sum _{j=k+1}^{N}\left({\hat{\sigma }}_{2,0}^{(k)}{\hat{\sigma }}_{1,0}^{(j)}+{\hat{\sigma }}_{1,0}^{(k)}{\hat{\sigma }}_{2,0}^{(j)}\right)\end{aligned}$	2k(k − 1)N	2⩽k⩽N
	[N − 2, 2]	(k, ℓ)	1121	$\begin{aligned}(\ell -3)& \left({\hat{\sigma }}_{1,0}^{(\ell )}\hat{O}\left\{{0}^{k-1}{2}^{1},[k-1,1];k,{2}^{1}\right\}+{\hat{\sigma }}_{2,0}^{(\ell )}\hat{O}\left\{{0}^{k-1}{1}^{1},[k-1,1];k,{1}^{1}\right\}\right)\\ -& \hat{O}\left\{{0}^{\ell -3}{1}^{1}{2}^{1},[\ell -2,1];k,{1}^{1}\right\}\end{aligned}$	2k(k − 1)(ℓ − 2)(ℓ − 3)	$\begin{matrix}\hfill 4{\leqslant}\ell {\leqslant}N\hfill \\ \hfill 2{\leqslant}k{< }\ell \hfill \end{matrix}$
	[N − 2, 1, 1]	(k; ℓ)	(11, 21)	$\begin{aligned}(\ell -2)& \left({\hat{\sigma }}_{2,0}^{(k)}{\hat{\sigma }}_{1,0}^{(\ell )}+{\hat{\sigma }}_{1,0}^{(k)}{\hat{\sigma }}_{2,0}^{(\ell )}\right)+\sum _{i=k+1}^{\ell -1}\left[\left({\hat{\sigma }}_{2,0}^{(k)}-{\hat{\sigma }}_{2,0}^{(\ell )}\right){\hat{\sigma }}_{1,0}^{(i)}+\left({\hat{\sigma }}_{1,0}^{(\ell )}-{\hat{\sigma }}_{1,0}^{(k)}\right){\hat{\sigma }}_{2,0}^{(i)}\right]\\ +& \left({\hat{\sigma }}_{2,0}^{(k)}-{\hat{\sigma }}_{2,0}^{(\ell )}\right)\hat{O}\left\{{0}^{k-1}{1}^{1},[k-1];0,0\right\}+\left({\hat{\sigma }}_{1,0}^{(\ell )}-{\hat{\sigma }}_{1,0}^{(k)}\right)\hat{O}\left\{{0}^{k-1}{2}^{1},[k-1];0,0\right\}\end{aligned}$	2ℓ(ℓ − 1)	$\begin{matrix}\hfill 3{\leqslant}\ell {\leqslant}N\hfill \\ \hfill 2{\leqslant}k{< }\ell \hfill \end{matrix}$

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 (39).

First, we introduce a by-hand method to compute individual instances of (48). Taking advantage of the fact that ${n}_{\tilde {\mu }}({v}^{{\ast}}+s)+{n}_{\tilde {\mu }}({v}^{{\ast}})={n}_{{\tilde {\mu }}^{\prime }}({v}^{{\ast}}+s)+{n}_{{\tilde {\mu }}^{\prime }}({v}^{{\ast}})$, we can write

$\begin{equation}{L}_{{\tilde {\mu }}^{\prime }\tilde {\mu }}^{(\boldsymbol{\lambda })}={\left\{\left[{n}_{\tilde {\mu }}({v}^{{\ast}})-{\rho }_{J}\right]\left[{n}_{{\tilde {\mu }}^{\prime }}({v}^{{\ast}}+s)-{\rho }_{J}\right]\right\}}^{1/2},\end{equation} \tag{ 4.1 }$

where ${\rho }_{J}=[{n}_{\tilde {\mu }}({v}^{{\ast}}+s)+{n}_{\tilde {\mu }}({v}^{{\ast}})]/2-{J}_{{\tilde {\mu }}^{\prime }\tilde {\mu }}^{(\boldsymbol{\lambda })}$. 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 $\tilde {\mu }$ available for excitation with the number of emitters in ${\tilde {\mu }}^{\prime }$ 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 D1 illustrates this calculation for the couplings between $\tilde {\mu }={0}^{N-2}{1}^{2}$ and ${\tilde {\mu }}^{\prime }={0}^{N-3}{1}^{3}$.

Table D1. Contribution of transitions in emitter space to the coupling coefficients between $\tilde {\mu }={0}^{N-2}{1}^{2}$ and ${\tilde {\mu }}^{\prime }={0}^{N-3}{1}^{3}$. Notice that ${n}_{\tilde {\mu }}(0)=N-2$ and ${n}_{{\tilde {\mu }}^{\prime }}(1)=3$.

λ	${J}_{{0}^{N-3}{1}^{3},{0}^{N-2}{1}^{2}}^{(\boldsymbol{\lambda })}$	ρJ	${L}_{{0}^{N-3}{1}^{3},{0}^{N-2}{1}^{2}}^{(\boldsymbol{\lambda })}$
[N]	$\frac{N}{2}$	0	$\sqrt{3(N-2)}$
[N − 1, 1]	$\frac{N}{2}-1$	1	$\sqrt{2(N-3)}$
[N − 2, 2]	$\frac{N}{2}-2$	2	$\sqrt{N-4}$
[N − 3, 3]	$\frac{N}{2}-3$	3	0

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

$\begin{equation}{\hat{J}}_{+}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}\left\vert \tilde {\mu },\boldsymbol{\lambda };\mathbf{y},\mathbf{w}\right\rangle ={L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}{C}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}({\mathbf{w}}^{\prime },\mathbf{w})\left\vert {\tilde {\mu }}^{\prime },\boldsymbol{\lambda };\mathbf{y},{\mathbf{w}}^{\prime }\right\rangle ,\end{equation} \tag{ 4.2 }$

it is possible to write

$\begin{equation}{\hat{J}}_{+}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}=\sum\limits _{\boldsymbol{\lambda }{\trianglerighteq}{\boldsymbol{\mu }}^{\prime }}\sum\limits _{i=1}^{\mathrm{dim}({S}^{\boldsymbol{\lambda }})}\sum\limits _{{j}^{\prime }=1}^{{K}_{{\boldsymbol{\lambda }\boldsymbol{\mu }}^{\prime }}}\sum\limits _{j=1}^{{K}_{\boldsymbol{\lambda }\boldsymbol{\mu }}}{L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}{C}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}({\mathbf{w}}^{\prime },\mathbf{w})\left\vert {\tilde {\mu }}^{\prime },\boldsymbol{\lambda };{\mathbf{y}}_{i},{\mathbf{w}}_{{j}^{\prime }}^{\prime }\right\rangle \left\langle \tilde {\mu },\boldsymbol{\lambda };{\mathbf{y}}_{i},{\mathbf{w}}_{j}\right\vert ,\end{equation} \tag{ 4.3 }$

and define

$\begin{equation}{\hat{U}}_{\tilde {\mu }}=\sum\limits _{\boldsymbol{\lambda }{\trianglerighteq}\boldsymbol{\mu }}\sum\limits _{i=1}^{\mathrm{dim}({S}^{\boldsymbol{\lambda }})}\sum\limits _{j=1}^{{K}_{\boldsymbol{\lambda }\boldsymbol{\mu }}}\left\vert \tilde {\mu },\boldsymbol{\lambda };{\mathbf{y}}_{i},{\mathbf{w}}_{j}\right\rangle \left\langle \tilde {\mu },\boldsymbol{\lambda };{\mathbf{y}}_{i},{\mathbf{w}}_{j}\right\vert .\end{equation} \tag{ 4.4 }$

The latter implies that, if the raising operator is written in the basis of the PBEs, the coefficients ${L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}{C}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda })}({\mathbf{w}}^{\prime },\mathbf{w})$ are its singular values, while the unitary operators ${\hat{U}}_{{\tilde {\mu }}^{\prime }}$ and ${\hat{U}}_{\tilde {\mu }}$ correspond to the matrices of singular vectors, which provide with the changes of basis between PBEs and SABEs. In short, singular value decomposition of ${\hat{J}}_{+}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}$ 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 (3) that there are symmetry labels such that λ ⊵ μ ' yet $\boldsymbol{\lambda }{\ntrianglerighteq}\boldsymbol{\mu }$ is solved by recognizing the set of states $\left\{\left\vert {\tilde {\mu }}^{\prime },\lambda ;y,{t}^{\prime }\right\rangle :\boldsymbol{\lambda }{\ntrianglerighteq}\boldsymbol{\mu }\right\}$ as the null-space of the raising operator, i.e. ${L}_{{\tilde {\mu }}^{\prime },\tilde {\mu }}^{(\boldsymbol{\lambda }{\ntrianglerighteq}\boldsymbol{\mu })}=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, between $\tilde {\mu }={0}^{3}{1}^{1}$ and ${\tilde {\mu }}^{\prime }={0}^{2}{1}^{2}$. The subspaces and couplings are given by

$\begin{equation}\left(\begin{matrix}\hfill \left\vert {1}_{1}\right\rangle \hfill \\ \hfill \left\vert {1}_{2}\right\rangle \hfill \\ \hfill \left\vert {1}_{3}\right\rangle \hfill \\ \hfill \left\vert {1}_{4}\right\rangle \hfill \end{matrix}\right)\underset{{\hat{J}}_{+}^{(1,0)}}{\overset{{\hat{J}}_{-}^{(1,0)}}{{\rightleftharpoons}}}\left(\begin{matrix}\hfill \left\vert {1}_{1}{1}_{2}\right\rangle \hfill \\ \hfill \left\vert {1}_{1}{1}_{3}\right\rangle \hfill \\ \hfill \left\vert {1}_{1}{1}_{4}\right\rangle \hfill \\ \hfill \left\vert {1}_{2}{1}_{3}\right\rangle \hfill \\ \hfill \left\vert {1}_{2}{1}_{4}\right\rangle \hfill \\ \hfill \left\vert {1}_{3}{1}_{4}\right\rangle \hfill \end{matrix}\right).\end{equation} \tag{ 4.5 }$

In this basis, the raising operator is

$\begin{equation}{\mathbf{J}}_{+}^{({0}^{2}{1}^{2},{0}^{3}{1}^{1})}=\left(\begin{matrix}\hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill \end{matrix}\right).\end{equation} \tag{ 4.6 }$

Its singular value decomposition yields

$\begin{equation}{\mathbf{U}}_{{0}^{3}{1}^{1}}=\frac{1}{2}\left(\begin{matrix}\hfill 1\hfill & \hfill -\sqrt{3}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill \frac{1}{\sqrt{3}}\hfill & \hfill -\sqrt{\frac{8}{3}}\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill \frac{1}{\sqrt{3}}\hfill & \hfill \sqrt{\frac{2}{3}}\hfill & \hfill -\sqrt{2}\hfill \\ \hfill 1\hfill & \hfill \frac{1}{\sqrt{3}}\hfill & \hfill \sqrt{\frac{2}{3}}\hfill & \hfill \sqrt{2}\hfill \end{matrix}\right),\end{equation} \tag{ 4.7 }$

$\begin{equation}{\mathbf{U}}_{{0}^{2}{1}^{2}}=\frac{1}{\sqrt{6}}\left(\begin{matrix}\hfill 1\hfill & \hfill -1\hfill & \hfill -\sqrt{2}\hfill & \hfill 0\hfill & \hfill -\frac{1}{\sqrt{2}}\hfill & \hfill \sqrt{\frac{3}{2}}\hfill \\ \hfill 1\hfill & \hfill -1\hfill & \hfill \frac{1}{\sqrt{2}}\hfill & \hfill -\sqrt{\frac{3}{2}}\hfill & \hfill \sqrt{2}\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill -1\hfill & \hfill \frac{1}{\sqrt{2}}\hfill & \hfill \sqrt{\frac{3}{2}}\hfill & \hfill -\frac{1}{\sqrt{2}}\hfill & \hfill -\sqrt{\frac{3}{2}}\hfill \\ \hfill 1\hfill & \hfill 1\hfill & \hfill -\frac{1}{\sqrt{2}}\hfill & \hfill -\sqrt{\frac{3}{2}}\hfill & \hfill -\frac{1}{\sqrt{2}}\hfill & \hfill -\sqrt{\frac{3}{2}}\hfill \\ \hfill 1\hfill & \hfill 1\hfill & \hfill -\frac{1}{\sqrt{2}}\hfill & \hfill \sqrt{\frac{3}{2}}\hfill & \hfill \sqrt{2}\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 1\hfill & \hfill \sqrt{2}\hfill & \hfill 0\hfill & \hfill -\frac{1}{\sqrt{2}}\hfill & \hfill \sqrt{\frac{3}{2}}\hfill \end{matrix}\right),\end{equation} \tag{ 4.8 }$

$\begin{equation}{\mathbf{S}}_{{0}^{3}{1}^{1}}^{{0}^{2}{1}^{2}}=\left(\begin{matrix}\hfill \sqrt{6}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \sqrt{2}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \sqrt{2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \sqrt{2}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation} \tag{ 4.9 }$

where ${\mathbf{S}}_{\tilde {\mu }}^{{\tilde {\mu }}^{\prime }}={\mathbf{U}}_{{\tilde {\mu }}^{\prime }}^{{\dagger}}{\mathbf{J}}_{+}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}{\mathbf{U}}_{\tilde {\mu }}$ is the matrix whose diagonal lists the singular values of ${\mathbf{J}}_{+}^{({\tilde {\mu }}^{\prime },\tilde {\mu })}$. Notice that these values correspond with the expected $\sqrt{2(N-1)}$ and $\sqrt{N-2}$ with degeneracy degrees of 1 and N − 1 respectively; furthermore, the null-space has dimension N(N − 3)/2.