Polariton Creation in Coupled Cavity Arrays with Spectrally Disordered Emitters

Integrated photonics has been a promising platform for analog quantum simulation of condensed matter phenomena in strongly correlated systems. To that end, we explore the implementation of all-photonic quantum simulators in coupled cavity arrays with integrated ensembles of spectrally disordered emitters. Our model is reflective of color center ensembles integrated into photonic crystal cavity arrays. Using the Quantum Master Equation and the Effective Hamiltonian approaches, we study energy band formation and wavefunction properties in the open quantum Tavis-Cummings-Hubbard framework. We find conditions for polariton creation and (de)localization under experimentally relevant values of disorder in emitter frequencies, cavity resonance frequencies, and emitter-cavity coupling rates. To quantify these properties, we introduce two metrics, the polaritonic and nodal participation ratios, that characterize the light-matter hybridization and the node delocalization of the wavefunction, respectively. These new metrics combined with the Effective Hamiltonian approach prove to be a powerful toolbox for cavity quantum electrodynamical engineering of solid-state systems.

Nanophotonic cavities with integrated quantum emitters have served as a rich playground for exploring quantum optics phenomena in solid-state systems.This includes demonstrations of weak [14] and strong [15] cavity quantum electrodynamical (QED) coupling, photon blockade and photon-induced tunneling [16], ultra-fast modulation of optical signals [17], and more.The large dipole moment of quantum emitters, paired with (sub)wavelength scale optical mode volumes in photonic crystal cavities, give rise to high optical nonlinearities and light-matter state hybridization that creates polaritons.Polaritonic interactions in nanophotonic systems can be several orders of magnitude higher than those achieved in atomic systems.Such strong interaction has been at the core of theoretical proposals for quantum state transfer [18,19], as well as for the photonic simulation [20,21] of Bose-Hubbard and fractional quantum Hall physics.Here, the system is made of an array of coupled cavities, each in the strong coupling regime of cavity QED, and described by the Jaynes-Cummings-Hubbard model.However, this model has been experimentally hard to achieve.
While progress toward the realization of coupled cavity arrays (CCAs) with embedded emitters has been made with quantum dots [22,23], the spectral disorder of these emitters has been a major roadblock to developing a large-scale resonant system.This problem is not present to such an extent with color center emitters, which are atomic defects in wide band gap materials.Recently, color center integration with nanocavities in diamond [24,25] and silicon carbide [26,27] has been demonstrated in the weak cavity QED coupling regime.
Though this regime is unsuitable for studies of polaritonic physics, proposals to demonstrate strong cavity QED regime have been presented with cavities integrating several (M ) emitters, as opposed to a single emitter.Additionally, there has been renewed interest in disordered cavity QED systems with the discovery of phenomena like collectively induced transparency [28].Such systems are described by the Tavis-Cummings, rather than the Jaynes-Cummings model.Here, the collective coupling of emitters to the cavity effectively boosts the lightmatter interaction rate by a factor of √ M .Due to the small, but nonzero, spectral disorder of color centers, the collective strong coupling is possible within the cavity protection regime, if its rate overcomes the spectral disorder ∆ of color centers [29,30], i.e. ∆ < g √ M .Such disordered multi-emitter cavity systems have been explored for applications in quantum light generation [31][32][33].
Here, we explore how all-photonic quantum simulators based on coupled cavity arrays can benefit from an increased interaction rate established in multi-emitter cavity QED.
We expand the Jaynes-Cummings-Hubbard approach to the spectrally disordered Tavis-Cummings-Hubbard model (TCHM) [34] and define conditions for polariton creation utilized in all-photonic quantum simulation, aided by the introduction of new localization metrics inspired by condensed matter approaches.Our model targets applications in technologically mature solid-state platforms and is reflective of the state-of-the-art parameters achieved in silicon carbide and diamond color center hosts.We find system limits that can be guiding for future experiments with polaritons in coupled cavity arrays: 1) polaritonic states are easier to create in systems where emitter/cavity interaction exceeds cavity hopping; 2) polariton creation in an array of lossy cavities can be achieved via integration of an increased number of emitters per cavity even in disordered ensembles; 3) as in single cavities, disordered emitter ensembles in coupled cavity arrays can create polaritons by increasing the number of emitters per cavity to reach the cavity protection condition; and 4) disorder in resonances in a coupled cavity array localizes polaritons if the difference in frequencies between neighboring cavities exceeds the cavity hopping rate.

II. THE CCA QED MODEL
Our CCA QED model captures the single-excitation regime of the spectrally disordered TCHM comprised of emitter-cavity localizing interactions and cavity-cavity delocalizing interactions: where N is the number of cavities in the array, M n is the number of emitters in the n-th cavity, ω c,n and a n represent the angular frequency and the annihilation operator of the n-th cavity, ω e,n,m , σ − n,m and g n,m correspond to the angular frequency, the lowering operator and the emitter-cavity coupling rate of the m-th emitter in the n-th cavity, J n,n+1 is the photon hopping rate between the enumerated neighboring cavities.In this work we will assume J n,n+1 = J as those parameters are set by the cavity design that is experimentally more controllable than the other parameters of the model [35].
A. Non-disordered CCA QED model Before examining the spectral disorder effects in CCA QED, let us first address the energy spectrum in the fully resonant system of a linear array of coupled cavities with identical emitters.Here, the eigenenergy spectrum features two CCA polariton bands with N states, each, and a degenerate set of N (M − 1) subradiant states as illustrated in Figure 1.The polariton band states are parameterized by discrete momenta k = k p = πp/(N + 1), (p = 1, 2, 3, .., N ) as ( Further discussion of the derivation of these equations for the fully resonant case is presented in Section 1 of the Supplementary Information.The origin of these spectral features can be decomposed to the QED and the CCA components.The resonant Tavis-Cummings model of M emitters in N = 1 cavity has the spectrum of two polaritons and M − 1 degenerate subradiant states, while a CCA of N > 1, M = 0 cavities has a single spectral band of N photonic states.The spectrum of the resonant TCHM is a product of these components as seen in Figure 1.This model has close analogs to the condensed matter models it aims to simulate.

B. Condensed matter analogs of TCHM
A more thorough discussion of these analogs and their limitations is available in the Supplementary Information, but we can begin by looking at the TCHM system with identical emitters (ω e,n,m = ω e ).Much of the derivation of eigenstates precisely parallels the calculation of the band structure of tight binding Hamiltonians commonly studied in condensed matter physics.For example, the case of no emitters (M n = 0) corresponds to the d = 1 (Bose) Hubbard Model (HM) [36], and the case with a single emitter in each cavity (M n = 1) to the Periodic Anderson Model (PAM) [37]: the hopping of photons between cavities is analogous to the conduction electrons which hybridize between different sites, and the photon-emitter coupling maps onto the hopping between the conduction electrons and localized orbitals which, like the emitters, do not have a direct intersite (intercavity) overlap.Thus, in the single excitation sector, and for one emitter per cavity, the models are identical.In practice, the polaritonicity (degree of light-matter hybridization) upon which we focus in the sections to follow, holds lessons for singlet formation in which local and itinerant electrons become tightly intertwined in the PAM and Kondo Lattice Model.

C. Experimentally-informed TCHM simulation parameters
The parameters of our model have been selected as representative of silicon carbide and diamond color center platforms.Recent demonstrations of emitter-cavity interaction in photonic crystal cavities support rates of approximately g/2π ∼ 2-7.3 GHz [26,38], therefore, we chose a constant value of g/2π = 5 GHz.While a variation in the coupling rate g among emitters is likely to occur due to their variable positioning inside the electromagnetic mode, our prior work indicates that the collective emitter-cavity coupling still takes place [29] at a well defined rate of g M = M m=1 g 2 m .Therefore, keeping g constant among the emitters should not take away from the overall phenomenology studied here.
The experimentally demonstrated cavity loss rates reach as low as κ/2π ∼ 15-50 GHz [26,38], while recent modeled designs could reduce these values by at least an order of magnitude [35].With a slight optimism, we chose cavity loss rate of κ/2π = 10 GHz.Our recent designs of photonic crystal molecules indicate that coupled cavity hopping rates can be straightforwardly designed in the range 1 GHz < J/2π < 200 GHz [35], thus spanning systems from the dominant cavity QED to the dominant photonic interaction character, represented in our choice of values J/g = 0.1, 1, 10.
Fabrication imperfections may yield drifts in cavity resonant frequencies and hopping rates.The effect of this issue was studied in another platform where GaAs coupled cavity arrays were integrated with quantum dots [23] and indicates that the coupling strength is an order of magnitude higher than the frequency and hopping rate perturbations.We apply this assumption in our model, maintaining that all cavities are mutually resonant and all rates J are constant.
The spectral inhomogeneity of emitters in fabricated devices, the main study of our model, has been characterized as ∆ ∼ 10 GHz for a variety of emitters in silicon carbide and diamond [39,40].We represent this parameter through its relation to the collective coupling rate g √ M in a cavity, spanning the spectral inhomogeneity across a range of values.The spectral disorder is implemented by sampling emitter angular frequency, ω e , from a Gaussian distribution } centered at ω c with a width of ∆ = 2σ.It is worth noting that the vibronic resonances are three orders of magnitude larger than the inhomogeneous broadening, for example 8.7 THz for the silicon vacancy in 4H-SiC [41,42], therefore the phonon side band is not expected to play a part in the collective emitter-cavity coupling process.Emitter lifetime in color centers is usually in the 1-15 ns range [43], we select the value of the spontaneous emission from the color center, γ/2π = 1/5.8GHz, as representative.Due to γ being the lowest rate in a color center-based cQED system, its minimal variations among emitters of the same species [40] affect the system only marginally, therefore we assume it has a constant value; this is representative of systems like atoms and color centers.

Lifetime-and nearly lifetime-limited emission of color centers has been demonstrated
upon photonic integration [40,44,45].Due to this experimental advance, our model does not consider the dephasing terms, though such analysis may prove valuable with further development of integrated coupled cavity arrays.
With these experimental constraints in mind, we believe our simulations will be directly relevant to future fabricated multi-emitter photonic crystal cavity chains.

III. EFFECTS OF DISORDER ON POLARITON FORMATION
It is in general computationally expensive to solve the Lindbladian master equation to obtain exact simulation results [46].As such, we are restricted to simulating only very small scale systems (∼ 6 elements total, an element being a single cavity or emitter) even in the low excitation regime.The results of these exact simulations are available in Sections 3 and 4 of the Supplementary Information; the matrix form of the Hamiltonian that describes the system simulated is available in Section 2 of the Supplementary Information.On the other hand, the effective Hamiltonian (H EFF ) uses the established non-hermitian effective approach to modeling Hamiltonians that are too resource intensive for the current state of the art classical computers to solve.Its approximation effectiveness is limited to the singleexcitation regime, which is suitable for our exploration.Taking the effective Hamiltonian approach a step further, we also introduce the nodal and polaritonic participation ratios (P N and P P , respectively).This method is derived from the condensed matter participation ratio metrics [47] used to quantify a system's localization properties.The P N and P P metrics are applied to eigenstates of the H EFF to quantify the delocalization and light-matter hybridization of the wavefunction, respectively.
To access modeling of larger systems we develop a software package in Python [48] that diagonalizes the Effective Hamiltonian in the approximate single-excitation regime thus reducing the computational complexity from exponential to polynomial (cubic) in N × (M + 1) for the single-excitation regime.With this approximate method we numerically solve systems with hundreds of elements compared to the several using the exact QME approach.Note that, in contrast to QME, this method does not contain a pump term, meaning it is agnostic to the starting cavity and its diagonalization will provide all possible states, regardless of their wavelength overlap with the initial cavity.
A. The Participation Ratio Approach: Metrics for characterizing disorder The node-by-node and element-by-element analysis required to examine each of the eigenstates found using H EFF in the previous sections is lengthy and not suitable for the much larger systems we will be exploring.In order to efficiently analyze these much larger systems, we develop new metrics for the characterization of TCHM wavefunctions, inspired by practices in Condensed Matter Physics.The phenomenon of Anderson localization describes the loss of mobility of quantum particles due to randomness [49].Originally studied in the context of non-interacting electrons hopping on a lattice with disordered site-energies, where all eigenstates were shown to be localized in spatial dimension less than or equal to two [50,51], Anderson localization has subsequently been extensively investigated in many further contexts, including the effect of interactions [52], correlations in the disorder [53], and importantly, new experimental realizations from cold atomic gases [54,55] to transport in photonic lattices [56].
A useful metric for quantifying the localization of a wavefunction v p , employed in these studies, is the participation ratio, P = p |v p | 4 −1 [57] and its generalizations [58].
Instead of measuring the participation ratio among all N (M + 1) vector components, we adapt P to two new metrics that measure the participation among N nodes (cavity-emitter sets), and two cavity-and emitter-like components.We define the nodal participation ratio, where N ph,n = a † n a n and N e,n = Mn m=1 σ + n,m σ − n,m are the usual number operators for each state v p representing cavity excitation and the sum of all emitter excitation in a cavity.Like the classic participation ratio, P N is at a minimum (maximum) when the wavefunction is localized (delocalized).Next, we define the polaritonic participation ratio, or polaritonicity, which is minimized when the wavefunction has completely cavity-like or completely emitterlike character and is maximized for an equal superposition of cavity-and emitter-like components.Here we assume the character is polaritonic when there is any type of light-matter hybridization whether or not it is coming from the same node and therefore note that a wavefunction can be polaritonic even when the cavity and emitter excitations do not belong to the same node.
These two new metrics allow us to seamlessly characterize multi-emitter CCAs.We normalize the metrics to 1 for easy comparison between the different model parameter cases.
To avoid numerical divide by zero errors, we set the identical emitter case of the leftmost column to have a small but nonzero value (∆ = ϵ ≈ 10 −7 ). of the eigenstates of Eq. 3. In this section we set ω c,n = ω c and g n,m = g.Figure 2 explores the regime where cavity QED dominates the photon hopping J/g = 0.1.For vanishing spectral disorder, ϵ, the eigenspectrum has the shape resembling the features of Figure 1: two highly polaritonic delocalized CCA bands with N = 5 states and N (M − 1) = 10 highly localized subradiant states, suitably characterized by the polaritonic and nodal participation ratio values.For moderate disorder, the polaritonic properties of eigenstates are maintained, while the nodal localization somewhat increases for polaritonic band states.The degeneracy of the subradiant states is lifted and the spectral gaps diminish as we move into the strong disorder regime wherein ∆ ≈ g M , which is usually considered a cutoff for cavity protection.
Most states become highly localized, demonstrated by the significant drop in P N value and the subradiant states gain a cavity component, as quantified by the increase of the P P value.While similar trends can be observed, the main difference is seen in the reduction of the number of polaritonic states in the CCA bands as the wavefunction obtains a higher cavity-like character.
This brings us to look into the formation of a polaritonic state in spectrally disordered CCA QED as a function of an increasing J/g ratio.Figure 3 shows the polaritonicity and localization of the lowest energy eigenstate.
When the photonic nature of the interaction increases, so does the cavity-like character of the wavefunction, reducing its level of polaritonicity.While this holds true for low and moderate values of disorder, in the case of high disorder, we observe an increase of P P , before the decline.This is an artifact of the disorder which randomly modifies the nature of the lowest eigenstate in the system to be more emitter-like, until the interaction value increases to a level that offsets the issue.This trend is paired with the increase in the delocalization metric P N as the wavefunction loses the dominant emitter-like characteristic.A decrease in polaritonicity and delocalization of the wavefunction take place for a range of system parameters.At low J/g there is less variance in P P for a larger ∆ compared to a higher photon hopping rate, suggesting that, as in the single node Tavis-Cummings model [29], the stronger cavity-emitter coupling compared to combined cavity losses (in the TCHM this includes cavity-cavity coupling) provides better cavity protection against the disorder in the TCHM.In the top plot, the coupling strength, g/2π, of each emitter is randomly taken from a Gaussian centered at 5 MHz with FWHM ∆ g that is limited to values between zero and ten.In the bottom plots ω c /2π of each cavity is pulled from a Gaussian centered at 0MHz = ω e with standard deviation ∆ c /2.The nodal localization effects of increasing ∆ g and ∆ c are comparably less than those seen in Figure 3 for increasing ∆ e .

C. Other experimentally relevant forms of disorder
While we expect disorder due to spectral inhomogeneity of the emitters to dominate the creation of polaritons in the TCHM, we explore the effects of other potential sources of disorder that will arise in experimental implementations of these systems.
One potential source is the disorder of the emitter-cavity coupling rate, ∆ g .This disorder is important experimentally to consider since the dipole direction and the exact placement A second experimentally significant source of disorder is the variation in the cavity frequency from one cavity to the next, ∆ c .Any minor variance in nanofabrication from one cavity to the next will alter ω c [13].The decrease in P N as disorder in ω c increases in Figure 4b (top) suggests that for an increasing cavity frequency disorder, the lowest energy eigen-states become more localized.In fact, Figure 5a (top) suggests that there will always be a drop in P N as ∆ c increases and that regardless of the strength of the cavity-cavity coupling, the rate at which P N decreases is the same.On the other hand, Figure 5b (bottom), like Figure 5a (bottom), suggests that the polaritonicity of the system is fairly tolerant of cavity fabrication errors.
Both ∆ c and ∆ g have similar effects on P P , the scale of which is set by the value of J.In contrast to variations in ω e that effect both P P and P N , variations in ω c and g affect only the localization properties of polaritons.In practice, this means that if neighboring cavities in the CCA have sudden jumps in ω c or g it can effectively cut the CCA into two smaller CCAs with independent TCHM physics from one another.

IV. DISCUSSION
In this work we characterize the CCA QED eigenstates described by the Tavis-Cummings-Hubbard model.Our goal is to provide a guiding tool for experimental implementations through the engineering of the CCA parameters.
Using the new participation ratio metrics, inspired by condensed matter physics studies of localization and band mixing, we confirm that highly polaritonic states can be formed in coupled cavity arrays despite the presence of spectral disorder in emitter ensembles and quantify the cavity protection effect.While the systems with a dominant cavity QED interaction, relative to the photon hopping rate, support creation of numerous polaritonic states, we find that other parts of the parameter space can also be utilized to study polaritonic physics.
We suggested approximate analogies between the case of M n = 1 emitter in each cavity with the periodic Anderson model where a single f orbital on each site hybridizes with a conduction band, and the Kondo lattice model where the local degree of freedom is spin-1/2.Condensed matter systems which connect to the multi-emitter case M n > 1 also have a long history, both in the investigation of multi-band materials and also as a theoretical tool providing an analytically tractable large-N limit [60,61].Indeed, large-N systems, realized for example by alkaline earth atoms in optical lattices, are also at the forefront of recent work in the atomic, molecular and optical physics community [62][63][64].In short, the TCHM offers a context to explore intertwined local and itinerant quantum degrees of freedom which, while distinct from condensed matter models, might still offer insight into their behavior.
These TCHM systems can ostensibly be realized in a number of photonic frameworks, from atoms in mirrored cavities, to quantum dots in nanophotonics.It is difficult, however, to experimentally create atom-based systems that couple multiple cavities together and to create large numbers of quantum dots that emit within the relatively modest range of disorder that we have shown will recreate polariton dynamics.As such, the most likely experimental realization of our systems will be in color center based nanophotonics.

II. MATRIX FORM OF H T CHM
As an illustrative example, the Hamiltonian of the largest system simulable on an M2 Macbook Pro with 32 GB of memory using the quantum master equation approach (2 cavities with 2 emitters in each cavity) is presented in matrix form: For convenience the equation form is reproduced here: The ω c,n a † n a n term gives the energy of the nth cavity, similarly the ω e,n,m σ † n,m σ n,m term gives the energy of the mth emitter in the nth cavity.The interaction between this emitter and cavity is described by the coupling term: g n,m a † n σ n,m + σ † n,m a .Finally, the interaction between nearest neighbor cavities is given by the last term, J n+1,n a † n+1 a n + a † n a n+1 .

III. EXACT SIMULATION RESULTS
We have developed a software package [1] that solves the Quantum Master Equation (QME).Our code uses the Quantum Toolbox in Python (QuTiP) [2] which solves the Lindbladian where is the cavity linewidth of the n-th cavity, γ n,m = γ is the emission rate of the m-th emitter in the n-th cavity, and P is the optical (laser) pumping term.The spectral intensity reported in Figure 2 is calculated as the Fourier Transform of the correlation function ⟨A † (t + τ )A(t)⟩: where A is replaced by the cavity annihilation and the emitter lowering operators, a n and The QME requires the use of the full density operator because of the non-number conserving term 2cρ(t)c † .
The density operator requirement makes solving this system highly resource-intensive (exponential in N × M ), as such we have restricted our exact calculations to small systems of six elements, or two coupled cavities with two emitters per cavity.Systems with more processing power may be able to simulate larger systems.Quantum trajectory theory is another potential method that may simulate larger systems using the quantum master equation approach.
We numerically solve the quantum master equation for a coupled 2 cavity system with each cavity also coupling to two emitters each.The results are shown in Figure 2. Along the top row, which shows the smallest inter-cavity coupling simulations, we see the most nodal localization regardless of emitter dispersion, as shown by only having peaks in Cavity 1 and Emitters 1 and 2 when pumping Cavity 1 and only having peaks in Cavity 2 and Emitters 3 and 4 when pumping Cavity 2. We can also see that the small inter-cavity coupling leads to highly polaritonic systems by comparing the peak heights of each cavity to those  of its emitters and finding the ratios are close to unity (0.78 -0.9).The largest inter-cavity coupling, J/g = 10 is heavily cavity-like as shown by cavity peak to emitter peak ratios of around 50.The middle row, J/g = 1 is partially polaritonic and partially photonic with peak ratios in the range 0.3-0.6.By introducing nonzero emitter detuning, we are able to access subradiant states for each of the three detuning values.These states are identified by the zero-width peaks found between the polariton peaks.Subradiant states are states that decay much slower than the polariton peaks and a number of proposals exist for their use in quantum information technologies including, light storage [3] and quantum light generation [4].

IV. BENCHMARKING OF THE EFFECTIVE HAMILTONIAN APPROXIMA-TION AND THE PARTICIPATION RATIO METRICS
To benchmark the approximate Effective Hamiltonian approach against the exact Quantum Master Equation solution, we simulate the systems with same parameters in both models.Figure 3 shows the lowest energy wavefunction calculated with the Effective Hamiltonian approach, corresponding to the excitations of the lowest energy peaks in Figure 2.
The case of the vanishing spectral disorder shows equal contribution of Cavities 1 and 2, which corresponds to the identical looking plots of the Cavity 1 and Cavity 2 excitations in the QME spectra when Cavity 1 and Cavity 2 are pumped, respectively.The asymmetry of the nodal occupations in Figure 3 for non-vanishing disorder (especially for J/g = 0.1) matches the non-identicality of the Cavity 1/2, as well as Emitter 1.1,1.2 and Emitter 2.1,2.2spectra in Figure 2. Localization of the wavefunction for an increasing disorder and J/g ≤ 1 follow the exact solution trends described in the previous section.
These parallels are closely described by the nodal P N and the polaritonic P P participation ratio shown in Figure 4.The value of P N (orange) is maximized for a fully node-delocalized wavefunction, and reduces as the wavefunction tends to increasingly excite one cavity and its emitters with an increasing disorder.The value of P P (green) is maximized for the wavefunctions that have equal excitation distribution between cavities and the emitters, and reduces with an increasing photonic interaction (high J/g value) as the cavities become predominantly excited.
We conclude that the Effective Hamiltonian and the participation ratio metrics are suitable for studies of the Tavis-Cummings-Hubbard model.

V. EFFECTS OF SPECTRAL DISORDER ON THE TCHM ENERGY SPEC-TRUM
We utilize the Effective Hamiltonian approach to study medium and large TCHM systems and are especially concerned with the influence of the spectral disorder on the polaritonic and localization properties of the model wavefunctions.The Figure 6 shows the energy spectra of an N = 5, M = 3 system for various regimes of cavity QED to photon hopping ratios.
For the vanishing disorder, we clearly see the three components of the spectrum illustrated in Figure ??: two polaritonic bands with N states and a subradiant band with N (M − 1) degenerate states.When cavity QED interaction is significant, a band gap opens between the polariton bands and the subradiant states.With an increasing spectral disorder, the band gap closes and the subradiant state degeneracy is lifted.
Further analysis is provided by the corresponding P N and P P values shown in Figure 6.Here, we observe that the polaritonic properties of the highly hybridized states in the polariton bands and the emitter-like states in the subradiant band shift significantly only for the high levels of spectral disorder, defined by the typical cavity protection cutoff ∆ = g √ M .The increasing localization trend (decreasing P N ) for most polariton band states with an increasing disorder is evident for all sets of parameters.In contrast, the subradiant states gain a cavity component with an increasing disorder and become more hybridized and delocalized.
While the polaritonicity of the lower (as well as the upper) polariton band reduces for most states with an increasing J/g ratio, the middle state of the band remains highly polaritonic.
This state, labeled the most polaritonic state (MPS) in our study, shows that even the systems with high hopping ratio can serve as testbeds for polaritonic physics explorations.

VI. A PROPOSED IMPLEMENTATION
A physical implementation of the described coupled cavity array system using atoms trapped in optical cavities is depicted in Figure 7.In such a system, the excitation photon would come from a narrow-band tunable laser with low enough power to only target the first excited rung of the Jaynes-Cumming ladder which is what is studied in this work.The ability to collect photons from each cavity will give information about state transfer and localization properties.[4] M. Radulaski, K. A. Fischer, K. G. Lagoudakis, J. L. Zhang, and J. Vučković, Photon blockade in two-emitter-cavity systems, Physical Review A 96, 011801 (2017).

FIG. 1 .
FIG. 1.The eigenspectra of non-disordered coupled cavity arrays.(Left) One cavity with M emitters has two polariton states and M − 1 subradiant states.(Middle) CCA of N cavities with no emitters has one CCA band of N states.(Right) CCA of N cavities and M emitters per cavity has two polariton bands, upper and lower, of N states each and N (M − 1) subradiant states.

[ 59 ]
FIG. 2. Effects of increasing disorder, ∆, on the eigenstates of a system of N = 5 cavities and M = 3 emitters per cavity with J/g = 0.1.(Top) Energy eigenvalues and (bottom) participation ratios P N , P P .Increasing ∆ lifts the degeneracy of the subradiant states (flat middle band) decreasing the gap between the polaritonic bands and the subradiant states.Increasing ∆ also causes an Anderson-like node localization in the polariton bands as shown by P N decreasing, but the polaritonic states remain mostly polaritonic.Mean of 100 random realizations; error bars are one standard deviation.g, κ, γ are detailed in Section II C.

FIG. 3 .
FIG. 3. The nodal and the polaritonic participation ratios for the lowest energy eigenstate of a CCA with N = 5 and M = 3 for increasing J/g.A small amount of disorder causes the state to node localize for small J/g.Mean of 100 random realizations; error bars are one standard deviation.N = 5, M = 3; g, κ, γ are detailed in Section II C.

FIG. 4 .
FIG. 4. The nodal and polaritonic participation values for the lowest energy eigenstate of a CCA with N = 5 and M = 10 showing for increasing J/g and increasing values of disorder in g (top) and ω c (bottom).Means of 100 runs; error bars are one standard deviation.In the top plots, g i /2π of each emitter is randomly taken from a Gaussian centered at 5 MHz with standard deviation ∆ g /2

FIG. 5 .
FIG. 5. Plot (a) shows the average P N and P P for the lower energy eigenstate over 100 runs for systems with increasing ∆ c .Plot (b) shows the same but for increasing ∆ g .Parameters are the same as those reported in Figure 4. Various values of J/g are denoted by the different plot markers in each plot.

FIG. 1 .
FIG. 1.The two polariton bands of a cavity-emitter system for representative parameters ω c = ω e ≡ ω 0 = 4, cavity hopping J = 1, and cavity-emitter coupling g = 0.2.If the number of emitters M > 1 the gap between the two bands is enhanced from ∆ = 2g to ∆ = 2 √ M g and there are M − 1 additional flat emitter bands at E α (k) = ω 0 in the gap between the two polariton bands.

FIG. 2 .
FIG. 2. The top plots show the pump operator acting on Cavity 1 and the bottom set show the pump term acting on Cavity 2. Each line within a subplot refers to the transmission spectrum acquired by probing a particular element of the systems, either Cavities or any of the Emitters.Cavity 1 (blue) is coupled to Emitters 1,1 and 1,2 (green and red respectively).Cavity 2 is coupled to Emitters 2,1 and 2,2 (purple and brown respectively).The rows correspond to different cavity-cavity to cavity-emitter coupling ratios, with J/g = 0.1 (top row) having the least coupling between cavities, and J/g = 10 (bottom row) having the most.Moving from left to right increases the amount of disorder introduced via emitter emission wavelength detuning from a nominally zero amount (ϵ), to detuning approximately equal to the cavity-emitter coupling (g).

FIG. 3 .
FIG. 3. Node occupancy of the lowest energy eigenstate of a system with N = 2 cavities and M = 2 emitters per cavity with random ω e sampled from a Gaussian distribution of width ∆.Single random realization for each ∆.g, κ, γ, and emitter frequencies are the same as those in

FIG. 4 .
FIG. 4. Participation ratios of the lowest energy eigenstate for a system with N = 2 cavities and M = 2 emitters per cavity with random ω e sampled from a Gaussian distribution of width ∆.Single random realization for each ∆.g = 5, κ = 10, γ = 1/5.8.The emitter frequencies are the same as those in Figure 2.

FIG. 7 .
FIG.7.An outline of a potential implementation of a polaritonic simulator for condensed matter systems.The excitation photon (left, red arrow) enters a mirror cavity (grey semi-ovals).The cavity is connected to a number of emitters (red circles).The cavities are coupled to each other with another two-way coupling constant (red double-headed arrows).After excitation, each cavity has a probability of emitting a photon (red vertical arrow) at a rate κ which then gets collected by an objective (blue) sent to a detector.