Non-equilibrium many-body effects in driven nonlinear resonator arrays

We study the non-equilibrium behavior of optically driven dissipative coupled resonator arrays. Assuming each resonator is coupled with a two-level system via a Jaynes-Cummings interaction, we calculate the many-body steady state behavior of the system under coherent pumping and dissipation. We propose and analyze the many-body phases using experimentally accessible quantities such as the total excitation number, the emitted photon spectra and photon coherence functions for different parameter regimes. In parallel, we also compare and contrast the expected behavior of this system assuming the local nonlinearity in the cavities is generated by a generic Kerr effect rather than a Jaynes-Cummings interaction. We find that the behavior of the experimentally accessible observables produced by the two models differs for realistic regimes of interactions even when the corresponding nonlinearities are of similar strength. We analyze in detail the extra features available in the Jaynes-Cummings-Hubbard (JCH) model originating from the mixed nature of the excitations and investigate the regimes where the Kerr approximation would faithfully match the JCH physics. We find that the latter is true for values of the light-matter coupling and losses beyond the reach of current technology. Throughout the study we operate in the weak pumping, fully quantum mechanical regime where approaches such as mean field theory fail, and instead use a combination of quantum trajectories and the time evolving block decimation algorithm to compute the relevant steady state observables. In our study we have assumed small to medium size arrays (from 3 up to 16 sites) and values of the ratio of coupling to dissipation rate $g/\gamma \sim 20$ which makes our results implementable with current designs in Circuit QED and with near future photonic crystal set ups.


Introduction
Coupled resonator arrays (CRAs) interacting with single two-level systems embedded in each resonator have recently emerged as an exciting new platform for the realization of novel quantum many-body effects. CRAs may offer several features complementary to those of the successful and well-established 'toolbox' of cold atoms in optical lattices [1,2], such as single-site addressibility and intrinsic non-equilibrium physics. This has led to proposals to realize a variety of phenomena of great interest in condensed matter physics such as Mott transitions [3,4,5], effective spin models [6] and fractional quantum Hall states [7], among others [8,9]. Though promising efforts are underway to build the first CRA systems, a number of technical challenges originating from the existence of strong dissipation must be overcome before equilibrium physics can be explored. On the other hand, the natural open nature of CRAs and the inherent ability to address and observe single resonators make this system an ideal test-bed to study many-body quantum lattice models out of equilibrium beyond the canonical Bose-Hubbard (BH) model and its realizations in optical lattices.
In this work we look for signatures of underlying many-body phenomena in CRAs by analyzing observables measured in the system's non-equilibrium steady state (NESS). These include optically accessible observables like photon spectra and correlation functions. We focus on small to medium size CRAs of a few sites implementable with current or near future experimental technologies. We study coupled single mode resonators each interacting with a two-level system, a setup known now as the Jaynes-Cummings-Hubbard (JCH) model [3]. This simple model, in contrast to proposals involving multi-level atomic systems and external fields [4], may be realized in a variety of technologies ranging from quantum dots [10,11] embedded in coupled defects in photonic crystals [12,13] to coupled superconducting transmission line resonators [14] interacting with superconducting qubits [15,16,17]. In our study we have assumed small to medium size arrays (from 3 up to 16 sites) and values of the ratio of coupling to dissipation rate g/γ ∼ 20 which makes our results implementable with current designs in Circuit QED [18] and with near future set ups involving fiber coupled cavities [19] and photonic crystals [20].
The JCH as implemented in CRAs (with only the photonic excitations allowed to hop between neighboring cavities) motivated parallels with the predictions of the BH model. The latter naturally emerges in CRAs when one assumes generic nonlinear resonator effects instead of a Jaynes-Cummings interaction [21,22,23,24]. Indeed, the JCH and BH share many similarities, both describing bosons hopping coherently between nearest neighbor sites, with local nonlinearities. An equilibrium quantum phase transition between Mott-insulating and superfluid-like phases exists for both models. However, the BH Hamiltonian involves only a single species of bosons (photons in this case), while the excitations of the JCH model have both photonic and atomic components. Consequently, the equilibrium physics of the JCH model is expected to be richer as shown in [25,26,27,28,29]. Studying the JCH out of equilibrium, as naturally implemented in open driven CRAs, is likely to highlight even more interesting differences with novel features beyond the realm of the driven BH model. In addition, most existing work to date in out of equilbrium CRAs has used mean-field theory to treat the system [30,22,31,32,33,34]. Going beyond mean-field theory is crucial to provide physically accurate insights, especially in the experimentally realistic fewresonator regime. This requires a faithful representation of the full Liouville space of the resonator system which in most cases beyond two resonators becomes challenging.
For similar reasons of complexity, existing works to date on non-equilibrium resonator arrays exploring correlations and anti-bunching effects were always limited to minimal systems of two resonators [35,36,37,38,39]. More recent work includes a study of the fluorescence spectrum of again two coupled Jaynes-Cummings resonators [40,41]. The possibility of simulating gauge fields in driven dissipative Jaynes-Cummings arrays was investigated in [42], and artificial gauge fields in multi-resonator arrays in Bragg reflector micro-cavities assuming a Kerr interaction have also been studied recently [43].
In the present work we simulate CRAs beyond the two resonator regime, always assuming the full JCH model, by exploiting a combination of the matrix product state time-evolving block decimation (TEBD) algorithm [44,45] and quantum trajectories. We propose and analyze the many-body phases using experimentally accessible quantities such as the total excitation number, the emitted photon spectra and photon coherence functions for different, experimentally feasible, parameter regimes. In parallel, we also compare and contrast the expected behavior of this system when assuming that the local nonlinearity in the cavities is generated by a Kerr effect rather than a Jaynes-Cummings interaction. We find that the physics predicted by the two models differs for realistic regimes of interactions even if one meaningfully matches the respective strengths of the model nonlinearities. We analyze in detail the extra features of the JCH model originating from the mixed nature of the excitations and investigate the regimes where the Kerr description would faithfully match the JCH predictions. We find that the latter is only possible for values of the light-matter coupling and losses beyond the reach of current technology.
The structure of this manuscript is as follows. In the first section we describe the system and the dissipation/driving mechanisms. In Section 2 we discuss the local eigenstates and outline a scheme for mapping the effective nonlinearities between the models. In Section 3 we investigate the steady state spectra of single driven dissipative resonators, and demonstrate that Kerr physics can be achieved as a limiting case of the Jaynes-Cummings model. Section 4 treats arrays of a few sites as could be implemented in near future experimental setups, and results for yet larger arrays are presented in Section 5. In Section 6 we investigate the behavior of the system in the strongly repulsive regime and look for signatures of photon fermionization and crystallization in the JCH model. We then conclude in Section 7.

The system
We study M coupled single-mode resonators (indexed by j) in a circular configuration, as shown schematically in Fig. 1 (a) for M = 3. The resonator frequency ω c is the same for each resonator j. Coherent photon hopping between nearest-neighbor oscillators at a rate J arises as a result of an overlap between neighboring resonator modes j, j . A large lattice spacing, relative to optical lattice setups, means that external driving lasers with amplitude Ω j and frequency ω L can independently excite and probe each resonator j. Assuming each resonator is coherently interacting with a two-level system, Here two-level systems are embedded in cavities in a 2D photonic crystal but alternative implementations involving Circuit QED architectures could be also envisaged [18]. Coherent photon hopping between cavities results from modal overlap between neighboring cavities and photons can be lost due to decay mechanisms. (b) The relevant frequency scales and detunings used to describe a driven Jaynes-Cummings resonator. (c) The low-lying eigen-levels for the Jaynes-Cummings interaction, showing the definition of the effective Kerr nonlinearity U eff . (d) The low-lying excitations of a Kerr-nonlinear resonator. Arrows schematically show that a driving laser resonantly exciting a single-particle mode from the ground state is detuned from higher modes. Solid horizontal lines represent multiples of the cavity frequency. Dotted horizontal lines are a guide to the eye. the physics is well described by the JCH model as defined below. We will analyze the non-equilibrium system response and in parallel contrast the results with the case where a generic Kerr nonlinearity is assumed in place of the Jaynes-Cummings interaction as described by the well known BH model.
The local physics for both cases is captured by the two Hamiltonians ‡ ‡ Assuming the spacing between neighboring modes within a resonator to be much larger than all other scales we only employ a single photon mode per resonator.
where we have seth = 1. The difference between the resonator frequency, and the atomic transition frequency is denoted by ∆ and illustrated in Fig. 1 (b). The operatorâ † creates a photon, while the operatorsσ ± denote the usual raising and lowering operations between the ground and excited states of the two-level system.
The eigenstates ofĥ JC are known as 'dressed' states, or polaritons. They are mixed atom-photon excitations, which are also eigenstates of the total excitation number operatorN =σ +σ− +â †â with eigenvalue n. The ground state is |0 ≡ |g, 0 , with energy E 0 = 0. For a given excitation number n = 0, there are two polaritonic modes, designated |n, ± , with associated frequencies in the bare frame ω ± n = nω c − (∆/2) ± (∆/2) 2 + ng 2 . The eigenstates can be written in the bare resonator-atom basis as |n, ± = α ± n |e, n − 1 + β ± n |g, n , with α ± n = (χ n ∓ ∆)/( [46]. Here χ n = √ ∆ 2 + 4ng 2 . It can be seen from these coefficients that setting a larger positive ∆ results in an enhanced atomic component of the '-' polaritons, while negative ∆ increases the photonic contribution, with the relative composition of polaritons reversed for the '+' species. The increased splitting of the polaritonic frequencies from the bare resonator frequency with increasing n gives rise to an effective interaction between incoming photons, leading to a well known phenomenon in Cavity QED, the photon blockade [47]. Figure 1 (c) shows the eigen-structure schematically for the particular case ∆/g = 0, illustrating the energy mismatch between adjacent polariton manifolds. The resonance condition for an nphoton excitation of the |n, ± polaritonic mode in a single resonator from the ground state |0 is to set the driving laser detuning ∆ c = ω L − ω c = ω ± n /n − ω c at The local Hamiltonian for the driven BH model,ĥ BH , describes a single-mode resonator with a Kerr-type nonlinearity in the particle number of strength U . The eigen-frequencies are ω n = nω c + U 2 n(n − 1), with corresponding eigenstates being the Fock number states |n . The eigen-structure of the three lowest-lying levels is shown in Fig. 1 (d). An n-photon excitation of the n th mode occurs for the laser detuning After transforming to a frame rotating at the driving laser frequency ω L [48], the Hamiltonian for both model systems can be written generically aŝ The first term describes the local physics in resonator j in the rotating frame, with X ∈ {JCH, BH}. The second term describes the action of a coherent driving laser on each resonator j, and the last term describes the photon hopping. The local contributionsĥ BH are identical in form toĥ JC ,ĥ BH , with the bare cavity frequency replaced by the detuning ω c → −∆ c .
We now describe the manner in which we compare the two models. As the two Hamiltonians are patently different, care must be taken when comparing their behavior. We map the Jaynes-Cummings nonlinearity to an effective Kerr interaction via a frequency mismatch argument, as considered in both equilibrium and nonequilibrium contexts [29,49]. In the following we focus on the '-' species of polaritons, as the driving laser frequency necessary to resonantly excite the |n, − mode increases with n, qualitatively similar to a Kerr nonlinearity. We define the effective Kerr nonlinearity U eff = U eff (g, ∆) as the energy penalty incurred in forming a two-particle polaritonic excitation in a resonator (with energy ω − 2 ) from two one-particle polaritons in neighboring resonators (with total energy 2ω − 1 ) Interestingly, an analogous definition for the '+' polaritonic branch yields an attractive Kerr-type interaction. We note that as the distribution of energy levels with particle number n is markedly different between the models, it is only sensible to compare them in the very weakly excited regime where the above definition is meaningful. While the above mapping 'matches' the nonlinearities of the two models by considering transitions between the one-and two-particle manifolds, we expect the single particle resonances of both systems to occur at different spectral locations. The dominant spectral features in weakly excited systems will therefore be observed for different driving laser parameter regimes. Additionally, the mapping involves a comparison only of diagonal elements of the two governing Hamiltonians. However the different natures of the excitations of the models are not fully accounted for in such a mapping, in particular the additional internal degree of freedom possessed by the JCH model. We therefore expect the corresponding distinct off-diagonal Hamiltonian terms to lead to different physical observables between the models even under this mapping.
We assume a finite photon loss rate γ p from each resonator for both models, and that the spontaneous emission rate from the excited level |e of the two-level systems in the JCH system is negligible, γ a = 0. Competition between coherent resonator driving, and photon loss leads to NESS conditions. We employ a master equation formalism [50] to describe the evolution of the system's density matrix ρ(t). The NESS density matrix ρ ss is given by the stationary point of the master equationρ = L[ρ] = 0, where the action of the Liouville super-operator is defined through with [·, ·] + denoting the anti-commutator operation. In general, solving Eq. (7) is a formidable task, owing to the exponential growth in the necessary size of a system's description with the number of cavities M . We exploit the permutational symmetry of a homogeneous minimally sized three-site cyclic resonator system to enable solution via an exact diagonalization scheme, and employ a more sophisticated TEBD based stochastic unraveling of the master equation for larger systems. Details of the latter approach are provided in the Appendix. Before discussing the solutions of Eq. (7) for an array of cavities, we first revise and compare the stationary behavior of a single driven resonator with both Jaynes-Cummings and Kerr-type nonlinearities to illustrate their intrinsic differences in the local physics regime.

Local resonator physics: Kerr versus Jaynes-Cummings nonlinearities
In Figures 2 (a) and (b) we show NESS particle numbers for both types of resonator as a function of the driving laser detuning from the bare cavity resonance. The most striking difference between the spectra is the existence of two 'wings' for the Jaynes-Cummings nonlinearity, symmetrically distributed about the cavity frequency for this special case ∆/g = 0. The driving laser frequency necessary to excite the '-' polaritons increases with the excitation number n (the converse is true for the '+' polaritons). Figure 2 (b) also demonstrates that the Kerr resonator always exhibits a single particle response at the bare cavity frequency, whereas the spectral location of the two single particle modes in a Jaynes-Cummings type system strongly depend on the atom-resonator coupling parameters.
The nonlinear response of the resonators can also be probed by analyzing the second order correlation function of the emitted photon from a certain resonator site. For later use, we define the generalized coherence function between any two resonators (j, k) as g (2) In the following we refer to the on-site coherence function as g (2) ≡ g (2) (j, j). As expected, in the strong coupling regime this quantity exhibits a dip below the coherent driving laser value g (2) = 1 leading to photon antibunching when the driving laser is tuned to the corresponding single-particle modes for both models. For the Jaynes-Cummings resonator, the strongest anti-bunching is expected for laser detunings ∆ c = ±g from the bare cavity frequency (for ∆ = 0). In contrast, the Kerr resonator demonstrates strongest anti-bunching when the driving laser is on resonance with the bare cavity mode, coinciding with the single particle mode at ∆ c = 0. The reasonable atom-cavity coupling and cavity loss rates we have chosen give g/γ = 20 (shown in Fig. 2), causing the correlation function to dip to a minimum g (2) ≈ 0.25 for the Jaynes-Cummings resonator, and to g (2) ≈ 0.15 for the Kerr resonator.
Dips in the coherence function also appear when the driving laser is resonant with higher underlying quantum resonances (at equally spaced laser frequency intervals U/2 for the Kerr resonator, and laser frequencies ∆ c = ± g √ n for the Jaynes-Cummings system). The magnitude of the correlation function at these resonances is highly sensitive to the magnitude of the driving laser strength and cavity loss rate.
Figures 2 (a) and (b) share several features common to driven dissipative nonlinear quantum systems. We see that, for both models spectral peaks corresponding to higher excitations n become successively narrower, because the monochromatic driving laser  (4), respectively. Parameters: Both systems share Ω/γ p = 2. Jaynes-Cummings coupling g/γ p = 20. For the Kerr resonator: U/γ p ≈ 11.7, set using Eq. (6). (c) The response of a far-detuned (∆/g = −10) Jaynes-Cummings resonator for identical atom-resonator coupling, driving and loss as above (d) Comparison of the Kerr response in (b) and that of a detuned Jaynes-Cummings resonator (∆/g = −10) with an ultra-strong nonlinearity g /γ p ≈ 1.6 × 10 4 . cannot be simultaneously resonant with all intermediate levels. This means all but the lowest (polaritonic or Fock) peaks involve off-resonant transitions. The peaks also become less bright with increasing n because in addition to being more difficult to populate, modes with more photons are more susceptible to photon loss. The Jaynes-Cummings system involves both atomic and photonic species, leading to different peak intensities between the models. This follows because we are only driving the photonic degree of freedom which in turn is coupled to the atomic degree of freedom in the Jaynes-Cummings resonator.
Thus the spectral signatures of driven dissipative JCH and BH systems for the case of vanishing photon hopping (single resonators) are in general quite different. However, the two models do possess broad qualitative similarities in that they both describe the interplay of coherent bosonic hopping with an on-site nonlinearity. Indeed, the BH Hamiltonian has been widely used as an approximation to treat CRAs in several works assuming a generic nonlinearity [21,22,23,24]. It is natural to ask if the two models are equivalent in some regime or even if the physics of the JCH Hamiltonian can be mapped on to an effective BH model. Such a mapping would enable the large body of existing knowledge about the BH model to be applied directly to CRAs and would additionally significantly simplify numerical simulation. As the nonlinearity acts locally inside each resonator, we address this question below by first investigating single resonators.
The photonic limit of the Jaynes-Cummings model (∆ −g) -Intuitively, we expect agreement between the models as the polaritonic JCH excitations are made more photonic in nature. For the '-' polariton species, this means setting −g/∆ 1. The two-level system in each resonator is then barely excited and the dispersive interaction available is capable of inducing a significant effective photon repulsion only for very large values of the coupling g/γ p 1. In this limit a Taylor expansion of the polaritonic eigen-energies ω ± n in the small quantity |g/∆| leads to frequencies quadratic in the particle number n = 0 as ω n,− (g, ∆) ≈ n (ω c + g 2 /∆) + g(g/|∆|) 3 n(n − 1). We see that the off-resonant interaction induces an effective shift in the bare resonator frequency ∆ shift = −g 2 /|∆|, and that the energy spectrum can be written in the canonical Kerr nonlinear form with effective repulsion U approx BH ≡ 2g (g/|∆|) 3 . Figure 2 (c) shows the spectrum for a far detuned (∆/g = −10) driven dissipative Jaynes-Cummings resonator for the same experimentally realistic atom-cavity coupling rate g/γ p as in Fig. 2 (a). We see that the response is essentially that of an empty (linear) driven dissipative resonator, albeit with a peak response at ∆ shift from the bare resonator frequency. The atomic excitation is an order of magnitude smaller, as expected. The effective Kerr nonlinearity is much smaller than the line-width (U approx BH /γ p = 1 50 ), so that nonlinear effects barely show up, apart from a slight asymmetry in the response. To observe sizable nonlinear effects for realistically achievable resonator parameters, we see that it is instead necessary to operate near the resonance point ∆/g ≈ 0.
For comparison, the physically unrealistic ultra-strong atom-resonator coupling limit is shown in Fig. 2 (d). The response of the Kerr-nonlinear resonator shown in Fig. 2 (b) is reproduced, and compared with the spectrum of a Jaynes-Cummings resonator again operating in the photonic regime (∆/g = −10). The coupling strength g/γ p > 10 4 is chosen by setting U approx BH = U in the above definition of the effective Kerr nonlinearity. We see good agreement between the photon number spectra, which worsens for higher excitation peaks as higher order terms in the Taylor expansion of the polaritonic eigen-frequencies become important. So while deep in the regime −g/∆ 1 the Jaynes-Cummings model can, with reasonable accuracy, be mapped on to a Kerr system, this regime is inaccessible with current technology.
The atomic limit of the Jaynes-Cummings model (∆ g) -In the opposite limit of ∆/g 1, the atomic component of the '-' polaritons is maximal. A Taylor expansion of the eigen-frequencies now yields: ω n,− (g, ∆) ≈ n (ω c − g 2 /∆) + g(g/∆) 3 n(n − 1) − ∆, for n = 0. The ground state |g, 0 still lies at the zero of energy. Therefore we see that to first order in g/∆, the spectrum becomes equally spaced with the same shift in the resonator frequency as for the limit ∆ −g, with the exception of the interval ω c − ∆ between the ground state and the first polaritonic mode |1, − . In this limit |1, − ≈ |e, 0 , with a small photonic component allowing transitions by the external driving laser. This may be summarized by the following effective Hamiltonian describing the 'vacuum shifted' |n, − ladder of states for a single resonator in the atomic limit: Thus, population of the two-particle state from the resonantly driven single-particle mode is strongly inhibited. Higher lying levels could, however, be near resonantly populated from the one-particle mode given an additional driving laser with frequency ω L = ω c . When strictly working within the low-excitation regime, it is then appropriate to assign a large effective nonlinearity describing this energy penalty to reach the twoexcitation manifold. However, care must be taken to distinguish between setups such as we consider, where only the lowest lying excitations are directly probed, and others which access the approximately harmonic ladder at larger n > 2.

Many-body signatures in steady state observables
Moving beyond the single resonator regime we will first analyze a minimal cyclic nonlinear CRA of M = 3 cavities. This case, though not truly many-body is interesting as this is where the first experimental implementations are likely to begin [18,19,20]. We consider a finite coherent photon tunneling J = 0 between adjacent cavities which splits the local polaritonic resonances into delocalized global modes. For the moment we drive our system homogeneously so that all external lasers are in-phase with Ω j = Ω, ∀j.
We obtain the NESS by diagonalizing the super-operator L after exploiting the permutational symmetry of the system to significantly decrease the number of unique density matrix elements and additionally retaingin only basis states allowing a maximum of P = 4 excitations in the system. § Figures 3 (a) and (b) show generalizations of the single resonator spectra discussed above for the relatively small rate of photon hopping J/γ p = 1. We expect the response of the systems to strongly resemble the 'local' physics in this regime, as the photon hopping is essentially a weak perturbation on top of the atom-cavity coupling.
The homogeneous nature of our chosen driving means that the total momentum of excitations in the NESS must be zero. The one-particle excitation in the BH system Figure 3. The steady state per-resonator photonic/atomic populations and g (2) function for homogeneously driven, cyclic three-site resonator systems described by (a) the JCH Hamiltonian (with ∆/γ p = 2J/γ p = 2) and (b) the BH model. Vertical dashed lines indicate the location of underlying modes of an isolated resonator (J = 0), while vertical dash-dotted lines indicate the positions of delocalized single particle modes. The atom-resonator coupling g/γ p = 20 for (a) and (b) is again as in Figs. 2, but there is now a finite photon hopping J/γ p = 1. (c) and (d) are generalizations of the single resonator 'photonic limit' spectra of Figs. 2, extended to three-resonator systems with the same photon hopping. The driving rate in (c) is set at Ω/γ p = 0.3 to enable numerical solution, otherwise all parameters are unchanged. Note in (d) the almost complete overlap of the response of the two many-body models now in this 'photonic' limit. This agreement, however, requires an unrealistic ultra-strong value of g/γ p ≈ 1.6 × 10 4 .
is then just the zero-momentum free-particle Bloch mode |k = 0 = 1 √ M M j=1â † j |0 excited at a laser detuning (∆ c ) |k=0> = −2J, as the Kerr nonlinearity does not influence a single particle. In contrast a single excitation in the JCH system can be shared between atomic and photonic degrees of freedom leading to two delocalized generalizations of the |1, ± polaritons, denoted |k = 0, ± = A ± |E + B ± |k=0 . Here |E = 1 √ M M j=1σ + j |0 is a delocalized atomic excitation. The coefficients A ± and B ± are identical in form to the coefficients α ± 1 , β ± 1 of the localized polaritons, after making the replacement ∆ → ∆ k=0 = ∆ − 2J. These are resonantly excited at laser detunings As only the cavities are directly coupled there is now an additional asymmetry between photons and atoms. The resonance point for delocalized polaritons is shifted to ∆ = 2J from ∆ = 0 found in the single resonator case. This is equivalent to setting the transition frequency of the atoms equal to the frequency of the lowest lying Bloch mode |k=0 .
Relative to Figs. 2 (a) and (b), we see new multi-particle modes appear between the delocalized generalizations of the dressed state resonances. By an N -particle mode we mean an eigenstate of the total excitation number operator for the whole system N = M j=1N j with eigenvalue N . For this particularly small photon hopping rate, equal to the cavity line-width, the delocalized single particle mode of the BH system is smeared into a broad hump along with delocalized two-and three-particle modes, while new features can also be discerned in the JCH spectra. Additionally, new modes appear in the JCH model at approximately the bare cavity frequency. These modes are symmetrized superpositions of '+' and '-' polaritons, and for our chosen driving strength are barely populated relative to the response of the BH model at the bare cavity frequency.
There is an asymmetry between the delocalized generalizations of the '-' and '+' wings of the JCH resonator array spectrum, despite setting ∆ = 2J. We note that for large photon hopping J g and ∆ = 2J, the resonator array spectrum is again symmetric about the Bloch mode frequency (not shown), though still with quantitative differences to the single resonator case. It is in the intermediate regime (J ≈ g) that asymmetries appear, as hopping brings multiple particle global excitations into the same spectral region. A signature of this hopping-induced behavior that is particularly amenable to experimental verification is a measurement of the photon correlation at the two delocalized single-particle resonance frequencies. For this particular driving and cavity loss rate, the values are g (2) ≈ 0.18, 0.35 for the |k = 0, ± modes respectively. Meanwhile g (2) reaches a minimum of ≈ 0.26 at the underlying single-particle mode of the BH system (indicated by the vertical dash-dotted line in Fig. 3 (b)). We note here again that if one assumes losses smaller than g/γ p = 20 used here then the correlation minimum approaches zero and clear anti-bunching should be achieved, as expected.
The spectral response of far-detuned finite-size CRAs operating in the photonic limit (−g/∆ 1) is shown in Figs. 3 (c) and (d). These spectra are direct generalizations of those shown in Fig. 2. Again, we see that for our assumed atomresonator couplings (i.e. g/γ p ≈ 20), the main effect of the atomic degree of freedom is to shift the free-particle mode by ∆ shift , in this case from the lowest Bloch mode. Figs. 3 (d) shows that (unrealistically) larger couplings 'matched' to a specific Kerr nonlinearity can reproduce the spectral features of a driven dissipative finite CRA with reasonable accuracy. . NESS particle numbers per resonator for sixteen-site Bose-Hubbard and Jaynes-Cummings CRAs driven at their single particle resonances. Parameters: photonic driving Ω/γ p = 2, atom-resonator coupling strength g/γ p = 20, photon hopping rate J/γ p = 1. Trajectory calculations retain p = 6 photons per resonator for the BH simulation, and p = 5 for the JCH. Each Jaynes-Cummings detuning ∆/g is mapped to an effective Kerr nonlinearity via Eq. (6) to enable a comparison.

Larger resonator arrays
So far we have investigated few sites JCH systems of a few sites operating either onresonance (i.e. ∆/g = 0), or in the photonic limit where the weak nonlinearity is captured by an effective BH interaction. We now demonstrate the effect of the changing nature of the polaritonic excitations and also investigate larger arrays of M = 16 sites as this parameter is varied about the strong-interaction regime ∆/g ≈ 0. Rather than construct a spectrum by varying the driving laser frequency, we instead selectively drive a particular spectral feature and vary ∆/g. Specifically, in Fig. 4 we show trajectory results for NESS particle numbers for larger (M = 16) JCH and BH arrays driven at their single particle resonances, corresponding to the driving laser detunings ∆ c = (∆ c ) |k=0,− and ∆ c = −2J respectively.
As expected, we see the changing nature of the system's excitations reflected in a larger atomic (and smaller photonic) occupation for increasing ∆/g. For strong positive detunings, the excitations are predominantly atomic, indirectly excited via the resonator field. As we have assumed lossless atoms, the steady state corresponds to oscillations between the ground and excited states, leading to an average half atomic occupancy. Under the mapping of Eq. (6), this regime corresponds to a BH system with large Kerr coefficient. Under the stronger driving we have chosen for this calculation, the BH system approximately oscillates coherently between zero and one photons per resonator, with occupation of higher Fock levels suppressed. Thus, while the total excitation number for the JCH asymptotically agrees with the NESS photon number for the BH in this limit, the underlying physics is very different.
In the opposite limit ∆/g < 0, the excitations of the JCH model become more photonic, and the effective Kerr nonlinearity decreases. We have already seen that agreement between the JCH and Kerr photon number is reached far into this limit.
In the intermediate regime, ∆/g ≈ 0, the total particle number for the JCH exhibits significant departures from the BH results. This is due both to the fact that only one component of the polaritons is driven in the JCH system, and also that this stronger driving allows access to higher-lying states (n > 2) outside the regime of validity for the effective Kerr strength definition. We therefore see non-trivial differences between the models both in spectra measured as a function of driving laser frequency, and as a function of the nonlinearity.

The strong nonlinearity limit: fermionization and crystallization of photons in CRAs
Having demonstrated that the NESS in resonator arrays governed by the JCH and BH Hamiltonians are in general different we now evaluate two non-equilibrium effects in the large nonlinearity regime recently studied assuming a BH description, with proposals to experimentally realise such effects in Jaynes-Cummings type CRA systems. Fermionization in coupled resonator arrays -Recently an exploration of the ultra-strong nonlinearity regime in a three-site driven dissipative Bose-Hubbard model was undertaken [21]. For a closed coherent system it was found that as the Kerr strength U approached infinity double occupancy of any resonator was completely suppressed allowing the system's bosonic wave-functions to be mapped to those of an equivalent fermionic system via a Jordan-Wigner transformation. In this interpretation double occupancy is prevented by fermionic statistics rather than a hard core bosonic interaction. The authors found that a cyclic system's N -particle eigenfunctions and corresponding energies can then be identified uniquely for a given total particle momentum. Subjecting the system to coherent driving and a finite particle loss rate from each resonator results in readily classifiable peaks in the resonator occupancy as a function of the driving laser frequency. Figure. 5 (a) reproduces such a spectrum for a homogeneously driven (Ω j = Ω, ∀j) minimal M = 3 resonator BH system. For the spectral range shown the only possible modes that can be observed in the hard core limit NESS are the zero momentum oneparticle mode, which for the BH model coincides with the Bloch mode |k = 0 of an empty resonator system, and a mode formed from two particles with opposite momenta k = ±2π/3. We show spectra for a strong nonlinearity U/γ p 1 (allowing for any number of photons per resonator in the calculation) and for the true fermionized limit U → ∞ (when strictly p = 1 photon per resonator is retained in calculations) illustrating the convergence. Also shown in Fig. 5 (a) are NESS photon numbers for a strongly detuned JCH system operating in the 'photonic regime', again computationally retaining multiple and single photons per resonator. As expected from the discussion in Sec. 3 good agreement with the BH results is observed, albeit for unrealistically large values of the light-matter coupling g. Figure 5. (a) The one and two particle particle peaks (right and left features, respectively) for a homogeneously driven three-resonator BH model in the strongly nonlinear regime (solid blue, U/J = 23.4), and the hard core limit (p = 1, dotted black). Spectra for larger U approach the hard core limit asymptotically. Parameters: photon hopping rate J/γ p = 20, driving Ω/γ = 0.5. Also shown are results for a strongly detuned JCH system operating in the photonic limit with ∆/g = −10, and an unrealistically large atom-resonator coupling g/γ p ≈ 2 × 10 5 . (b) Analogous photon number spectra for the JCH model with g/γ p ≈ 800 (chosen according to Eq. 6), and the photon driving and tunneling rates as in (a), for two different atom-resonator detunings ∆/g (dashed and dash-dotted). Again, the hard-core spectra are recovered in the limit g → ∞. The horizontal axis is taken relative to the single particle spectral position (∆ c ) |k=0,− , different for each ∆. (c) The auto-and cross-photon correlation functions measured at the two particle peak as a function of the nonlinearity in both the BH and JCH (∆ = 2J) Hamiltonians. JCH systems with different couplings g are compared with Bose-Hubbard systems whose strength U = U (g) (see upper horizontal axis) is obtained from Eq. 6. Note these calculations are performed at the lower driving Ω/γ p = 0.25, to avoid truncation errors at low nonlinearity.
Moving beyond from the photonic limit we now explore the novel physics of ultrastrong nonlinearity in JCH systems operating in the regime of strongest interaction between resonators and atoms, i.e. ∆/g ≈ 0. In this regime correspondence between the models is less clear and the fully polaritonic nature of excitations must be taken into account. To recover physics which resembles hard-core photons we operate in the limit g/γ p 1, along with g ω c to satisfy the rotating wave approximation, and drive at ω L ≈ ω − 1 near the |1, − polaritonic resonance of isolated resonators. Under these conditions the only relevant states in each resonator are the ground state and |1, − polaritonic state. Figure 5 (b) shows hard-core JCH model spectra for two different atom-resonator detunings ∆ = 0 and ∆ = g. Increasing the coupling g while holding the ratio ∆/g fixed leads to a qualitatively similar evolution of spectral features as in the BH model. A two-particle excitation splits from the single particle resonance and asymptotically approaches a splitting that only depends on J and ∆/g. For the coupling g/γ p ∼ 800 used good agreement is seen in Fig. 5 (b) between the full JCH model and the spectrum obtained when retaining only the p = 1 photons per resonator necessary to describe the hard-core polariton limit. We see that even in the ultra-strong non-linearity limit of the JCH model the internal structure of polaritons (governed by the parameter ∆/g) still plays a crucial role in determining the location and amplitude of spectral features. This behavior is a result of the effective hopping, driving and loss rates for polaritons being different from the bare photon parameters. Indeed the effective polaritonic hopping rate is J pol = |β − 1 | 2 J, the effective driving rate is Ω pol = (β − 1 )Ω, while the loss rate is γ pol = |β − 1 | 2 γ p , where β − 1 is the coefficient controlling the photonic component of the |1, − polaritons defined in Sec. 2. This reflects that the hopping, driving and loss processes involve only the photonic component of the polaritons. Thus, on resonance (∆/g = 0 ⇒ β − 1 = 1 √ 2 ) we see that the peak separation is J pol = J/2 as seen in Fig. 5 (b). Focusing on experimentally measurable photonic quantities, in Fig. 5 (c) we track for both models as a function of their nonlinearity the photon density-density correlations measured at the two-particle peak on a single site via g (2) (j, j) and between neighboring sites via g (2) (j, j + 1). The spectral location of the zero-momentum twoparticle modes for the JCH and BH systems, are found using the results in Refs. [51] and [52] respectively. As outlined in Ref. [21] at small nonlinearities the two-particle peak resides within the one-particle spectral feature and so correlations inherit Poissonian statistics from the driving laser. Larger nonlinearities split the two-particle resonance from the one-particle peak, as shown in Figs. 5 (a) and (b), and lead to strong bunching. For very strong nonlinearities on-site anti-bunching is expected as the two excitations are distributed in such a way that the pair of photons are never in the same resonator resulting in the nearest neighbor correlations becoming large. Such considerations, being consequences of generic nonlinear behavior, lead to qualitatively similar correlation functions for both the driven dissipative JCH and BH models. Thus fermionized photons are a feature of the on-resonance JCH model as well, once the different spectral frequencies for correlation measurements are taken into account.
Polariton crystallization -Another intriguing phenomenon of interacting photons in resonator systems, photon crystallization, was recently predicted to occur in a one dimensional ring of optical cavities with Kerr-type nonlinearity [23]. Driving lasers with a phase difference of π/2 between each site k, i.e. Ω k = Ω exp(ikπ/2), create a flow of bosons around the system. The contact interaction energy U was found to result in a 'crystallization' of particles as they flowed around the ring, even in the presence of dissipation. The signature of this effect was identified in the particle density-density correlations g (2) (j, k) between cavities j and k, measured in the system's steady state.
On-site anti-bunching is accompanied by nearest-neighbor density-density correlations stronger than correlations between more distant cavities. The conclusion drawn in Ref. [23] was that particles form 'dimers' of light, which flow around the system. We demonstrate here that this dimerization can also be seen in a JCH system outside the photonic limit. The additional degree of freedom in the atom-resonator detuning ∆ allows the strength of this effect to be adjusted on demand.
We study a system of M = 16 coupled Jaynes-Cummings resonators under periodic boundary conditions. In addition to the relative phases, we must choose the frequency of our driving lasers. For a BH-type system, a laser detuning ∆ c = 0 directly drives the single particle k = π/2 mode. For the JCH model we choose to drive the delocalized generalization of the single-particle polaritonic mode |1, − k=π/2 to produce an analogous effect. Details of the parameters used for this numerically demanding calculation are given in the Appendix.
We see in Fig. 6 (a) that the signatures of crystallization show up in our polaritonic system -there is a larger probability of finding photons in neighboring cavities than in cavities further apart. We have additional control over the atom-resonator detuning ∆, which allows us to tune the nature of the system's excitations either more photonic ∆/g < 0, or atomic ∆/g > 0. We see that the crystallization effect becomes weaker for larger ∆ as the atoms and resonators are tuned away from resonance, while on-site anti-bunching is enhanced. This is a consequence of resonator photons being strongly coupled to the two-level systems, which hold most of the excitation for large ∆, but longer-range correlations tend to unity. Figure 6 (b) shows the actual photonic and atomic expectation values per resonator, for the three detunings considered illustrating the changing make-up and associated enhanced photon anti-bunching. Also shown for comparison in Figure 6 are results for the driven dissipative Kerr system like that considered in Ref. [23]. We conclude that the 'dimerization' predicted there for the BH model persists in CRA calculations using the full atom-resonator Hamiltonian without approximations, and that some degree of tunability should be observable by bringing the cavities in and out of resonance while maintaining significant resonator populations.

Discussion
We have proposed and analyzed excitation number and photon coherence spectra for small to medium sized weakly driven dissipative coupled resonator systems described by the Jaynes-Cummings-Hubbard Hamiltonian. Simulations were at all times performed using the full JCH Hamiltonian without approximations, using model parameters realizable by current state-of-the-art technology. We have also presented analogous spectra for resonator arrays governed by the Bose Hubbard model, drawing attention to Figure 6. (a) Steady state photon density-density correlations for a cyclic 16 site JCH system, driven by lasers exciting the π/2 momentum mode of the system. Thin lines between markers are drawn to guide the eye. Parameters: hopping J/γ p = 2, atom-resonator coupling g/γ p = 10, driving strength |Ω k |/γ p = 2. Also shown (solid green) are density-density correlations for a M = 16 site system with 'matched' Kerr nonlinearity U/γ p ≈ 6 in each resonator, and otherwise identical parameters. (b) Closer view of the correlation functions. (c) The relative atomic and photonic population in each resonator in the steady state for the three different values of atom-resonator detuning chosen, as well as the on-site photon correlation function. the differences in experimentally accessible observables between the two models. These differences arise primarily because of the composite nature of the elementary excitations of the JCH Hamiltonian, leading to readily identifiable unique spectral characteristics. Therefore we conclude that when simulating coupled resonator arrays the full JCH Hamiltonian must be used in calculations to properly account for the underlying physics.
Generalizations of two bosonic interaction-induced non-equilibrium phenomena were observed in the NESS of CRAs modeled retaining the richer JCH physics, with additional tunability. The fingerprints of the polaritonic equivalents of 'fermionic' photons and photon crystallization exhibit subtle departures from the pure bosonic case, detectable via experimentally accessible spectral quantities.
Finally, this work has shown that the combination of the quantum trajectory method with the TEBD algorithm is an especially powerful tool in finding the NESS of open driven dissipative many-body systems. Such a numerical approach can be readily applied to calculate more complex out of equilibrium properties of coupled resonator arrays where an analytic approach is unfeasible. This might include important open problems such as determining the transport properties of linear coupled resonator arrays in the presence of photon loss, or exploring novel one-dimensional quantum states of light that are robust to experimentally realistic decay processes. improvement in accuracy afforded by the time averaging procedure, a relatively modest matrix dimension in the MPS system description is sufficient for accurate results. This is a consequence of long range correlations in the system being constantly broken up by the local incoherent processes. A quantitative analysis of the relationship between correlations in individual trajectory wave-functions, and those in the steady state density matrix, will be presented elsewhere [59].
Polariton crystallization calculations -To obtain the steady states used in Fig. 6, we began time averaging at g(T start ) = 500. We found that N T = 5000 time-steps g(∆t) = 0.5, with around 100 trajectories retaining around 40 states in the matrix product state representation, were sufficient to obtain acceptable statistical fluctuations (≈ 0.1%) in the NESS expectation values we are interested in. We found that retaining three or four photons in the bare resonator basis was sufficient for systems in the 'atomic' regime (∆ ≥ 0), but that more were required to properly evaluate steady state quantities for ∆ < 0 (results not shown).