Integrability, multifractality, and two-photon dynamics in disordered Tavis-Cummings models

The Tavis-Cummings model is a paradigmatic central-mode model where a set of two-level quantum emitters (spins) are coupled to a collective cavity mode. Here we study the eigenstate spectrum, its localization properties and the effect on dynamics, focusing on the two-excitation sector relevant for nonlinear photonics. These models admit two sources of disorder: in the coupling between the spins and the cavity and in the energy shifts of the individual spins. While this model was known to be exactly solvable in the limit of a homogeneous coupling and inhomogeneous energy shifts, we here establish the solvability in the opposite limit of a homogeneous energy shift and inhomogeneous coupling, presenting the exact solution and corresponding conserved quantities. We identify three different classes of eigenstates, exhibiting different degrees of multifractality and semilocalization closely tied to the integrable points, and study their stability to perturbations away from these solvable points. The dynamics of the cavity occupation number away from equilibrium, exhibiting boson bunching and a two-photon blockade, is explicitly related to the localization properties of the eigenstates and illustrates how these models support a collective spin description despite the presence of disorder.


I. INTRODUCTION
In the past years cavity quantum electrodynamics (QED) has emerged as a new platform for the quantum simulation of light-matter interactions [1].In such cavity setups the electromagnetic cavity mode is coherently coupled to a set of quantum emitters behaving as pseudospins, allowing for polaritons arising as a strong hybridization between the photonic mode of the electromagnetic field and the emitter modes.Furthermore, the cavity-mediated interactions in these models are highly tunable, allowing for systematic studies of e.g. the role of disorder in a controlled way [2,3].Motivated by these experimental developments, there have been a wealth of studies on the interplay between strong light-matter coupling and disorder, specifically on the consequences on the localization properties of the eigenstates [4][5][6][7].Cavity models can exhibit 'bright' polaritonic eigenstates with strong light-matter hybridization and 'dark' states with minimal entanglement between the cavity and the emitters, and different classes of states can have different localization propreties that respond differently to the presence of disorder, which is in turn reflected in the dynamical properties of these models [8].Ref. [4] argued that dark states exhibit 'semilocalization', being localized on multiple noncontiguous sites.Bright and dark states were also shown to exhibit multifractality, being neither fully localized nor delocalized, in a single-excitation setup described by arrowhead matrices [5].
However, the presence of disorder is not the only factor that needs to be taken into account: the many-body dynamics and eigenstate properties can be modified by integrability.Various cavity models are either integrable or close to an integrable point, where the model can be exactly solved using Bethe ansatz and the dynamics is constrained by an extensive amount of conservation laws.Variants of the Tavis-Cummings model, the main focus of this work, were originally shown to be integrable by Gaudin [9] and later a set of conserved charges were identified by Dukelsky et al. [10].Such cavity models exhibiting 'one-to-all' interactions are part of a family of integrable Richardson-Gaudin Hamiltonians [11][12][13].Crucially, these models remain integrable in the presence of disorder, and their eigenstates have been shown to exhibit anomalous localization properties [14].Indeed, the effect of integrability was subsequently studied in Ref. [7], where it was shown that integrability-breaking restores the ergodicity of the eigenstates in the thermodynamic limit for a finite density of excitations.
One paradigmatic model in cavity QED is the Tavis-Cummings model [15,16], describing the interaction of  pseudospins with a central bosonic mode under the rotatingwave approximation [17].Its Hamiltonian can be written as where â † ( â) are the bosonic creation (annihilation) operators for the cavity mode and Ŝ  spin operators describing pseudospin .We will assume that these are spin-1/2 pseudospins, representing two-level quantum emitters, but this is not a necessary assumption.This model can be realized in both cavity and circuit QED [3,8,18,19].The first term gives the energy of the bosonic mode and the second term that of the quantum emitters (each with bare Zeeman energy   ).The final term describes the inhomogeneous coupling between the bosonic mode and the pseudospins with interactions strengths   .Here we have allowed for disorder in both the bare energies and the couplings between the cavity mode and the atoms.The Hamiltonian commutes with the total excitation number, , such that we can restrict ourselves to sectors with a fixed number of excitations.For any sector with a finite number of excitations, i.e. the number of excitations does not scale with , the factor 1/

√
in the interaction term is required to have an finite energy in the thermodynamic limit of infinite system size  → ∞.In this way this factor plays the same role as the Kac factor in systems with long-range interac-tions [20], an analogy that will be elaborated in the remainder of this work.
Theoretical studies of the model in the one excitation sector,  = 1, have already been carried out in Ref. [4], illustrating the multi-fractal structure of the eigenstates along with their semilocalization properties.The dynamics in this sector was also studied in Ref. [21], where the presence of disorder gave rise to a variety of complex behaviors.However, nonlinear photonics require extending the Hilbert space to include multiple excitation [22].Here, we work within the two excitation sector,  = 2, to study the localization properties of the eigenstates and their effect on the dynamics.While the model has long been known to be integrable for homogeneous coupling strengths and inhomogeneous bare energies, we show that the opposite limit is also integrable.This exact solution is remarkably transparent and allows us to identify three different classes of eigenstates.These two integrable limits conspire to result in multifractal eigenstates where the localization properties are surprisingly robust to integrability-breaking perturbations.
The eigenstates properties are directly reflected in the dynamics, and we use the nearby integrable limit to present exact results for the short-time dynamics of disordered Tavis-Cummings models.We note that the structure of these cavity models is closely related to central spin models, where a set of noninteracting spins interact with a central spin, similar to the cavity mode in these setups.Such central spin models similarly support both bright and dark states [23][24][25][26][27][28][29][30][31], now depending on the hybridization between the central spin and the environment, and recent studies have shown that these also exhibit the combination of Richardson-Gaudin integrability and semilocalization and multifractality [32].
This work is structured as follows.In Sec.II we consider the Tavis-Cummings Hamiltonian for two excitations in the absence of disorder, highlighting the existence of three classes of eigenstates.In Sec.III we present the exact solution for the model in the presence of disorder in the interaction strenghts, reducing the diagonalization of the Hamiltonian to solving a single nonlinear equation with different classes of solutions, and present the exact conserved charges.For interaction strengths scaling as in Eq. ( 1), we show how the dynamics are further constrained due to the approximate conservation of permutation operators and derive the corresponding relaxation times in Sec.IV.For completeness, the exact solution of the model for inhomogeneous bare energies is presented in Sec.V.The localization properties of the three different classes of eigenstates are discussed in Sec.VI, and these are related to different examples of the dynamics of the cavity mode occupation number in Sec.VII.Sec.VIII presents our conclusions.

II. THE HOMOGENEOUS TAVIS-CUMMINGS HAMILTONIAN
In order to understand the spectrum of the two-excitation Tavis-Cummings Hamiltonian and distinguish different classes of eigenstates, it is instructive to first consider the homogeneous limit.In the absence of disorder, when all the bare en- where we have set   = , ∀ and   = , ∀.The Hamiltonian can now be expressed by introducing total spin states | tot ,   tot ⟩ = |,   ⟩, with total spin  and total spin projection   .
The restriction to the two-excitation sector requires that there can only be up to two spin excitations, such that   can only take the values −/2, −/2 + 1 and /2 + 2 (in what follows, we assume that  ≥ 4).The Hamiltonian additionally conserves total spin quantum number, which can take a maximal value of  = /2, such that the Hamiltonian decomposes in three blocks with different total spin.Either  = /2 − 2, for which there is a single possible state or  = /2 − 1, resulting in two states or  = /2, resulting in three states In this basis of multiplet states the Hamiltonian takes the form Here we have taken the limit of  ≫ 1 in all matrix elements and neglected subleading corrections in 1/, which however do not change the structure of this matrix.This Hamiltonian clearly returns three classes of eigenstates, corresponding to the three different blocks: the singlet states from Eq. ( 3) return eigenstates |0⟩⊗|/2 − 2, −/2 + 2⟩ where the spin degrees of freedom do not hybridize with the cavity mode.These states act as dark states.The doublet states from Eq. ( 4) hybridize states with no photon in the cavity and with a single photon in the cavity.We will refer to these states as doublet polaritons.Finally, the triplet states (5) result in states that are a linear combination of zero, one, or two photonic excitations, such that we will refer to these as triplet polaritons.While in the literature the term polaritons is usually reserved for the strong-coupling regime, we will always refer to these states as polaritons for clarity (also following Ref.[5]).
The degeneracy of these eigenvalues is given by the total number of ways in which  spin-1/2 particles can be coupled to total spin .This number of ways follows from Catalan's triangle, resulting in a total degeneracy of  ( − 3)/2 for the singlet state, a total degeneracy of ( − 1) for each of the doublet states, and the triplet states are nondegenerate since there is only a single way of coupling  spin-1/2 particles to total spin  = /2.

III. DISORDERED COUPLINGS
The eigenspectrum of the Tavis-Cummings model is illustrated in Fig. 1 for varying Δ, both in the absence and presence of disorder.In the presence of disorder the highly degenerate dark states and the doublet polaritons split into bands of states, with the number of states corresponding to the degeneracy in the absence of disorder, whereas the triplet polaritons remain isolated states.For sufficiently weak disorder strengths the different bands do not overlap, which we will take to be the case in the remainder of this work.
In the following we will show that the presence of weak disorder in the coupling between the cavity mode and the pseudospins indeed preserves the three classes of eigenstates.In this limit the Hamiltonian can be written as with For convenience we have defined   =   / √ .Since the Hamiltonian commutes with the total excitation number , we can set  to zero without loss of generality.We first present an explicit construction of the eigenstates, proving the exact solvability of the Tavis-Cummings model (in the two-excitation sector) for disordered couplings and in the absence of disorder in the bare energies.This exact solution then also allows us to construct an extensive set of conserved charges commuting with the Hamiltonian.The integrability of this model was conjectured in Ref. [32], and we here prove a limited version of this conjecture, showing that the integrability holds if the model is restricted to the two-excitation sector.

Dark states
In the presence of a disordered coupling the model still allows exact dark states.These are eigenstates of the Hamiltonian (7) of the form |0⟩ ⊗ |D⟩, where the spin state satisfies These states are adiabatically connected to the singlet states from Eq. (3), since in the limit of a homogeneous interaction strength Ĝ− ∝ Ŝ− tot , and are well studied in the literature on central mode models [23][24][25][26][27][28][29][30][31].The total number of dark states is  ( − 3)/2, and these states are annihilated by the interaction part in the Hamiltonian, such that these are again (highly degenerate) eigenstates of the Hamiltonian (7) with eigenvalue 2, since by construction

Polaritons
For the polaritons, the exact diagonalization of the Hamiltonian (7) can be reduced to solving a nonlinear equation for the eigenvalues.In the spirit of integrability, we can consider an ansatz for the polariton eigenstates with a single free parameter, in such a way that the eigenvalue equation reduces to a nonlinear equation for this parameter.
Specifically, we consider an ansatz for the (unnormalized) eigenstates of the form where  is a free variable and |  ⟩ is a wave function for the  spin degrees of freedom with a single spin excitation, dependent on .Acting with the Hamiltonian (7) on the ansatz (11), it is straightforward to show that this state is an eigenstate with eigenvalue  + Δ provided This equation is a self-consistent eigenvalue equation since the spin Hamiltonian Ĥ () explicitly depends on the eigenvalue .Note also that for a normalized state |  ⟩ the norm of the polariton wave functions (11) can be calculated from the Hellmann-Feynman theorem as Crucially, the Hamiltonian in Eq. ( 12) is integrable for every choice of  and its eigenvalues and eigenstates can be explicitly constructed.As shown in Appendix IX, using this exact solution the self-consistent eigenvalue equation can be recast as a secular equation.Defining ḡ2 =  =1  2  , Δ = Δ/2, a state |  ⟩ can be written as with |∅⟩ = |↓, . . ., ↓⟩.The eigenvalue equation for the polariton energy  then reduces to solving the following equation for Remarkably, different solutions to this equation can be directly identified with either the doublet or triplet polaritons.In Fig. 2 the structure of this equation is made explicit by plotting both sides as a function of , and the intersections between these curves correspond to the solutions of Eq. ( 15).
The left-hand side has a vertical asymptote at  = Δ/3 and a horizontal asymptote at 1/3, and is monotonically increasing everywhere.The right-hand side is an even function of  and has a series of poles at with each pair of poles corresponding to a coupling strength   .The function is monotonically decreasing (increasing) for  positive (negative) and goes to zero for || → ∞.There are two classes of solutions to these equations: the solutions corresponding to the doublet states lie in between two poles, leading to two sets of ( − 1) states.The three triplet states correspond to the remaining isolated solutions to these equations away from the poles.Taking these together, we find the two curves generally have 2 + 1 intersections, exhausting the available polariton Hilbert space.Note that for a sufficiently asymmetric distribution of interaction strenghts a pair of poles can vanish, but the secular equation will still have 2 + 1 solutions, as discussed in Appendix XI.We emphasize that, while diagonalizing the exact Hamiltonian rapidly becomes unfeasible for large system sizes, the secular equation ( 15) can be efficiently solved for arbitrarily large system sizes using e.g. an intersection method, where specific triplet or doublet states can be systematically targeted since their relative position w.r.t. the poles is known.
In the limit where  is far away from the poles, the secular equation approximately reduces to a cubic equation with three solutions and corresponding energies  =  + Δ given by: These eigenvalues now reduce to the eigenvalues for the triplet states from Eq. ( 6).We note that the secular equation ( 15) closely resembles the dispersion equation within the random phase approximation (RPA), which aims to construct approximate particle-hole excitations on top of a reference state [33].This similarity arises more generally within the theory of Richardson-Gaudin integrability (see e.g.Ref. [34]): for single-excitations states the Bethe equations typically reduce to the dispersion relation for the Tamm-Dancoff approximation (TDA), which aims to construct particle-like excitations on top of a reference state, and the RPA can be seen as the two-excitation generalization of the TDA [33].

Conservation laws
The conserved charges for the Hamiltonian ( 7) can be constructed from those of the factorizable Richardson-Gaudin Hamiltonians [9,11,13].These conserved charges are easiest to represent in a block matrix representation, similar to Eq. ( 6), and are of the form where Q  are operators acting on the spin degrees of freedom as That they commute with the Hamiltonian (7) for two excitations can be verified by direct calculation and is shown in Appendix X.The existence of these conserved quantities guarantees the integrability of the model.We note that this calculation is highly similar to a recent calculation of the conserved charges in a spin-1 central spin model [32], a related central mode model where a similar block matrix structure appears.
For this reason, we will defer from discussings these conserved charges in more detail and refer the reader to Ref. [32].We only note that this previous work conjectured the integrability of the Tavis-Cummings model with disordered couplings (for an arbitrary number of excitations), and this work establishes its integrability in the limiting case of two excitations.The integrability for a higher number of excitations remains an open question.

IV. RELAXATION TIMES FOR PERMUTATION SYMMETRY
The Bethe states in the previous derivation can be adiabatically connected to collective spin states.This correspondence suggests that, at least for not too strong disorder, the considered model can be described in terms of the collective spin operators of the homogeneous model.Here we show that for the Hamiltonian (1) such a description of the dynamics is justified up to times scaling as  ∝ √  in the absence of disorder in the bare energies.More specifically, collective spin states are indicative of an underlying spin permutation symmetry, and we show that this permutation symmetry is preserved up to relaxation times scaling as

√
in the new integrable model.This derivation builds on a similar argument for systems with sufficiently long-range interactions [35].Such models similarly support a description in terms of collective spin states in the presence of a nearby integrable (but quasiclassical) limit [35][36][37][38].These relaxation times then present an additional similarity between the current model and lattice systems with sufficiently long-range interactions.In such systems the presence of an integrable semiclassical limit was recently argued to be crucial [38] for such a description, and we here show the stability of the collective spin description near the integrable point.
In the homogeneous models all spin modes are identical, i.e. the Hamiltonian is invariant under any permutation of the spins.
The Hamiltonian commutes with permutation operators P  , permuting spins  and , and e.g. for spin-1/2 we can write these permutations in terms of Pauli matrices as These operators are exactly conserved in the dynamics of the homogeneous model, for which the eigenstates in terms of collective spins are similarly eigenstates of the permutation operators.While the permutations operators are no longer exactly conserved in the inhomogeneous case, it is possible to derive the inequality for systems with homogeneous bare energies, indicating that permutation symmetry is conserved up to a time scale As such, for any initial state that satisfied this permutation symmetry, the dynamics up until the time scale    can be accurately modeled using the collective spin operators, since these are exactly the operators that do not take the states out of the initial symmetry sector.In the case of homogeneous bare energies, |  −   | scales (at worst) as 1/

√
, resulting in time scales    ∝ √ , indicating that any such description becomes increasingly accurate with increasing system size.Note that this is a lower bound -the scaling as  −1/2 holds for the extremal values of   and   , whereas the typical closest distance scales as  −3/2 , indicating much longer-lived conservations for the corresponding spin permutation operator.
The proof is straightforward and follows a similar proof for systems with long-range interactions from Ref. [35].We have that where we have used that where we have bounded the expectation value in terms of the operator norm and used that || P  || = 1.The operator P  Ĥ P  is the original Hamiltonian with spins  and  exchanged, such that we find While the operator norm of â and â † is in general unbounded, here we can make use of the restriction to the two-excitation subspace, for which The operator norm of the above operator norm is bounded by the sum of the operator norms of each term, resulting in Plugging this bound in Eq. ( 24), we obtain the proposed bound from Eq. ( 22).While the integrability of the model guarantees nonergodicity due to the presence of conservation laws, this result further constrains the dynamics.For any initial state that is fully symmetric in the spins, as is e.g. the case for a state with no spin excitations and 2 photonic excitations, the dynamics can be restricted to the Dicke manifold, i.e. spin states that are fully symmetric under spin exchange, up to a time scale scaling as √ .This space is also known as the totally symmetric subspace (TSS), and the role of the TSS in spin dynamics has been the subject of active study [39].Only after this relaxation time scale can the system move out of the TSS, and we find that the TSS is increasingly stable for increasing system sizes.
Note also that the above relaxation time only depends on the restriction to the two-excitation sector by setting For the -excitation sector, in the above derivation the factor √ 2 only needs to be replaced by √ .In any sector where the number of excitations does not scale with system size, we hence expect relaxation times for the permutation symmetry scaling as If the number of excitations  scales with , i.e.  ∝ , then the Kac factor in the Hamiltonian ( 7) also needs to be modified to 1/, i.e.
Repeating the derivation above directly results in the bound again indicating that permutation symmetry is preserved up to time scales scaling as √ .We obtain the general result that permutation symmetry is preserved up to time scales scaling as √ , irrespective of the number of the excitations, provided that the Hamiltonian is defined in such a way that energy is finite.This argument fails for a disorder in the bare energies   , since then all obtained time scales would be  (1).In this limit, however, it appears that integrability again stabilizes the collective spin description.
In the context of long-range systems, we note that the scaling with system size of the corresponding relaxation times follows directly from the Kac factor fixing the extensivity of the energy, similar to how the time scale in our context requires the correct normalization of the interaction strengths.

V. DISORDERED BARE ENERGIES
For completeness, we reiterate the exact solution of the Hamiltonian with homogeneous couplings and disorder in the bare energies, It is known that in this limit the model is integrable [10] and can be solved using the Bethe ansatz [40] (see also Refs.[41][42][43][44]).
The eigenstates can then be written as Bethe states, characterised by two variables  1 and  2 : expressed in terms of generalized raising operators acting on the vacuum state |∅⟩ = |0⟩ ⊗ |↓ . . .↓⟩.These states are eigenstates of the partially homogeneous Hamiltonian with total energy  =  1 +  2 , provided the variables satisfy the set of Bethe equations [41]: The two variables  1,2 are also referred to as 'quasi-energies' due to their role in the energy of the eigenstate they describe.These equations have been well studied in the literature [40,[45][46][47].For our results it is relevant that there exist different classes of solutions, depending on the position of the variables  1,2 with respect to the poles   in the Bethe equations: either both  1 and  2 are both far away from the poles, or one variable is away from the poles and the other is 'trapped' between a pair of poles, or both variables 'trapped' between pairs of poles.Following our discussion for the limit of inhomogeneous couplings, these solutions can be identified with triplet polaritons, doublet polaritons, and singlet dark states respectively.
The Hamiltonian again supports an extensive set of conserved quantities, one for each spin in the system, where now These form a set of mutually commuting conserved charges satisfying [ Ĥ, Q ] = [ Q , Q  ] = 0, ∀, .

VI. SEMILOCALIZATION AND MULTIFRACTALITY
The different classes of eigenstates do not just differ in that they belong to different bands, but they also have different localization properties.Localization for eigenstates in cavity models have recently gained attention because of their anomalous localization properties, which in turn directly translate to a lack of thermalization for a local perturbation [7].For a single excitation dark states eigenstates were argued to be 'semilocalized', i.e. being localized on multiple noncontiguous sites, in Ref. [4].In a follow-up work, it was argued that the polaritons in such a model exhibit multifractality: the eigenstates are extended, i.e. not localized, but yet non-ergodic and not fully delocalized [5].Multifractality was similarly observed in the integrable Tavis-Cummings model with a finite excitation density [7], and arguably already appeared in earlier studies of the (closely related) integrable Richardson model [14].This multifractality was similarly observed in Ref. [32] for an integrable central spin model where the central mode is a spin-1 particle, where the dynamics of the central spin mode can serve as a probe for the multifractality.
In order to probe the multifractality of a state |⟩, we consider the -dependent inverse participation ratio, defined as: where  are product states in the full (photonic and spin) dimensional Hilbert space.The -dependent IPR quantifies the distribution of the components of an eigenstate in a product state basis, with  acting as the equivalent of the order in the Rényi entropies.For a delocalized eigenstate all coefficients scale as 1/ √ , resulting in an IPR scaling with dimension of the Hilbert space as  1− .A change in this scaling as  is varied is a signature of multifractality in the eigenstate [48,49]. The

Ŝ+
|∅⟩ and introducing a corresponding notation for the eigenstate components.The IPR then reads In the following, we show that it is useful to consider the different contributions to the IPR separately and introduce 'restricted' versions of the IPR, where the summations are restricted to the sectors with a fixed number of photons.Specifically, we write The scaling of the full IPR will be determined by the scaling of the largest term of these three.
In Fig. 3 we present numerical result for these three components as a function of system size  in the presence of disorder in the couplings as well as (weak) disorder in the bare energies.The disorder in the bare energies is chosen to be sufficiently weak such that the different bands from Fig. 1 do not mix, allowing different states to be identified based on their position in the spectrum.We consider a single disorder realization, since the fluctuations over the different eigenstates within the bands are smoothened by averaging the IPR within each class of eigenstates.In this limit we find that all three classes of states exhibit multifractality, with different scaling exponents.These scaling exponents are numerically observed to be identical to the exponents obtained in the integrable limit with a homogeneous coupling in all sectors.For inhomogeneous couplings and homogeneous bare energies, care needs to be taken.First of all, the dark states are exactly degenerate, such that it is not meaningful to consider the localization properties of a single dark states.Second, we observe that the scaling exponents in the 0and 1photon sector are identical to the exponents in the presence of both sources of disorder, but the single amplitude of the 2-photon component can have different scaling.However, this amplitude does not contribute to the full IPR, such that the scaling of the total IPR will be identical in both the two integrable limits and in the non-integrable case.The scaling of the IPR for the Bethe states is derived in Appendix XII, and we here focus on the scaling away from these integrable limits.
Triplet polaritons.First focusing on the three triplet polariton states, we observe the following the scaling in the large- limit: independent of .However, the different scaling of these contributions indicate that the total IPR will exhibit a change in scaling as  in varied, with either the 0-photon term being dominant ( < 1) or the 2-photon term being dominant ( > 1).
The IPR for the triplet polaritons will then scale as ,  > 1. ( This result has a direct interpretation: the total weight of the triplet states within each photon sector is  (1), i.e. the probability of observing  photons in these eigenstates is  (1) for all values of , as quantified in the restricted IPR's for  = 1: but within each sector these states are fully delocalized.Recall that a delocalized eigenstate in a Hilbert space of dimension  results in an IPR scaling as  1− for all values of .The 1-photon states span a Hilbert space of dimension , and the delocalization in this Hilbert space leads to a scaling  (1−) in Eq. ( 41).The 0-photon states span a Hilbert space of dimension  ( − 1)/2 =  ( 2 ) and the delocalization in this space leads to the scaling  2(1−) in Eq. ( 40).The single 2-photon state has a weight  (1), leading to the observed scaling from Eq. ( 42).While the states are delocalized within the separate photon sectors, they are not delocalized within the full Hilbert space due to the total  (1) weight of the states within each sector.For full delocalization within the Hilbert space the states within each photon sector would not be normalized [up to an  (1) factor] and the restricted IPRs would include additional scaling factors.Delocalization in the full Hilbert space would e.g.predict that the total weight | 0 | 2 of the 2-photon sector vanishes as the relative dimension of this sector, i.e.  ( −2 ), to be contrasted with the observed  (1) scaling.As such, it is the relative weights of the different sectors that give rise to the change in scaling of the IPR as  is varied, indicating multifractality.Doublet polaritons.For the doublet polaritons, we find that Remarkably, even within the subsections of the Hilbert space with a fixed number of photons, the eigenstates exhibit multifractality and are not fully delocalized.These scalings reflect underlying semilocalization in the 1-photon sector, i.e. there is an  (1) number of components dominating the scaling for  > 1/2, whereas in the 0-photon sector an  () number of states dominate.Interestingly, in the 0-photon sector the eigenstates are hence localized within a vanishing fraction of the Hilbert space [ () components in a  ( 2 ) Hilbert space].While the states in both the 0and 1-photon sector are distributed over a vanishing fraction over the Hilbert space, within the 1-photon sector these are localized on non-contiguous sites, whereas in the 0-photon sector these are distributed over noncontiguous regions in the Hilbert space (see App. XI).
The IPR for the doublet polaritons reflects the scaling of the dominant components in the restricted IPR, changing from the scaling of the 0-photon sector to the scaling of the 1-photon sector as Singlet dark states.The dark states similarly exhibit multifractality within each of the three different sectors, with the restricted IPR scaling as Taking these results together, the scaling of the IPR follows the 0-photon IPR, as could be expected for dark states: These scalings now reflect semilocalization in both the 0 and 1photon sector, where in both cases  (1) components dominate the IPR scaling for  large enough.
All presented scalings can be clearly numerically observed in Fig. 3, where the numerically obtained scaling exponents are close to the theoretically predicted values.However, we emphasize that these results are limited to weak disorder and system sizes up to  = 200.While it is possible that these scalings break down for larger system sizes, there is no indicator of this happening in our numerics.
These scalings are analytically derived for the integrable limit in Appendix XII.In the same way that the different classes of eigenstates can be connected to the relative position of the Bethe root to the poles in the secular equation ( 15), these different scalings can be directly related to the Bethe states.This connection is detailed in Appendix XII.Any Bethe root that lies close to a pole will lead to large contributions in the eigenstate component related to this pole and corresponding semilocalization on the corresponding basis states.Different localization properties are hence expected for states characterized by a different relative position of the Bethe root to the poles.These different scalings and the connection with the pole structure can be visualized by considering the eigenstate components within e.g. the single-spin excitation basis states.If the set of bare energies is ordered, these components will be a smooth functions of   , with a possible divergence at  1 and  2 if these lie in between two poles.This behavior is illustrated in Fig. 4, and it is clear that the presence of (weak) integrability-breaking terms does not change the multifractal character of the eigenstates.For the triplet states all components exhibit the same scaling, where either one or two peaks appear for the doublet and singlet states respectively.

VII. DYNAMICS AND PHOTON BUNCHING
These localization properties directly translate to the dynamical behavior of initial states with a fixed photon number.While it is customary to introduce leakage and describe the dynamics in terms of open systems, we here focus on closed system dynamics in order to keep the connection with the previous results.In this sense these dynamics is expected to be reflective of the short-time dynamics of realistic cavity models with dissipation, i.e. the dynamics on times shorter than the dissipation time scale.
Since the restricted IPR indicated different eigenstates localization properties depending on the number of photons, we consider the dynamics of the probability of observing a fixed number of photons in the cavity.In Fig. 5 we first consider an initial state with two photons, i.e. |( = 0)⟩ = |2⟩ ⊗ |∅⟩.The photon numbers exhibit coherent oscillations with near-perfect revivals and only a slow dephasing.In the large -limit the initial state only has a nonvanishing overlap with the triplet polaritons, following our previous discussion, such that the dynamics can effectively be treated as a three-level system.The period for revivals directly follows from Eq. ( 18) as At integer multiples of the period the system is to good approximation in a 2-photon states, whereas at half-integer multiples the system is close to a 0-photon state, reminiscent of the boson bunching in the Hong-Ou-Mandel effect [50].These coherent dynamics are now a direct consequence of the multifractality of the triplet polaritons: for delocalized eigenstates all eigenstates would be involved in the dynamics and all coherences would rapidly decay and thermalize, but since the initial (product) state only has  (1) overlap with the triplet polariton the system behaves as a three-level system with long-lived coherences.Following Eq. ( 43) and the surrounding discussion, this  (1) overlap directly relates to the change in IPR as  is varied and hence the multifractality.Second, we consider an initial state with a single photon excitation and a single spin excitation, |( = 0)⟩ = |1⟩ ⊗ Ŝ †  |∅⟩.We first observe that there are again coherent oscillations, now between the 1-photon and the 0-photon sector.These oscillations are now a direct consequence of the semilocalization of the doublet states in the 1-photon sector: the initial state will have  (1) overlap with an  (1) number of doublet polaritons, such that we can again restrict the dynamics to an  (1) number of doublet states.Additionally, the initial state has a vanishing overlap with both the dark states and the triplet polaritons: the former because of the vanishing weight of the dark state in the 1-photon sector, and the latter because of the delocalization of the triplet polaritons in the one-photon sector, leading to a vanishing overlap between the initial (localized) state and the delocalized eigenstate in this sector.Again, for purely delocalized eigenstates no such coherent dynamics would be observed and the system would rapidly thermalize.
We emphasize that these results conform to the expected cavity dynamics in the homogeneous limit, but here allow for a direct interpretation in terms of the localization properties of the different classes of eigenstates: In the presence of disorder, such coherent oscillations require nontrivial localization properties of the eigenstates within the different -photon sectors.The stability of the eigenstate localization properties to disorder indicates a stability of the dynamics of the homogeneous and integrable limits to the presence of disorder.
The delocalization of the triplet polaritons within the 1photon sector can also be directly observed in the vanishing probability of observing 2 photons in the cavity in Fig. 5.b), similar to the two-photon blockade [51][52][53].The overlap of the initial state with the triplet polaritons will scale as  ( −1/2 ) due to the delocalization within the -dimensional 1-photon Hilbert space.Since only these states have a nonvanishing contribution to the probability of observing 2 photons, the probability of observing a 2-photon state scales as  ( −1 ).Photon blockade can be attributed to the fact that the states with appreciable two-photon character are delocalized within each photon sector and thus two photons cannot be observed unless they already exist in the initial state (restricting to the reasonable conditions where atoms can be excited only locally).
This argument can be extended to probe the delocalization properties more generally.The time-averaged probability of observing two photons in the cavity for a generic initial state will be due to the overlap of the initial state with the triplet polaritons.For an initial state with 2 photon excitations we have already argued that this overlap is  (1), such that this time-averaged probability will be  (1).For an initial state with 1 photon excitation and a single spin excitation this probability scales as  ( −1 ), and for an initial state 0 photon excitations and two spin excitations this probability scales as  ( −2 ).These three different scales directly reflect the delocalization within different subspaces, and are illustrated in Fig. 6.In order to avoid averaging, we consider a so-called 'picket-fence' model of evenly spaced interaction strengths   .The different scalings can be clearly observed, relating the different scalings of the restricted IPRs for the triplet polaritons to a physical observable.

VIII. CONCLUSION AND DISCUSSION
We considered the Tavis-Cummings model in the presence of disorder in both the bare energies and the interaction strengths, focusing on the sector with two excitations.In the absence of disorder in the bare energies but for disorder in the interaction strengths, we derived an exact solution for the eigenstates and eigenvalues of the model.The model supports dark states, where the cavity mode does not hybdridize with the spins, and polariton states, where it does so.The dark states are known, and for the polaritons we show how the diagonalization of the full Hamiltonian can be reduced to solving a single secular equation, which can be numerically done in a straightforward way.Every eigenvalue corresponds to a solution of this equation, and we can identify different classes of eigenstates with different localization properties.
The main advances of this work are that we: (i) introduce a new solvable limit of the disordered Tavis-Cummings model, relevant for nonlinear photonics, presenting its exact solution and conservation laws, (ii) strengthen the connection between integrability and multifractality by analytically showing how the two integrable limits of this model exhibit multifractality, (iii) illustrate how this multifractality can be stable in the presence of disorder away from these integrable limits, and (iv) show how the multifractality can be better understood by introducing a restricted version of the inverse participation ratio, restricted to -photon sectors, which is in turn reflected in the dynamics of the photonic mode and is apparent in photon bunching and the two-photon blockade.
These different results indicate that, despite the presence of the disorder, the Hamiltonian can be effectively described in terms of collective spin operators.Such collective spin dynamics is expected for identical spins and naturally appears in models with sufficiently long-range interactions [35][36][37][38].While the spins in our setup are not identical due to the disorder, we showed for the new integrable limit that they can be treated as such up until relaxation time scales scaling as √ , with  the number of spin modes, provided that the energy is bounded in the limit  → ∞.Remarkably, the numerical results on the dynamics indicate that such a description remains accurate even in the presence of disordered bare energies.It would be interesting to further investigate when collective spin descriptions hold in the presence of disorder for cavity models and clarify the role of integrability, following similar studies for systems with long-range interactions [38].
The accessibility of exact eigenstates and eigenenergies has led to various studies of the dynamics in Richardson-Gaudin models [54][55][56][57][58][59][60], and this work opens up avenues to further study the dynamics in disordered Tavis-Cummings models.A natural extension of this work is to consider open systems.Here we only note that the model with inhomogeneous interaction strengths remains solvable if we choose Δ to be complex, leading to non-Hermitian and hence dissipative dynamics, since the calculation of the eigenstates did not depend on Δ being real.The dynamics generated by a non-Hermitian Richardson-Gaudin Hamiltonian can be directly studied, as e.g.done in Refs.[61,62], with only minimal modifications of the presented framework.An additional extension is to consider the model with an arbitrary number of excitations, where it is expected that many of the results presented in this work hold more generally.
Applying the Hamiltonian (54) to the state (55) and using the fact that  − |0⟩ = 0 gives: We will now focus on finding the commutator present above and begin by using the  (2) commutation relations of the spin operators to calculate: The commutator in (58) can now be rewritten as Applying this result to the vacuum state, substituting the definition of Ĝ+ , and relabeling dummy variables in summations, we end up with the result Using the definition of |⟩, our initial equation (58) becomes In order for |⟩ to be an eigenstate of the Hamiltonian, we require Ĥ () |⟩ =  |⟩, where  is the eigenvalue of the eigenstates.Imposing this on the above expression, we obtain the correct Bethe equation along with the expected eigenvalue for our initial state.

X. DERIVATION OF THE CONSERVED CHARGES
In this Appendix we explicitly derive the commutation relations of the conserved charges.The derivation is analogous to a similar derivation for the conserved charges of a spin-1 central spin model as presented in Ref. [32].These conserved charges are again constructed by using the properties of the integrable Richardson-Gaudin models (54), where Ĥ () has conserved charges with Q  defined in Eq. ( 20) in the main text.These charges satisfy for all ,  = 1, . . ., .The conserved charges from the main text are given in block matrix representation by where the different blocks correspond to different photon number.The Tavis-Cummings Hamiltonian can similarly be represented as a block-diagonal matrix Using only these block matrix representations, the commutator of the Hamiltonian with these charges can be evaluated as where we have not yet made use of any properties of the operators.The commutator on the diagonal vanishes since which can be checked either from direct calculation or by noting that Ĝ+ Ĥ ( = −1) = Ĥ ( = 1) Ĝ+ .This identity relates the eigenstates of Ĥ ( = −1) and Ĥ ( = 1), as also discussed in Refs.[28,32], and since the conserved charges share a common set of eigenstates with these Hamiltonians this identity should also hold on the level of the conserved charges, and we can rewrite the above equation as Since all matrix elements of the commutator vanish in the block matrix representation, we hence find that the Hamiltonian commutes with all proposed conserved charges.

XI. SECULAR EQUATION WITH VANISHING POLES
The pole structure of the secular equation ( 15), can change abruptly if the inhomogeneous interaction strengths are strongly asymmetrically distributed.Here ḡ2 =  =1  2

𝑖
and Δ = Δ/2.It is now possible for a pair of poles to vanish whenever which can be rewritten as This situation occurs whenever a single   is sufficiently large compared to the remaining interaction strengths.This equation can clearly only be satisfied for a single   , such that, in addition to the discussion of the main text, the only additional case that needs to be considered is the one where a single pair of poles vanishes.Both sides of the resulting secular equation are plotted in Fig. 7.There are again 2 + 1 intersections between both sides, indicating 2 + 1 solutions, such that we obtain the corect number of states.While the total number of solutions remains the same, the direct identification between solutions to the secular equation and eigenstates in the homogeneous limit now breaks down.The pole structure results in 2 − 4 solutions in between pairs of poles (marked by circles), and there are 5 additional solutions.The triplet polariton solutions can again be identified near the vertical asymptote of the lefthand side and for large  (as marked by squares), but two additional solutions now appear in the middle interval (marked by crosses).These additional solutions always appear, since the right-hand side of Eq. ( 69) lies above the horizontal asymptote of the right-hand side for  around zero.
Let us focus on the center interval where the two additional poles appear.The number of intersection will be determined by the behavior of the right-hand side for  → 0. In the case from the main text no poles vanish and it directly follows that for  → 0 the right-hand side will always be negative and hence below the horizontal asymptote 1/3 of the left-hand side, such that there are no additional intersections with the left-hand side.If a pole vanishes, then for  → 0 the right-hand side can be shown to be larger than 1 and hence above the horizontal asymptote 1/3, introducing two additional intersections.
That the right-hand side is always larger than 1 for  = 0 can be directly checked.Assuming that the pole corresponding to  1 vanishes, the right-hand side for  = 0 can be written as where  2 = Δ2 +  =2  2  such that  2 1 ≥  2 (in order for the pole to vanish) and  2  1 ≥  2  .For  2 1 → ∞ this expression approaches 1, for  2  1 →  2 from above this expression diverges to +∞.The derivative w.r.t. 2  1 can easily be checked to be negative for all values of  2 1 in between these two limits, such that this expression is monotonous between these limits.

XII. IPR SCALINGS FROM THE BETHE ANSATZ
The multifractality scalings can be directly obtained from the Bethe state (30) and reflect the different positions of the solutions w.r.t. the poles.From Eq. (30) we have that where ψ equals  up to a global normalization factor, since the Bethe states are unnormalized.
We can now identify different scaling behaviors depending on the relative positions of the variables  1,2 w.r.t. the poles   ,  = 1 . . ., in the same way that the relative position of the root to the secular equation (15) w.r.t. the poles   ,  = 1 . . . allowed us to distinguish different classes of eigenstates.The derivations for both integrable limits are highly similar, but since the IPR cannot be defined for the degenerate dark states in the limit of homogeneous bare energies, we first focus on the limit of homogeneous interaction strength and inhomogeneous bare energies.
Triplet polaritons.Consider a situation where the two Bethe roots  1,2 are a distance  (1) away from all poles   ,  = 1 . . ..In this scenario all terms 1/(  −   ) will be  (1), such that the scaling of the amplitudes is purely set by the Each first term in a product is the number of components and the second term is the scaling of the individual components.Crucially, we find that P 0 ( = 1) = P 1 ( = 1) = P 2 ( = 1) =  (1), such that ψ and  have the same scaling.These results then reproduce the observed scalings from the main text.Doublet polaritons.In this case a single Bethe root, e.g. 1 , lies in between two poles, where the typical distance between

FIG. 1 .
FIG. 1. Eigenspectrum of the clean and disordered Tavis-Cummings Hamiltonian for two excitations.In the absence of disorder three classes of eigenstates can be observed resulting in dark states (singlet states) and polaritons (doublet and triplet states).For disordered systems (gray lines) the dark states and the doublet polaritons split into bands of states, while the triplet polaritons remain isolated.Parameters:  = 40,   ∈ [−0.1, 0.1] and   ∈ [1, 2] uniformly distributed for the disordered model and  =   and  2 =  2 for the homogeneous model.

2 FIG. 4 .
FIG. 4. Illustration of the one-photon wave function components |  | 2 for representative polariton (triplet and doublet) and dark (singlet) states.For homogeneous couplings these components are a smooth function of the corresponding bare energy   , with different states characterized by the number of poles where the amplitudes scale as |  | 2 ∝ 1/(  −   ) 2 .In the presence of disorder in the couplings this overall structure is preserved and in turn reflected in the IPR.Parameters:  = 100, Δ = 1,   uniformly distributed in [−0.1, 0.1] and   uniformly distributed in [1, 2]/ √ .

FIG. 5 .
FIG. 5. Dynamics of the probability of observing  photons in the cavity for an initial state with 2 photons (a) and an initial state with 1 photon and a single cavity excitation (b).Parameters:  = 40,   uniformly distributed in [1, 3]/ √  and   uniformly distributed in [−0.2, 0.2].

4 FIG. 7 .
FIG. 7. Graphical illustration of the secular equation in the case where a pair of poles vanishes.Each intersection between the left-hand side (red line) and right-hand side (blue line) returns the eigenvalue of a polariton state, leading to different classes of solutions.Parameters: Δ = 1,  = 0,  = 5 and  2  ∈ [2, 4, 6, 8, 32].