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

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 addressability 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–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–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–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–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 equilibrium CRAs has used mean-field theory to treat the system [22, 30–33]. Going beyond mean-field theory is crucial to provide physically accurate insights, especially in the experimentally realistic few-resonator 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 [34–38]. More recent work includes a study of the fluorescence spectrum of again two coupled Jaynes–Cummings resonators [39, 40]. The possibility of simulating gauge fields in driven dissipative Jaynes–Cummings arrays was investigated in [41], and artificial gauge fields in multi-resonator arrays in Bragg reflector micro-cavities assuming a Kerr interaction have also been studied recently [42].

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 (MPS) time-evolving block decimation (TEBD) algorithm [43, 44] 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 (modeled using the BH Hamiltonian) 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 section 1, 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.

We study M coupled single-mode resonators (indexed by j) in a circular configuration, as shown schematically in figure 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, 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.

Zoom In Zoom Out Reset image size

The local physics for both cases is captured by the two Hamiltonians⁵

$$\begin{eqnarray} \hat{h}'_{\mathrm{JCH}} & = & \omega_{\mathrm{c}} \hat{a}^{\dagger} \hat{a} + \left (\omega_{\mathrm{c}} - \Delta \right) \hat{\sigma}^+ \hat{\sigma}^- + g \left (\hat{a}^{\dagger} \hat{\sigma}^- + \hat{a} \hat{\sigma}^+ \right), \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \hat{h}'_{\mathrm{BH}} & = & \omega_{\mathrm{c}} \hat{a}^{\dagger} \hat{a} + \frac{U}{2} \hat{a}^{\dagger} \hat{a}^{\dagger} \hat{a} \hat{a}, \end{eqnarray} \tag{ 2 }$$

where we have set ℏ = 1. The difference between the resonator frequency, and the atomic transition frequency is denoted by Δ and illustrated in figure 1(b). The operator $\hat {a}^{\dagger}$ creates a photon, while the operators $\hat {\sigma }^\pm$ denote the usual raising and lowering operations between the ground and excited states of the two-level system.

The eigenstates of $\hat {h}'_{\mathrm {JCH}}$ are known as 'dressed' states, or polaritons. They are mixed atom–photon excitations, which are also eigenstates of the total excitation number operator $\hat {\mathcal {N}} = \hat {\sigma }^+ \hat {\sigma }^- + \hat {a}^{\dagger} \hat {a}$ with eigenvalue n. The ground state is |0〉 ≡ |g,0〉, with energy E₀ = 0. For a given excitation number n ≠ 0, there are two polaritonic modes, designated |n,±〉, with associated frequencies in the bare frame $\omega _n^\pm = n \omega _{\mathrm {c}} - (\Delta / 2) \pm \sqrt {\left (\Delta / 2 \right)^2 + n g^2}$. The eigenstates can be written in the bare resonator-atom basis as |n,±〉 = α^±_n|e,n − 1〉 + β^±_n|g,n〉, with $\alpha _n^\pm = (\chi _n \mp \Delta) / (\sqrt {2}\sqrt {\chi _n^2 \mp \Delta \chi _n})$ and $\beta _n^\pm = \pm 2\,g \sqrt {n}/(\sqrt {2}\sqrt {\chi _n^2 \mp \Delta \chi _n})$ [45]. Here $\chi _n = \sqrt {\Delta ^2 + 4 n g^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 [46]. 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 n-photon 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

$$\begin{equation} \left (\Delta_{\mathrm{c}} \right)_n^\pm = \frac{1}{2 n} \left (\Delta \pm \chi_n \right). \end{equation} \tag{ 3 }$$

The local Hamiltonian for the driven BH model, $\hat {h}_{\mathrm { BH}}'$, describes a single-mode resonator with a Kerr-type nonlinearity in the particle number of strength U. The eigen-frequencies are $\omega _n = n \omega _{\mathrm {c}} + \frac {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 figure 1(d). An n-photon excitation of the nth mode occurs for the laser detuning Δ_c = ω_L − ω_c = (ω_n/n) − ω_c

$$\begin{equation} (\Delta_{\mathrm{c}})_n = \frac{U}{2}(n-1). \end{equation} \tag{ 4 }$$

After transforming to a frame rotating at the driving laser frequency ω_L [47], the Hamiltonian for both model systems can be written generically as

$$\begin{equation} \hat{H}_X = \sum_{j=1}^M \hat{h}_X^{(j)} + \sum_{j=1}^M \Omega_j \left (\hat{a}_j^{\dagger} + \hat{a}_j \right) - J \sum_{\langle j,j'\rangle} \hat{a}_j^{\dagger} \hat{a}_{j'}. \end{equation} \tag{ 5 }$$

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 $\hat {h}_{\mathrm {JCH}}^{(j)},\hat {h}_{\mathrm {BH}}^{(j)}$ are identical in form to $\hat {h}_{\mathrm {JCH}}',\hat {h}_{\mathrm {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. A one-to-one mapping between models is not possible for arbitrary Hamiltonian parameters, and we must therefore 'match' some specific features of the models. In this work we probe only the lowest-lying excitations of our lattice systems. Much of the physics in this 'dilute' regime follows from the effective repulsion between just two quasi-particles per resonator, which in turn follows from the energy cost incurred when bringing these two excitations together. We therefore map the Jaynes–Cummings nonlinearity to an effective Kerr interaction via a frequency mismatch argument, as considered in both equilibrium and non-equilibrium contexts [29, 48, 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 ω⁻₂) from two one-particle polaritons in neighboring resonators (with total energy 2ω⁻₁)

$$\begin{equation} \frac{U_{\rm eff}}{g} = \frac{\omega_2^- - 2 \omega_1^-}{g} = \frac{\Delta}{2g} + 2 \sqrt{\left (\frac{\Delta}{2g} \right)^2 + 1} - \sqrt{\left (\frac{\Delta}{2g} \right)^2 + 2} . \end{equation} \tag{ 6 }$$

Figure 1(c) shows these low-lying polaritonic eigenstates and the definition of U_eff pictorially, while figure 1(d) shows the level structure for a Kerr-nonlinear resonator. Interestingly, an analogous definition for the ' + ' polaritonic branch yields an attractive Kerr-type interaction. 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. We note that several early works on the equilibrium physics of the JCH discuss the different nonlinear behavior expected of Jaynes–Cummings and Kerr-interactions, specifically the fact that the Jaynes–Cummings nonlinearity saturates with increasing excitation number while the Kerr nonlinearity does not [5, 26, 29, 50].

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 [51] 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 $\dot {\rho } = \mathcal {L}[\rho ] = 0$, where the action of the Liouville super-operator is defined through

$$\begin{equation} \mathcal{L}_X [\rho_{\rm ss}] = \frac{1}{\mathrm{i}} [\hat{H}_X, \rho_{\rm ss}] + \sum_{j=1}^M \frac{\gamma_{\mathrm{p}}}{2} \left (2 \hat{a}_j \rho_{\rm ss} \hat{a}^{\dagger}_j - [\hat{a}^{\dagger}_j \hat{a}_j, \rho_{\rm ss}]_+ \right) = 0, \end{equation} \tag{ 7 }$$

with [·,·]₊ denoting the anti-commutator operation. In general, solving equation (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 equation (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.

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 depends on the atom–resonator coupling parameters.

Zoom In Zoom Out Reset image size

The nonlinear response of the resonators can also be probed by analyzing the second order correlation function of the emitted photons from a certain resonator site. For later use, we define the generalized coherence function between any two resonators (j,k) as g⁽²⁾(j,k) = 〈a^†_ja^†_ka_ja_k〉/〈a^†_ja_j〉〈a^†_ka_k〉. In the following we refer to the on-site coherence function as g⁽²⁾ ≡ g⁽²⁾(j,j). As expected, in the strong coupling regime this quantity exhibits a dip below the coherent driving laser value g⁽²⁾ = 1 leading to photon anti-bunching 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 figure 2), causing the correlation function to dip to a minimum g⁽²⁾ ≈ 0.25 for the Jaynes–Cummings resonator, and to g⁽²⁾ ≈ 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 $\Delta _{\mathrm {c}} = \pm \frac {g}{\sqrt {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 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–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.

3.1. The photonic limit of the Jaynes–Cummings model (Δ ≪ − g)

3.2. 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: $\omega _{n,-}(g,\Delta) \approx n \left (\omega _{\mathrm {c}} - g^2 / \Delta \right) + g (g / \Delta)^3 n (n-1) - \Delta$, 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:

$$\begin{equation} \hat{h}_{\Delta \gg g} = (\omega_{\mathrm{c}} - g^2 / \Delta) \hat{a}^{\dagger} \hat{a} + g \left (\frac{g}{\Delta} \right)^3 \hat{a}^{\dagger} \hat{a}^{\dagger} \hat{a} \hat{a} + \Delta |g,0 \rangle \langle g,0| - \Delta. \end{equation} \tag{ 8 }$$

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 two-excitation 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.

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–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 $\mathcal {L}$ after exploiting the permutational symmetry of the system to significantly decrease the number of unique density matrix elements and additionally retaining 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 system 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.

Zoom In Zoom Out Reset image size

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 is then just the zero-momentum free-particle Bloch mode $|k = 0 \rangle = \frac {1}{\sqrt {M}} \sum _{j=1}^M \hat {a}^{\dagger} _j |0 \rangle$ 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 \rangle = \frac {1}{\sqrt {M}}\sum _{j=1}^M \hat {\sigma }^+_j |0 \rangle$ is a delocalized atomic excitation. The coefficients A^± and B^± are identical in form to the coefficients α^±₁, β^±₁ of the localized polaritons, after making the replacement Δ → Δ_k=0 = Δ − 2J. These are resonantly excited at laser detunings (Δ_c)_|k=0,±〉 = −2J − Δ_k=0/2 ± χ⁽¹⁾_k=0, where $\chi _{k=0}^{(1)} = \sqrt {g^2 + \Delta _{k=0}^2/4}$. 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 figures 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 $\hat {\mathcal {N}} = \sum _{j=1}^M \hat {\mathcal {N}}_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⁽²⁾ ≈ 0.18,0.35 for the |k = 0,±〉 modes respectively. Meanwhile g⁽²⁾ reaches a minimum of ≈0.26 at the underlying single-particle mode of the BH system (indicated by the vertical dash-dotted line in figure 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 figures 3(c) and (d). These spectra are direct generalizations of those shown in figure 2. Again, we see that for our assumed atom–resonator 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. Figure 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.

So far we have investigated JCH systems of a few sites operating either on-resonance (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 figure 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.

Zoom In Zoom Out Reset image size

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 equation (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 photon 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.

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 realize such effects in Jaynes–Cummings type CRA systems.

6.1. Fermionization in coupled resonator arrays

Recently an exploration of the ultra-strong nonlinearity regime in a three-site driven dissipative BH 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 one-particle 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 figure 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 section 3 good agreement with the BH results is observed, albeit for unrealistically large values of the light–matter coupling g.

Zoom In Zoom Out Reset image size

Moving beyond the photonic limit we now explore the novel physics of ultra-strong 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 ≈ ω⁻₁ 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 figure 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 nonlinearity 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 = |β⁻₁|²J, the effective driving rate is Ω_pol = (β⁻₁)Ω, while the loss rate is γ_pol = |β⁻₁|²γ_p, where β⁻₁ is the coefficient controlling the photonic component of the |1,−〉 polaritons defined in section 2. This reflects that the hopping, driving and loss processes involve only the photonic component of the polaritons. Thus, on resonance ($\Delta / g = 0 \Rightarrow \beta _1^- = \frac {1}{\sqrt {2}}$) we see that the peak separation is J_pol = J/2 as seen in figure 5(b).

Focusing on experimentally measurable photonic quantities, in figure 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⁽²⁾(j,j) and between neighboring sites via g⁽²⁾(j,j + 1). The spectral location of the zero-momentum two-particle modes for the JCH and BH systems are found using the results in [52, 53] respectively. We set a weaker driving Ω/γ = 0.25 than that in figures 5(a) and (b) so as to avoid truncation error at low nonlinearities, where the steady state exhibits a long-tailed Poissonian number distribution characteristic of coherent states. As outlined in [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 figures 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.

6.2. 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. $\Omega _k = \Omega \exp (\mathrm {i} k \pi / 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⁽²⁾(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 [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 figure 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(c) 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 [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.

Zoom In Zoom Out Reset image size

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 JCH Hamiltonian. Simulations were at all times performed using the full JCH Hamiltonian and fully accounted for losses. Current state-of-the-art setups involving transmission-line resonators coupled to superconducting qubits can access a wide range of Hamiltonian coupling rates and photon leakage rates, by varying the geometry defined by lithographic patterning. The qubit–resonator coupling rate g can be set in the range 1–300 MHz, while the resonator tunneling rate J and loss rate γ_p can be independently varied in the range 50 kHz–50 MHz, so that our simulations model physically realizable systems. We have also presented analogous spectra for resonator arrays governed by the BH model, drawing attention to 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.

Here we give a short summary of the numerical components involved in our NESS calculations. We describe the quantum trajectory method, the ansatz used for the stochastic wave-functions, and the algorithm to evolve the system description.

A.1. The quantum trajectory algorithm

A.2. The matrix product state representation

Time evolving the stochastic wave-functions is itself a non-trivial task, as the dimension of Hilbert space is still exponentially large for a representation in the bare basis, as seen by writing the state of an M-site lattice generically as

$$\begin{equation} | \psi \rangle=\sum_{j_1=1}^{\mathbf{d}}\cdots\sum_{j_N=1}^{\mathbf{d}}c_{j_1 \cdots j_N} | j_1,\dots,j_N \rangle \nonumber, \end{equation} \tag{ A.1 }$$

where c_j1···jN is a set of d^N complex amplitudes and |j_m〉 is a basis spanning the local d-dimensional Hilbert space of site m. An alternative representation known as a Schmidt decomposition follows by breaking the system into two blocks to the left {1,...,n} and right {n + 1,...,N} of site n:

$$\begin{equation} | \psi \rangle=\sum_{\alpha=1}^{\chi_n}\lambda^{[n]}_{\alpha} | L^{n}_{\alpha} \rangle | R^{n}_{\alpha} \rangle. \nonumber \end{equation} \tag{ A.2 }$$

Here λ^[n]_α are non-negative numbers called Schmidt coefficients and |Lⁿ_α〉 and |Rⁿ_α〉 span the left and right subsystems, respectively. The number of Schmidt coefficients χ_n is called the Schmidt rank. Importantly, the entropy of the probabilities (λ^[n]_α)² is a measure of the entanglement of the state. It is known that for one-dimensional systems with a Hamiltonian composed of nearest-neighbour terms, the ground and low-lying excited states have a very small amount of entanglement resulting in their Schmidt coefficients λ^[n]_α decaying rapidly with α when arranged in descending order, allowing us to heavily truncate the description at a much smaller fixed χ while still providing a near exact description.

Following [43] the state in equation (A.1) can be brought into the canonical form

$$\begin{equation} c_{j_1 j_2 \cdots j_N}=\sum_{\{\alpha\}=1}^{\{\chi\}}\Gamma^{[1]j_1}_{\alpha_1}\lambda^{[1]}_{\alpha_1} \Gamma^{[2]j_2}_{\alpha_1\alpha_2}\lambda^{[2]}_{\alpha_2}\cdots\Gamma^{[N]j_N}_{\alpha_{N-1}}, \end{equation} \tag{ A.3 }$$

where {α} = {α₁,...,α_N−1} and {χ} = {χ₁,...,χ_N−1}. Here Γ^[n]j_n is a set of χ_n−1 × χ_n matrices indexed by the local state j_n, while λ^[n] is a diagonal matrix containing the Schmidt coefficients for the nth bipartition. This is a matrix product description since once all the local configurations j₁j₂···j_N are specified the above sum over all the α's implements a sequence of matrix multiplications. By fixing the Schmidt ranks χ_n ⩽ χ for all bipartitions n the number of parameters needed to describe |ψ〉 scales as O(dNχ²). We note that a relatively modest matrix dimension χ is sufficient for accurate trajectory 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.

A.3. The time evolving block decimation method

A.4. Polariton crystallization calculations