Group-theoretic treatment of strong light–matter coupling with an arbitrary number of excitations

The TC Hamiltonian describing a system of $N$ two-level systems (TLSs) coupled uniformly to a single electromagnetic mode can be written under the RWA as $H_{\mathrm{TC}}=H_{0}+H_{I}$, where

with $\hat{S}:=\sum_{n=1}^{N}\hat{\sigma}_{-}^{(n)}$ being the collective spin lowering operator for $N$ TLSs. Under the RWA, the Hamiltonian commutes with the total excitation number operator $\hat{X}:=\hat{S}^{\dagger}\hat{S}+\hat{a}^{\dagger}\hat{a}$, meaning that the total excitation number $X$ is conserved. This allows one to separate the Hilbert space into an infinite sum of finite-dimensional excitation manifolds according to the eigenvalue $X$, as

where in reality the occupied Hilbert space is limited by some $X_{\,\mathrm{max}}$ as the energy of a physical system must be finite. This corresponds to a block-diagonal form of the Hamiltonian, known to have corresponding Hilbert space dimensions

up to the maximum excitation number considered [7].

The system can be given more structure by taking note of the symmetries of the Hamiltonian. The Hilbert space can be decomposed into photonic (Fock space) and material (TLS) parts as $\mathcal{H}=F\otimes\left(\mathbb{C}^{2}\right)^{\otimes N}$. This allows for a natural action of the Lie algebra $\mathfrak{su}(2)$ on the individual TLS spaces, and of the symmetric group $S_{N}$ on the tensor power of TLSs. One can then decompose the system into irreps of the symmetric group. With uniform coupling constant $g$, the Hamiltonian is invariant under the action of $S_{N}$ permuting the TLSs, making the closed dynamics respect this decomposition. Since the actions of the two groups commute, one can apply the well-known Schur-Weyl duality from representation theory, resulting in these irreps combining with the ones of $SU(2)$ under joint labels $\lambda$, which correspond to Young diagrams of $S_{N}$ or equivalently weights (spins) $S$ of the Lie group $SU(2)$. The cooperation number or total spin $S$ is the usual eigenvalue of the Casimir operator $\hat{S}^{2}$ of $\mathfrak{su}(2)$ [25]. The “2” of $SU(2)$ limits the number of rows of the considered Young diagrams (Appendix A) so that we can write $\lambda=(N-m,m)$, where $m$ is the number of cells on the second row in English notation [37]. Because a Young diagram cannot have more cells on the second row than the first row, we find that $0\leq m\leq N/2$ (for odd $N$, $0\leq m\leq\left\lfloor N/2\right\rfloor$). Each excitation manifold then further decomposes into irreps as

where $\pi_{m}$ is an irrep of $S_{N}$, $W_{m}$ is an irrep of $SU(2)$, $m=\frac{N}{2}-S$ connects the Young diagrams and weights of $SU(2)$, and $d(N,m)$ is the multiplicity of the irrep $\mathcal{H}_{m}^{(X)}$, given by the dimension of the corresponding irrep of the symmetric group $\pi_{m}$ [37]. The block-diagonal TC Hamiltonian of each manifold can then be written as

where each submatrix $H_{m}^{(X)}$ is repeated $d(N,m)$ times. Due to the combinatorial nature of the symmetric group, the integer $m$ is further limited by the number of excited TLSs being permuted among ground-state TLSs. The maximum value of $m$ within an excitation manifold is therefore

The dimensions of the irreps $\mathcal{H}_{m}^{(X)}$ follow the integer dimension of the weight spaces of $SU(2)$, given by $X$ and $m$ as

where $N-2m+1=2S+1$, using the composite spin quantum number $S$. If $m>X$, $\dim{\mathcal{H}_{m}^{(X)}}:=0$. Since $N-2m+1=2S+1$ is the dimension of the $SU(2)$ irrep of weight $S$, it must be the maximal dimension of each $\mathcal{H}_{m}^{(X)}$. The irreps are then repeated $d(N,m)=\dim{\pi_{m}}$ times. The dimension of the $S_{N}$ irrep $\pi_{m}$ is famously given by the Hook length formula, a combinatorial formula calculated using “hook lengths” of a Young diagram of shape $\lambda$ (see Appendix A). For two rows, this formula can be written analytically as

Combining the last two equations and writing out the binomial coefficient, one gets the dimension formula

The excitation number $X$ limits the dimension from above to $X-m+1$, resulting in the minimum between the two values. This dimensional structure is presented in Table 1 for $N=9$. In addition to the above, we see the symmetric nature of the binomial coefficient in the total manifold dimension Eq. (5). This “locking” of the irrep dimensions will affect later results in Section 4.

Fig. 2 shows the dimensions given by Eq. (11) as functions of the symmetry index $m$ for $N=1000$ and $X\in\{150,300,400,500\}$. Panels (a–d) show the high-$m$ regions in linear scale, while the entire range of $m$ for $X=500$ is considered in (e), with the dimension plotted in logarithmic scale. We see that the irrep maximal in $m$ holds the highest multiplicity up to $X\approx N/3$, after which the dimensions approach a Poissonian-like distribution near the high-$m$ side of the distribution. We also see that the low-$m$ irreps become comparably insignificant in multiplicity, and that the multiplicity of the $m=X$ irrep approaches zero as $X\rightarrow\frac{N}{2}$.

The structure can be easily understood by comparing to the usual $X=1$ case. The representation $\mathcal{H}_{m}^{(1)}$ is two-dimensional with multiplicity $1$, and $\mathcal{H}_{1}^{(1)}$ is one-dimensional with multiplicity $N-1$. This is the well-known case of the upper polariton (UP), lower polariton (LP), and $N-1$ dark states [24]. The $X=2$ case has also been previously studied [44, 11]. As we will show in the next section, one can identify the subspace $V_{0}^{(X)}$ as multipolaritons [7] and $V_{X}^{(X)}$ as dark states ($X<N/2$). The representations in-between correspond to dark polaritons, recognized by their photonic content being between these two extremal regimes [12, 14]. Fig. 2 then shows that the dominant role of dark states is overtaken by dark polaritons and that the well-known problem of the large number of dark states becomes more subtle at higher excitation numbers. Thus, radiant processes from high-$m$ irreps might become statistically relevant [39, 28].

Let us present the diagonalization of the first two excitation manifolds to set what we seek to generalize. The first manifold Hamiltonian can be decomposed according to Eq. (7) into

leading to the UP, LP, and $N-1$ dark states

with the (first-order) Hopfield coefficients and eigenvalues satisfying

The dark states naturally carry the eigenvalue $E_{TLS}$ and we have denoted $g_{N}:=\sqrt{N}g$. For the rest of this paper, we focus on the resonant case of $E_{TLS}=E_{c}$. With this zero detuning, the expressions in Eqs. (17)–(19) simplify in a clear way. This shows us that the TLS-cavity energy difference acts as a small disturbance to a more balanced situation. We focus on the resonant case (unless otherwise specified), which allows for further analytic solvability, with this perturbative picture in mind.

Compared to the $X=1$ manifold, a new irrep becomes available for $X=2$, leading to the Hamiltonian decomposition [11]

with the values of $m$ corresponding to multipolaritons, dark polaritons, and dark states, respectively. These correspond to eigenstates of the form

where $\lvert E^{(X,m)}_{k,p}\rangle$ refers to the $k$th eigenstate (in increasing energy) within the $X$th manifold and $m$th irrep, with $p$ indexing equivalent states of given multiplicity. These states correspond to the eigenenergies and (second-order) Hopfield coefficients

accordingly for $i=0,1,2$. Note that the $m=1$ eigenstates have the first-order Hopfield coefficients Eqs. (17)–(18).

In what follows, we generalize the examples above. We find proper bases for all the irreps for given $N$ and $X$ to write the matrix elements of the Hamiltonian in symmetric form. We then use this form to derive a variety of properties for the systems considered by the model and to solve the eigensystem.

Take the unnormalized basis of symmetric states of the TC system [15]. Add irrep information carrying complex-valued symmetric coefficients $c_{\mathcal{I}}^{(\mathcal{J})}$, where $\mathcal{I}$ and $\mathcal{J}$ are index sets. These are determined by the action of the Young symmetrizer on the $N$-TLS states [37, 48, 18]. Let us denote arbitrary vector form states of $V_{m}^{(X)}$ in this basis as

where $\alpha_{i}$ is the amplitude of a symmetric state with $n$ photonic excitations, for example

By orthogonality of the irreps, one finds that the symmetric coefficients obey

where $j\in\mathcal{J},j^{\prime}\in\mathcal{J}^{\prime}$, and $\mathcal{I},\mathcal{J}$, and $\mathcal{J}^{\prime}$ are index sets of finite order, with $|\mathcal{J}|=m$ and $|\mathcal{I}|=X_{\,\mathrm{TLS}}$, the latter being the number of atomic excitations present. We further explicitly define $c_{\mathcal{I}}^{(\emptyset)}:=1$ for completeness. It should be noted that it is justified to take this basis into use, as these are not yet physical states. The key point is that we can use these as ansatz to find the true normalized bases through similarity transformations.

The interaction Hamiltonian only allows transitions between basis states of adjacent photonic content, and so Eqs. (1)–(2) expanded in the symmetric basis is a tridiagonal matrix. Under the resonance condition $E_{TLS}=E_{c}$, the free Hamiltonian is a scalar matrix $XE_{TLS}1\kern-2.5pt\text{l}$, and so the eigenvalues and eigenvectors of the system depend only on those of the rescaled interaction matrix

where $(H_{I})_{X-m+1}=gA$. If the original eigenvalue equation of the TC Hamiltonian reads $H_{\text{TC}}\lvert V^{(X)}_{m}\rangle=E\lvert V^{(X)}_{m}\rangle$, we have reduced the problem to solving the eigenvalues of $A$, $A\vec{\alpha}=\lambda\vec{\alpha}$, where $\lambda:=\frac{E-XE_{TLS}}{g}$.

Let us then find the matrix elements of $A$. Consider the two terms of the interaction Hamiltonian, $H_{I}^{[1]}:=\sum_{n=1}^{N}\lvert e_{n}\rangle\langle g_{n}\rvert\hat{a}$ and $H_{I}^{[2]}:=\sum_{n=1}^{N}\lvert g_{n}\rangle\langle e_{n}\rvert\hat{a}^{\dagger}$, acting on bare states of the type $\lvert V\rangle=\lvert e^{k}\rangle\lvert l\rangle$, where $\lvert e^{k}\rangle=\lvert e_{n_{1}}e_{n_{2}}\cdots e_{n_{k}}\rangle$ , and $k+l=X$. The cavity operators $\hat{a}\penalty\ (\hat{a}^{\dagger})$ operating on the bare states simply correspond to factors of $\sqrt{l}\left(\sqrt{l+1}\right)$. For the $H_{I}^{[1]}$ input state $\lvert e^{k-1}\rangle$ and output state $\lvert e^{k}\rangle$, $\binom{k}{k-1}=k$ states can transform into $\lvert e^{k}\rangle$, so $\sum_{n=1}^{N}\lvert e_{n}\rangle\langle g_{n}\rvert$ leaves a factor of $k$. Similarly for $H_{I}^{[2]}$ input state $\lvert e^{k+1}\rangle$ and output state $\lvert e^{k}\rangle$, $N-k$ of states $\lvert e^{k+1}\rangle$ can transform into $\lvert e^{k}\rangle$, leaving a factor of $N-k$. The symmetric coefficients change according to the corresponding irrep. We summarize these results in Table 2.

The procedure we want to generalize is as follows: Find the matrix elements in the symmetric basis and transform them into the properly normalized basis, for general $N$, $X$, and $m$. If we apply the results of Table 2 on the symmetric basis and note that the off-diagonal sequences contain $X-m$ terms, we find that the matrix elements are $a_{l}^{(C)}\cdot a_{l}^{\downarrow}$ and $a_{l}^{(C)}\cdot a_{l}^{\uparrow}$ for the sub- and superdiagonal sequences with values given in Table 3.

The Hamiltonian is diagonalizable and Hermitian and therefore symmetric and real in some basis [$a_{l}=a_{l}^{\prime}\in\mathbb{R}$ in Eq. (3.1)]. This coincides with the properly normalized basis, and the proper matrix $A$ can be found recursively: The matrix $P$ carrying the basis transformation $PAP^{-1}$ is diagonal, and the first element is $P_{00}=1$. The second element is $P_{11}=\frac{\sqrt{X}}{\sqrt{N-X+1}}$, comprising of the square roots of the first terms in the sub- and superdiagonal sequences. The rest are then determined recursively as $P_{ll}=P_{l-1,l-1}\cdot\frac{\sqrt{c^{\uparrow}_{l}}}{\sqrt{c^{\downarrow}_{l}}}$. The new matrix element coefficients are of the form $c_{l}:=c_{l}^{(N)}\cdot c_{l}^{(C)}\cdot c_{l}^{(M)}$, where the cofactors $c_{l}^{(i)}$ represent normalization ($N$), cavity ($C$), and TLS (matter, $M$) contributions. The matrix elements $(A)_{ij}=c_{l}\delta^{|i-j|}_{1}\delta^{\min{\{i,j\}}}_{l}$ are presented in Table 4.

As a real symmetric tridiagonal matrix, $A$ has many useful properties. An immediate consequence is that all the eigenvalues are simple and real. This means that the representation-theoretic multiplicity then determines also the multiplicities of the eigenvalues. If $D:=\,\mathrm{diag}(\dots,-1,1,-1,1)$, it is easy to see that $A$ and $D$ anticommute, so the spectrum of $A$ is symmetric about zero (see Appendix B.1). A direct corollary is that $0$ is an eigenvalue for $X-m\equiv 0\mod{2}$, and that one has to solve only one half of the spectrum, allowing for significant numerical optimization. As $m$ restricts only the length of the matrix element sequence above—within a given excitation manifold—the matrix corresponding to $m^{\prime}=m-1$ is a principal submatrix of that of $m$. It follows from the Cauchy interlacing theorem that the eigenvalues of $H^{(X)}_{m}$ and $H^{(X)}_{m-1}$ alternate [38]. Starting recursively from $H^{(X)}_{X}$, it follows that $H^{(X)}_{0}$ has the extremal eigenvalues of an excitation manifold. These extremal values have a bound proportional to the matrix elements according to the Geršgorin disc theorem [21]. By estimating the values of Table 4 upward for $m=0$ by $\sqrt{l\cdot(X-l)}\leq X/2$ and $c^{(N)}_{l}\leq\sqrt{N}$, a Geršgorin bound of $X\sqrt{N}$ is found. Later calculations show that this is very close to the true extremal values.

With the bases and matrix elements derived above, we are able to calculate the eigenvalues $E^{(X,m)}_{k}=XE_{TLS}+\lambda^{(X,m)}_{k}g$, the eigenvectors $\lvert E^{(X,m)}_{k,p}\rangle$, their probability amplitudes squared $|\alpha_{i}|^{2}$, average photonic content $\langle n_{c}\rangle$, and average TLS excitation content $\langle n_{\,\mathrm{TLS}}\rangle$. These quantities are presented in Table 5 for $X=1,2,3$, where we omit writing the indices $(X,m)$ where it does not cause confusion. The vector components $\alpha_{i}$ have convergent $N$-dependence, strong already at realistic values of $N>10^{6}$, so we have calculated $\lim_{N\rightarrow\infty}\alpha_{i}$ for the tabled values (more on the limit in Section 4). We will continue by inspecting each quantity at a time.

We see from Table 5 that the eigenvalues tend to follow a pattern of nested square roots containing polynomials of $N$. While it may be possible, an analytical expression for general $N$, $X$, and $m$ seems to be highly non-trivial to find. The properties of $A$ described above allow for heavy optimization with numerical methods. For the following calculations, we used the NumPy library for Python. Fig. 3 shows the eigenvalues $\lambda$ of $A$ against their multiplicities for $N=1000$ and $X\in\{3,10,75,350,450,500\}$. The color encodes the symmetry index $m$. Panels (a–c) are in logarithmic scale and $0\leq m\leq X$. Panels (d–f) are in linear scale and $0.8X\leq m\leq X$. Recall that $m=0$ corresponds to multipolaritons and $m=X$ to dark states, with maximal and minimal $\langle n_{c}\rangle$, respectively.

We see that the statistical relevance of differing energies centers near the dark states. This effect is countered by the fact that lower-$m$ irreps are more dense in energy states due to higher dimensions and the spectrum being bounded. The eigenvalues are mirror-symmetric around zero. We also recognize the effects on the relevance of dark states seen already in Fig. 2. Furthermore, photonic energy states appear at the dark-state energy for even $X-m$ due to the symmetric spectrum, although this is perturbed by the detuning, i.e., $E_{TLS}\not=E_{c}$ [44].

Note that the multiplicities grow considerably between panels. This growth is stunted after $X>\frac{N}{2}$ as no new irreps appear and the energy levels only get denser (Table 1). This can be understood as approaching a continuous limit for the energy spectrum. As the $S_{N}$ irreps define the TLS excitation content of the eigenstates, dark states must disappear after this limit, with new excitations going to the cavity.