Quantum simulation with interacting photons

In cavity arrays, the individual cavities are coupled via tunnelling of photons between them due to the overlap of the spacial profile of their resonance modes, c.f. figure 1. Following reference [259] an array of cavities can be described by a periodic dielectric constant, $\epsilon (\vec{r})=\epsilon (\vec{r}+R\vec{n})$, where $\vec{r}$ is the three dimensional position vector, R the lattice constant (i.e. the distance between the centres of adjacent cavities) and $\vec{n}$ a vector of integers. Expanding the electromagnetic field in terms of Wannier functions ${w}_{\vec{R}}$, localised in the cavities at locations $\vec{R}=R\vec{n}$, the tunnelling rate of photons between neighbouring cavities can be expressed as

$$\begin{eqnarray}&&{J}_{{phot}}=2{\omega }_{C}\int {d}^{3}r[{\epsilon }_{\vec{R}}(\vec{r})\,-\,\epsilon (\vec{r})]\,{\vec{w}}_{\vec{R}}^{\star }{\vec{w}}_{\vec{R}^{\prime} },\end{eqnarray} \tag{ 1 }$$

where $| \vec{R}-\vec{R}^{\prime} | =R$, ${\omega }_{C}$ is the resonance frequency of the considered cavity mode and the dielectric function ${\epsilon }_{\vec{R}}(\vec{r})$ describes a single cavity surrounded by bulk material only.

Zoom In Zoom Out Reset image size

Introducing creation and annihilation operators ${a}_{\vec{R}}^{\dagger }$ and ${a}_{\vec{R}}$ for the Wannier modes, the Hamiltonian of the field can thus be written as

$$\begin{eqnarray}&&{ \mathcal H }={\omega }_{C}\displaystyle \sum _{\vec{R}}{a}_{\vec{R}}^{\dagger }{a}_{\vec{R}}+{J}_{{phot}}\displaystyle \sum _{\langle \vec{R},\vec{R}^{\prime} \rangle }{a}_{\vec{R}}^{\dagger }{a}_{\vec{R}^{\prime} },\end{eqnarray} \tag{ 2 }$$

where ${\sum }_{\langle \vec{R},\vec{R}^{\prime} \rangle }$ is a sum over nearest neighbours. Since the tunnelling rate is typically much less than the photon frequency, ${J}_{{phot}}\ll {\omega }_{C}$, a rotating wave approximation has been applied. Equation (2) assumes that all the cavities have the same resonance frequency and that the tunnelling rate is the same for all cavity-cavity interactions. In practise there will always be some disorder in the array and the resonance frequencies, ${\omega }_{C}(\vec{R})$, and tunnelling rates, ${J}_{{phot}}(\vec{R},\vec{R}^{\prime} )$ will differ from cavity to cavity. This disorder in the array is a significant challenge for experimental realisations, see section 2.1.7 for further discussion. On the other hand it can even give rise to interesting effects such as the emergence of glassy phases, see section 4.

Although the Hamiltonian (2) has applications for guiding light-modes in 'coupled resonator optical waveguides', see [259], it is of limited interest in the context of quantum simulation. Indeed, the model is fully harmonic and can thus be solved efficiently by for example deriving equations of motion for the first and second order moments of the creation and annihilation operators in the Heisenberg picture. Such non-interacting models thus do not feature a significant quantum complexity and no quantum simulators are required to obtain understanding of them. Nonetheless several interaction free models have recently received significant interest for experiments as they can model artificial gauge fields and give rise to peculiar band structures and chiral edge modes. We discuss such models in more detail in section 5.

The complexity of a model in turn grows dramatically if interactions between excitations play a role and usually no exact solutions are known in such cases. Here we in particular consider on-site interactions as their interplay with the tunnelling leads to interesting many-body physics. A paradigm for such a situation is the Bose–Hubbard Hamiltonian,

$$\begin{eqnarray}{H}_{{BH}} & = & \mu \displaystyle \sum _{\vec{R}}{p}_{\vec{R}}^{\dagger }{p}_{\vec{R}}-J\displaystyle \sum _{\langle \vec{R},\vec{R}^{\prime} \rangle }({p}_{\vec{R}}^{\dagger }{p}_{\vec{R}^{\prime} }+{\rm{H.c.}})\\ & & +\displaystyle \frac{U}{2}\displaystyle \sum _{\vec{R}}{p}_{\vec{R}}^{\dagger }{p}_{\vec{R}}({p}_{\vec{R}}^{\dagger }{p}_{\vec{R}}-1),\end{eqnarray} \tag{ 3 }$$

where ${p}_{\vec{R}}^{\dagger }$ creates a boson at site $\vec{R}$, J is the hopping rate, U the on-site interaction strength, and μ the chemical potential. Motivated by its application to Josephson junction arrays [78], this model received significant interest in the 1990s and its ground state phase diagram has been discussed by Fisher et al [85]. It is characterised by two different phases at zero temperature, an incompressible Mott insulating phase in the interaction dominated regime, $U\gt J$, with commensurate filling and a superfluid phase elsewhere. The phase transition between both phases has been observed in a seminal experiment with ultra-cold atoms in an optical lattice [28, 29, 94].

Approaches to generate an effective Bose–Hubbard model for quantum simulation purposes with photonic excitations consider polaritons that are formed by photons which either interact with atoms [108] or with quantum well excitons [248]. We first discuss setups involving atom-photon interactions.

A Bose–Hubbard model can be generated in an array of cavities that are filled with atoms of a specific four-level structure [108] and driven by a laser in the same manner as in Electromagnetically Induced Transparency (EIT) [86], see figure 2. The transitions between levels $| 2\rangle$ and $| 3\rangle$ couple to the laser field whereas the transitions $| 2\rangle \leftrightarrow | 4\rangle$ and $| 1\rangle \leftrightarrow | 3\rangle$ couple to the cavity resonance mode. Levels $| 1\rangle$ and $| 2\rangle$ are assumed to be metastable with negligible spontaneous emission rates.

Zoom In Zoom Out Reset image size

As has been shown by Imamŏglu and co-workers, this atom cavity system can exhibit a very large optical nonlinearity [125, 257] leading to the phenomenon of photon blockade, where an input drive that is resonant to the first excitation in the system can only generate one excitation which needs to decay before a subsequent excitation can be generated [125, 257]. Photon blockade has so far been experimentally shown with coherent [26, 32, 75, 77, 152, 207] as well as incoherent driving [114].

For one cavity filled with N atoms of the level structure sketched in figure 2 the Hamiltonian describing the atom-photon interactions reads

$$\begin{eqnarray}{H}^{I} & = & \displaystyle \sum _{j=1}^{N}(\varepsilon {\sigma }_{j}^{22}+\delta {\sigma }_{j}^{33}+({\rm{\Delta }}+\varepsilon ){\sigma }_{j}^{44})\\ & & +\displaystyle \sum _{j=1}^{N}({\rm{\Omega }}\,{\sigma }_{j}^{23}\,+\,{g}_{13}\,{\sigma }_{j}^{13}\,{a}^{\dagger }\,+\,{g}_{24}\,{\sigma }_{j}^{24}\,{a}^{\dagger }\,+\,{\rm{H.c.}}),\end{eqnarray} \tag{ 4 }$$

in a rotating frame with respect to ${H}_{0}={\omega }_{C}\left({a}^{\dagger }a+\tfrac{1}{2}\right)+{\sum }_{j=1}^{N}({\omega }_{C}{\sigma }_{j}^{22}+{\omega }_{C}{\sigma }_{j}^{33}+2{\omega }_{C}{\sigma }_{j}^{44})$. Here ${\sigma }_{j}^{{kl}}=| {k}_{j}\rangle \langle {l}_{j}|$ is the transition operator between levels $| l\rangle$ and $| k\rangle$ of atom j, ${\omega }_{C}$ is the frequency of the cavity mode, Ω is the Rabi frequency of the laser drive and g₁₃ (g₂₄) are the coupling strengths of the cavity mode to the atomic transitions $| 2\rangle \leftrightarrow | 4\rangle$ ($| 1\rangle \leftrightarrow | 3\rangle$).

A cavity array as described in section 2.1.1, where each cavity is doped with N four-level atoms as depicted in figure 2, can form a quantum simulator for a Bose–Hubbard model, c.f. equation (3), if all atoms interact in the same way with the cavity mode and the number of excitations is significantly lower than the number of atoms in each cavity. In this regime, the dynamics generated by the Hamiltonian (4) can be described in terms of polaritons, hybridised light-matter quasi-particles that obey bosonic statistics. One species of these polaritons does only occupy the atomic levels which do not have a direct dipole transition to the atomic ground state. These polaritons have been considered by Fleischhauer et al [86] and are called dark state polaritons as they do not lead to emission of radiation and are therefore long lived. The latter property is obviously rather beneficial for quantum simulation applications. The creation operators of the dark state polaritons read

$$\begin{eqnarray}&&{p}^{\dagger }=\displaystyle \frac{1}{B}({{gS}}_{12}^{\dagger }-{\rm{\Omega }}{a}^{\dagger })\end{eqnarray} \tag{ 5 }$$

where $g=\sqrt{N}{g}_{13}$ is a collective coupling rate, ${S}_{12}^{\dagger }=\tfrac{1}{\sqrt{N}}{\sum }_{j=1}^{N}{\sigma }_{j}^{21}$ creates a spin wave in the metastable atomic levels and $B=\sqrt{{g}^{2}+{{\rm{\Omega }}}^{2}}$. The dynamics of the dark state polaritons decouples from the remaining species for a suitable parameter regime, where the frequencies of the species are sufficiently separated.

The coupling of the dark state polaritons to level $| 4\rangle$ of the atoms induces an effective interaction between all dark state polaritons in one cavity. For $| {g}_{24}g{\rm{\Omega }}/{B}^{2}| \ll | {\rm{\Delta }}|$ the strength of this interaction can be calculated in a perturbative manner [257],

$$\begin{eqnarray}&&{U}_{{BH}}=\displaystyle \frac{2{g}_{24}^{2}}{{\rm{\Delta }}}\,\displaystyle \frac{{{Ng}}_{13}^{2}\,{{\rm{\Omega }}}^{2}}{{({{Ng}}_{13}^{2}+{{\rm{\Omega }}}^{2})}^{2}}.\end{eqnarray} \tag{ 6 }$$

Similarly, the two photon detuning ε leads to an energy shift of ${\mu }_{{BH}}=\epsilon {g}^{2}/{B}^{2}$ for the polaritons that plays a similar role as a chemical potential, see also section 4.1.1.

Provided the hopping rate of photons between cavities is small compared to the frequency separation between the polariton species, the hopping of photons translates into a hopping of dark state polaritons at a rate [108]

$$\begin{eqnarray}&&{J}_{{BH}}=\displaystyle \frac{2{{\rm{\Omega }}}^{2}}{{{Ng}}_{13}^{2}+{{\rm{\Omega }}}^{2}}{J}_{{phot}}\,\end{eqnarray} \tag{ 7 }$$

where J_phot is defined in equation (1).

As a consequence the Hamiltonian for the dark state polaritons takes on the form of a Bose–Hubbard model as introduced in equation (3), with $J={J}_{{BH}}$ (c.f. equation (7)), $U={U}_{{BH}}$ (c.f. equation (6)) and $\mu ={\mu }_{{BH}}$. Note that both, J_BH and U_BH can be tuned by varying the intensity of the laser driving ${{\rm{\Omega }}}^{2}$. See also [106] for a generalisation to two polariton species and [33, 188, 210] for alternative nonlinearities.

Concepts for quantum simulators of Bose–Hubbard models have also been put forward for a rather different experimental platform based on semiconductor structures. We discuss these in the next section.

An interaction of the form as in the Bose–Hubbard model, see equation (3), was also shown to exists and analysed for polaritons formed by photons and excitons of a quantum well [248]. The Hamiltonian describing the interactions between the photons of one cavity mode and the excitons of the matching wave-vector here reads,

$$\begin{eqnarray}&&H={\omega }_{X}{b}^{\dagger }b+\displaystyle \frac{{U}_{X}}{2}{b}^{\dagger }{b}^{\dagger }{bb}+{\omega }_{C}{a}^{\dagger }a+g({{ab}}^{\dagger }+{a}^{\dagger }b),\end{eqnarray} \tag{ 8 }$$

where ${a}^{\dagger }(a)$ creates (annihilates) a cavity photon and ${b}^{\dagger }(b)$ creates (annihilates) an exciton. The excitons have transition frequencies ${\omega }_{X}$, interact with each other at a strength U_X and couple to photons at a strength g. An anharmonic exciton-photon coupling that depends on the exciton oscillator strength saturation density can typically be neglected compared to the interaction between excitons [248].

In a so called 'strong coupling regime' where g is larger than the linewidth of the cavity and the excitons, polaritons, i.e. hybridised light-matter excitations, become the elementary excitations of the system. These have annihilation operators

$$\begin{eqnarray}{p}_{+} & = & \cos (\theta )a+\sin (\theta )b\\ {p}_{-} & = & \sin (\theta )a-\cos (\theta )b,\end{eqnarray} \tag{ 9 }$$

and frequencies

$$\begin{eqnarray}&&{\omega }_{\pm }=\displaystyle \frac{{\omega }_{X}+{\omega }_{C}}{2}\pm \displaystyle \frac{\sqrt{{({\omega }_{C}-{\omega }_{X})}^{2}+{g}^{2}}}{2}\end{eqnarray} \tag{ 10 }$$

with $\sin (\theta )=g/\sqrt{{g}^{2}+{\widetilde{{\rm{\Delta }}}}^{2}}$, $\cos (\theta )=\widetilde{{\rm{\Delta }}}/\sqrt{{g}^{2}+{\widetilde{{\rm{\Delta }}}}^{2}}$ and $\widetilde{{\rm{\Delta }}}={\omega }_{C}-{\omega }_{X}+\sqrt{{({\omega }_{C}-{\omega }_{X})}^{2}\,+\,{g}^{2}}$.

Since $g\gg {U}_{X}$ the dynamics of both polariton species decouples from each other and the interaction between excitons leads to an interaction $({U}_{p-}/2){p}_{-}^{\dagger }{p}_{-}^{\dagger }{p}_{-}{p}_{-}$ with ${U}_{p-}={U}_{X}{\cos }^{4}(\theta )$ for the lower polaritons ${p}_{-}$. These interactions can have a strength comparable or even larger than their linewidth due to the strong interactions of the excitonic components of the polaritons. Ways to enhance these polariton-polariton interactions via Feshbach type resonances between polariton pairs and bi-excitons have been explored by Carusotto et al [41].

By confining the light modes in a periodic structure, a lattice model as given in equation (3) can be generated [37]. The periodic confinement for the photons can thereby be either a photonic crystal or a structure of connected micro-pillars. Such polariton-polariton interactions and their interplay with polariton tunnelling between lattice sites has recently been investigated in a self trapping experiment [1], see section 2.1.8 for a more detailed discussion.

Besides Bose–Hubbard physics, most of the research effort on lattice models with optical photons has considered scenarios where the photons couple to one two-level system in each cavity. We review these approaches in the next section.

An effective repulsive interaction between photons or polaritons can also be generated by doping a cavity with one atom of which only one internal transition couples to the cavity resonance mode. This situation of a two-level emitter coupled to a single cavity mode is described by the celebrated Jaynes–Cummings model [127],

$$\begin{eqnarray}&&{H}^{{JC}}={\omega }_{C}{a}^{\dagger }a+{\omega }_{0}| e\rangle \langle e| +g({a}^{\dagger }| g\rangle \langle e| +a| e\rangle \langle g| ),\end{eqnarray} \tag{ 11 }$$

where ${\omega }_{C}$ and ${\omega }_{0}$ are the frequencies of the cavity mode and the atomic transition, g is Jaynes–Cummings light-matter coupling, ${a}^{\dagger }$ is the creation operator of a photon in the resonant cavity mode, and $| g\rangle (| e\rangle )$ are the ground (excited) states of the two level system. Since the energy of a photon ${\omega }_{C}$ and the atomic transition energy ${\omega }_{0}$ are much greater than the coupling g, the number of excitations is conserved by the Hamiltonian (11). Hence it can be diagonalised for each manifold with a fixed number of excitations n separately. The energy eigenvalues E_n for n excitations read ${E}_{0}=0$ and ${E}_{n}^{\pm }=n{\omega }_{C}+\tfrac{{\rm{\Delta }}}{2}\pm \sqrt{{{ng}}^{2}+\tfrac{{{\rm{\Delta }}}^{2}}{4}}$ for $n\geqslant 1$, where ${\rm{\Delta }}={\omega }_{0}-{\omega }_{C}$. As a consequence, the energy of the lowest state with two excitations is not twice the energy of a single excitation state and their difference

$$\begin{eqnarray}&&{U}_{{JC}}={E}_{2}^{-}-2{E}_{1}^{-}=2\sqrt{{g}^{2}+\displaystyle \frac{{{\rm{\Delta }}}^{2}}{4}}-\sqrt{2{g}^{2}+\displaystyle \frac{{{\rm{\Delta }}}^{2}}{4}}-\displaystyle \frac{{\rm{\Delta }}}{2}\end{eqnarray} \tag{ 12 }$$

plays the role of an effective on-site repulsion that can be tuned via the detuning Δ. Photon blockade due to the effective repulsion (12) has been observed for an individual atom [26], quantum dots in photonic crystals [75, 77, 207] and a circuit quantum electrodynamics (circuit QED) setup [32, 152], see also section 3 for further details.

The concept of photon blockade in a Jaynes–Cummings model can also be explored for higher excitation numbers, where the difference between the energy of $n+1$ excitations, ${E}_{n+1}^{-}$, and the energies of n excitations, ${E}_{n}^{-}$, and a single excitation, ${E}_{1}^{-}$, should be considered. For ${\rm{\Delta }}=0$ one finds

$$\begin{eqnarray}&&{U}_{{JC};n}={E}_{n+1}^{-}-({E}_{n}^{-}+{E}_{1}^{-})=g\,(\sqrt{n}+1-\sqrt{n+1}),\end{eqnarray} \tag{ 13 }$$

so that the interaction energy per excitation decreases with growing excitation number n, ${U}_{{JC};n}/n\to g/n$ for $n\gg 1$. The degrading of the anharmonicity of the spectrum of the Jaynes–Cummings model at high excitation numbers was observed in a self-trapping experiment [204], see section 2.1.8, and the resulting breakdown of the photon blockade effect was investigated by Carmichael [38] as an example for a dissipative quantum phase transition. He found the transition to be first order, as indicated by a bi-modality of the Q-function, except for a critical value of drive strength for zero drive detuning, where it is a continuous second order transition.

In an array of cavities where each cavity is doped with a single two-level emitter that is coupled to the light mode as described by equation (11), the effective repulsion of equation (12) will suppress the mobility of excitations in a similar manner as the on-site interactions in the Bose–Hubbard model since the lowest energy for two excitations in one cavity is ${E}_{2}^{-}$ and moving one additional excitation to a cavity requires an extra energy of U_JC. This setup has first been investigated by Angelakis et al [4] and Greentree et al [93]. The full effective many-body model in a cavity array that is based on the effective interaction (12) has been coined Jaynes–Cummings-Hubbard model and reads,

$$\begin{eqnarray}&&{H}_{{JCH}}=\displaystyle \sum _{\vec{R}}{H}_{\vec{R}}^{{JC}}-{J}_{{JC}}\displaystyle \sum _{\langle \vec{R},\vec{R}^{\prime} \rangle }({a}_{\vec{R}}^{\dagger }{a}_{\vec{R}^{\prime} }+{\rm{H.c.}}),\end{eqnarray} \tag{ 14 }$$

where $\vec{R}$ labels the site of a cavity. Each cavity contains one atom interacting via the Jaynes–Cummings interaction with the cavity mode and photons tunnel between neighbouring cavities at a rate ${J}_{{JC}}={J}_{{phot}}$. The dispersive regime of the Jaynes–Cummings-Hubbard and Rabi-Hubbard model was investigated in [262].

Multiple two-level atoms per cavity. Effective many-body physics based on the Hamiltonian (14) can not only be observed for a single two level system in each cavity, as assumed in equation (11), but also for setups with several two level systems per cavity. Such a model can describe photonic crystal micro-cavities doped with substitutional donor or acceptor impurities. This approach to implementing effective many-body models, which can have suitable parameters, has been proposed in [186]. The phase transitions of a model with several two level systems in each cavity have also been studied in [157, 212], see section 4 for further discussions.

Although it is not the main topic of this review, we note here that coupled cavity arrays have also been considered for the simulation of spin lattice Hamiltonians [4, 8, 92, 109, 112]. Here the internal levels of the atoms in the cavities represent the spin degrees of freedom and interactions between spins can be mediated via off-resonant couplings to collective photon modes [8, 92, 109, 112, 262].

After having discussed the main strands of the theory work on simulating quantum lattice Hamiltonians with optical photons, we now take a closer look at the requirements and challenges for implementing sthese approaches in experiments.

Electromagnetic excitations, including polaritons and photons, inevitably couple to the electromagnetic vacuum that can not be excluded from any experiment. These excitations will therefore always by limited to a finite lifetime or trapping-time. In order to be able to explore their dynamics, they need to be kept in the experimental sample for a time that exceeds the timescales associated to the kinetic and interaction energies. Denoting the rate of photon losses from the cavities by κ and the rate of spontaneous emission for the atoms by γ, one needs

$$\begin{eqnarray}&&{U}_{{BH}},{J}_{{BH}}\gt \kappa ,\gamma \quad {\rm{or}}\quad {U}_{{JC}},{J}_{{JC}}\gt \kappa ,\gamma \end{eqnarray} \tag{ 15 }$$

for effective Bose–Hubbard or Jaynes–Cummings-Hubbard models. In terms of the light matter couplings, meeting the conditions (15) requires both transitions to operate at high cooperativity for single photons, ${g}_{13}^{2}\gg \kappa \gamma$ and ${g}_{24}^{2}\gg \kappa \gamma$, for the Bose–Hubbard model, or a strong coupling regime, $g\gg \kappa ,\gamma$ for the Jaynes–Cummings-Hubbard model. For the approach to realise a Bose–Hubbard model as explained in section 2.1.2 the condition ${\rm{\Delta }}\gt \gamma$ is also needed. The fact that for this model only the product of the two decay rates κ and γ needs to be bounded can be understood by realising that the relative contributions of the photonic and atomic components in the dark state polaritons, c.f. equation (5) can be varied to avoid the faster decay channel and optimise their lifetime.

An alternative approach to the photon blockade effect was discovered in a two-site Bose–Hubbard system where a resonantly driven nonlinear cavity is tunnel-coupled to an auxiliary cavity. Remarkably this leads to anti-bunched output photons even if $U\ll \kappa$ [167] but requires low rates of pure dephasing ${\gamma }^{* }$ [80]. The effect is due to a destructive interference of two excitation paths as can be seen from an expansion in excitation numbers at low drive intensities [17]. The direct path to generate a second excitation in cavity one via the coherent input, $| 1,0\rangle \to | 2,0\rangle$, destructively interferes with the excitation path where the first excitation tunnels to the second cavity, $| 1,0\rangle \to | 0,1\rangle$, a second excitation is generated in cavity one via the drive, $| 0,1\rangle \to | 1,1\rangle$, and the excitation in cavity two tunnels back to cavity one. Its origin in an interference of excitation paths explains why this photon blockade effect requires that the nonlinearity exceeds the rate of pure dephasing $U\gg {\gamma }^{* }$. Despite these requirements for an experimentally suitable technology, there has already been significant progress as we discuss next.

Demonstrating coherent coupling between high finesse optical resonators with low enough disorder and light-matter interactions in the strong coupling regime at the same time remains an experimental challenge at the time of writing. Yet, coupled arrays of photonic band-gap cavities in photonic crystal structures that are suitable for quantum simulation applications have been built with 10-20 cavities [126, 175]. A further possibility to engineer a tunnel coupling between adjacent cavities is to connect these via waveguides that come close to each other at a coupling point [162]. In this way a controllable coupling can be combine with a small mode volume but open cavity that allows to strongly couple laser trapped atoms to its resonance modes.

Realization of the Jaynes–Cummings-Hubbard model with trapped ions. A minimal version of the Jaynes–Cummings-Hubbard model with two lattice sites has been realised with trapped ions [244]. Here the motion of the ions was employed to implement the harmonic degree of freedom thus replacing the cavity modes and the internal levels of the ions represent the two-level systems. Hence the creation and annihilation operators ${a}_{j}^{\dagger }$ and a_j here describe phonons rather than photons. In the experiment, an adiabatic sweep from a regime dominated by the on-site repulsion U_JC to a regime dominated by the tunnelling J_JC was performed via varying the detuning Δ by tuning the frequency of a laser that drives a red side-band transition for the ions. The mobility of the polaritons in the tunnelling dominated regime and the suppression of their mobility in the interaction dominated regime was then shown by measuring their on-site number fluctuations.

Self-trapping experiments. The interplay of interactions and tunnelling between adjacent lattice sites has also been investigated in two experiments with two coupled resonators [1, 204]. These experiments explored self-trapping effects in regimes of high excitation densities, an interaction phenomenon that already becomes accessible for moderate interactions between individual excitations. The effect appears for two coupled resonators with a strong imbalance of excitation numbers, where the interaction energies per excitation differ between both resonators by an amount that exceeds the rate of inter-resonator tunnelling.

For a nonlinearity of the form ${Un}(n-1)$ [n is the number of particles], which has been realised in an experiment with exciton polaritons in two coupled Bragg stack micro-pillars [1], the interaction energy per particle is $U(n-1)$ and self trapping occurs for high particle densities but ceases as the particle number decays. In their experiment, Abbrachi et al [1] thus observed a transition from a self-trapped regime to a regime of excitation oscillations between the resontors as the particle number decreased over time due to dissipation, see figures 3(a)–(c).

Zoom In Zoom Out Reset image size

For interactions as present in the Jaynes–Cummings model, the interaction energy per excitation degrades as the number of excitations grows, see equation (13). One thus observes oscillations of excitations between both resonators for a strong initial excitation number imbalance [229]. This has been seen by Raftery et al in an experiment with two coupled superconducting coplanar waveguide resonators that each interact with a transmon qubit [204], c.f. Section 3. As the excitation number decreased a transition to the self-trapped regime was observed, see figures 3(d)–(e).

After reviewing work on the simulation of quantum lattice models, we now turn to discuss efforts towards quantum simulators for continuum models. Importantly, for photons at optical frequencies these approaches face less experimental challenges.

A realisation of the model (16) with dark-state or slow-light polaritons has been proposed by Chang et al [45]. The approach generalises the concept as explained in section 2.1.2 to continuous models and combines it with a so called 'stationary light' regime [15] generated by two counter-propagating control fields. A possible implementation could be a hollow-core photonic crystal fibre filled with a gas of Doppler cooled atoms [14], see figure 4, or atoms trapped in the evanescent field of an optical fibre that is tapered down to a diameter comparable or below the wavelength of the employed light [249].

Zoom In Zoom Out Reset image size

An approach to exploiting unitary as well as dissipative interactions has been introduced by Kiffner et al [138] by decomposing the field ψ into momentum modes. The operator that excites a dark-state polariton in momentum k is here defined as

$$\begin{eqnarray}&&{\psi }_{k}=\cos \theta \displaystyle \frac{{a}_{{k}_{c}+k}+{a}_{-{k}_{c}+k}}{\sqrt{2}}-\displaystyle \frac{\sin \theta }{\sqrt{N}}\displaystyle \sum _{\mu =1}^{N}{\sigma }_{\mu }^{12}{e}^{-{{ikz}}_{\mu }},\end{eqnarray} \tag{ 18 }$$

where k_c is the momentum of the control fields, $\sin \theta =\sqrt{N}{g}_{13}/\sqrt{{{Ng}}_{13}^{2}+2{{\rm{\Omega }}}^{2}}$, $\cos \theta =\sqrt{2}{{\rm{\Omega }}}_{c}$ $/\sqrt{{{Ng}}_{13}^{2}+2{{\rm{\Omega }}}^{2}}$ and ${N}^{-1/2}{\sum }_{\mu =1}^{N}{\sigma }_{\mu }^{12}{e}^{-{{ikz}}_{\mu }}$ describes a spin coherence. The dynamics of the dark state polaritons can be described in terms of their density matrix ${\rho }_{D}$. Neglecting single particle dissipation which can be sufficiently suppressed for suitable parameters, ${\rho }_{D}$ obeys the equation of motion [138, 139],

$$\begin{eqnarray}&&{\hslash }{\dot{\rho }}_{D}=-{{iH}}_{{\rm{eff}}}{\rho }_{D}+i{\rho }_{D}{H}_{{\rm{eff}}}^{\dagger }+{ \mathcal I }{\rho }_{D},\end{eqnarray} \tag{ 19 }$$

where ${H}_{{\rm{eff}}}={H}_{{LL}}$ with effective mass ${m}_{{\rm{eff}}}=-{\hslash }({{Ng}}_{13}^{2}+2{{\rm{\Omega }}}^{2})/(2\delta {c}^{2}{\cos }^{2}\theta )$ and a complex interaction constant $\tilde{g}=2{\hslash }{{Lg}}_{24}^{2}{\cos }^{2}\theta /({\rm{\Delta }}-{\cos }^{2}\theta {\rm{\Delta }}\omega +i{\gamma }_{42}/2)$. Here c is the speed of light. The term

$$\begin{eqnarray}&&{ \mathcal I }{\rho }_{D}=-{\rm{Im}}(\tilde{g}){\int }_{0}^{L}{dz}{\psi }^{2}{\rho }_{D}{\psi }^{\dagger 2}.\end{eqnarray} \tag{ 20 }$$

describes correlated decay processes, where always two polaritons are simultaneously lost. Equation (19) thus has the form of a Lindblad master equation where the jump operator describes correlated decays of polariton pairs.

For the realisation as described by equation (19), the absolute value of the Lieb-Liniger Parameter is

$$\begin{eqnarray}&&| G| =\displaystyle \frac{{g}_{13}^{2}{g}_{24}^{2}{L}^{2}N}{{c}^{2}| \delta | \sqrt{{{\rm{\Delta }}}^{2}+{\gamma }_{42}^{2}/4}\,{N}_{p}}\end{eqnarray} \tag{ 21 }$$

and is maximal for purely dissipative (${\rm{\Delta }}=0$) interactions between the polaritons [138], see also [102].

The strongly correlated regime with $| G|$ larger than unity thus becomes accessible for an optical depth per atom that exceeds 160. For the density-density correlations, one finds ${g}^{(2)}(z,z)=(1-1/{N}_{{\rm{ph}}}^{2})4{\pi }^{2}/(3| G{| }^{2})$ for $z=z^{\prime}$ which vanishes in the limit $| G| \to \infty$. Moreover, in this strongly correlated regime, the ground state of this generalised Lieb-Liniger model is the same as in the original model with repulsive interaction [70], indicating a crystallisation of the polaritons.

A generalisation of the above approaches to a situation with two atomic species filling the hollow core fibre, c.f. figure 4, was considered by Angelakis et al [5]. These conditions give rise to two polariton species ${\psi }_{1}$ and ${\psi }_{2}$ that are each described by a Lieb-Liniger Hamiltonian as in equation (16) and in addition are subject to a density-density interaction of the form

$$\begin{eqnarray}&&{\int }_{0}^{L}{{dzV}}_{12}{\psi }_{1}^{\dagger }{\psi }_{1}{\psi }_{2}^{\dagger }{\psi }_{2},\end{eqnarray} \tag{ 22 }$$

with strength V₁₂ [5]. The scenario can thus emulate Luttinger liquid behaviour and allows for exploring an analogue of spin-charge separation due to the mapping between hard-core bosons and free fermions in one dimension [91]. Subsequently applications of this approach to simulate Cooper pairing [122] and relativistic field theories [6] have been investigated, see also [190].

For the above approaches to Lieb-Liniger physics with dark state polaritons, the effect of dissipative interactions and the build-up of Friedel oscillations was studied numerically in [140]. Moreover photonic transport in a fibre as sketched in figure 4 was studied in [100, 102] via a Bethe ansatz solution for the driven Lieb-Liniger model where the setup was shown to act as a single photon switch for repulsive interactions and able to support two-photon bound states for attractive interactions. As we will discuss next, there has already been significant progress in experimental observations of the physics of such continuous one-dimensional models.

There has been substantial progress towards realising quantum many-body systems of strongly interacting photons with optical fibres, although the simulation of a Lieb-Liniger model appears to still pose challenges. Laser-cooled ⁸⁷Rb atoms have been loaded int a hollow core photonic crystal fibre by Bajcsy et al [14] to demonstrate all optical switching. Via a switching beam coupled to the transition $2\leftrightarrow 4$, see figure 2, the initial transmission of a beam coupled to the transition $1\leftrightarrow 3$ was reduced by 50%. Another approach trapped laser-cooled neutral atoms with a multicolour evanescent field surrounding an optical nano-fibre and showed appreciable coupling of the atoms to fibre guided light modes [249]. The concept was then developed further to demonstrate switching of optical signals between two optical fibres [194] and nanophotonic optical isolators [224]. As achieving sufficiently strong optical nonlinearities for simulating strongly correlated quantum many-body systems with optical photons remains a challenge in these devices, atomic Rydberg media are now being considered more intensely.

In contrast to optical photons interacting with solid state emitters or atoms, microwave photons in superconducting circuits can be made to interact non-locally [128]. That is a photon in one resonator can scatter off a photon in a second, even spatially distant, resonator. Such processes are described by cross-Kerr interactions of the form

$$\begin{eqnarray}&&V\,{a}_{1}^{\dagger }{a}_{1}\,{a}_{2}^{\dagger }{a}_{2},\end{eqnarray} \tag{ 26 }$$

where ${a}_{1}({a}_{2})$ annihilates a photon in resonator 1(2). The circuit that generates this type of interactions together with correlated tunnelling processes is sketched in figure 8. The dc-SQUID connecting the two resonators gives rise to an energy contribution of $-{E}_{J}\cos ({\varphi }_{1}-{\varphi }_{2})$, where ${\varphi }_{1(2)}\propto {a}_{1(2)}+{a}_{1(2)}^{\dagger }$, and an expansion up to 4th order in ${\varphi }_{1}-{\varphi }_{2}$ leads to

$$\begin{eqnarray}{H}_{{cK}} & = & {\hslash }\widetilde{\alpha }\omega ({a}_{1}^{\dagger }{a}_{2}+{a}_{1}{a}_{2}^{\dagger })-2\alpha {E}_{C}{a}_{1}^{\dagger }{a}_{1}{a}_{2}^{\dagger }{a}_{2}\\ & & +\alpha {E}_{C}({a}_{i}{a}_{j}^{\dagger }{a}_{j}^{\dagger }{a}_{j}+{a}_{i}^{\dagger }{a}_{i}^{\dagger }{a}_{i}{a}_{j}\\ & & \left.-\displaystyle \frac{1}{2}{a}_{i}^{\dagger }{a}_{i}^{\dagger }{a}_{j}{a}_{j}+{\rm{H.c.}}\right),\end{eqnarray} \tag{ 27 }$$

where ${E}_{C}={e}^{2}/[2(C+2{C}_{J})]$ and $\alpha ={C}_{J}/(C+2{C}_{J})$ with C the capacitance in each resonator and C_J the capacitance that shunts the dc-SQUID. The coefficient $\widetilde{\alpha }$ can be made vanishingly small by tuning the SQUID to the point, where photon tunnelling via the SQUID and via the shunt capacitance interfere destructively and cancel each other [187]. When all couplings between neighbouring resonators in a large lattice are built as in figure 8, one arrives at a Bose–Hubbard Hamiltonian augmented by the cross-Kerr interactions (26) and correlated tunnelling processes [128]. This scenario can give rise to a density wave type ordering, see also section 4.

Zoom In Zoom Out Reset image size

Lattice elements coupled by Josephson junctions have furthermore been considered for the simulation of Anderson and Kondo lattices [88] and tunable coupling elements [23].

Bose–Hubbard physics in superconducting architectures has been investigated in the early 1990 already, well before the realisations of the model in optical lattices. The investigated structures consisted of Josephson junction arrays where Cooper pairs could tunnel between superconducting islands through the junctions and interact via Coulomb forces on each island, see [78] for a review. The practical difference of the new generation of networks discussed here is that the individual network nodes are separated by larger distances on the chip and can therefore be individually addressed via control lines. The high precision control over the individual nodes allows to suppress disorder much better than in Josephson junction arrays, where it was a significant limitation to experiments. Yet circuit-QED lattices also allow to emulate further many-body models, such as the Jaynes–Cummings-Hubbard model.

A dimer of two lumped element resonators has been empoyed by Eichler et al to show quantum-limited amplification and entanglement [72]. In their setup, both resonators are nonlinear because their inductors are formed by a series of dc-SQUIDs and are mutually coupled via a capacitance. Here a moderate nonlinearity is desirable to allow for an appreciable bandwidth of the amplification mechanism.

A chain of three capacitively coupled transmon qubits has been realised inside a microwave cavity by Hacohen-Gourgy et al [96]. Here all qubits couple dispersively to a common resonance mode of the cavity and engineered cooling and heating processes are generated by driving the cavity resonance with red or blue detuned input tones. These processes act as a quantum bath that can exchange energy and entropy with the Bose–Hubbard chain but conserves its excitation number, see also section 4.1.1. An extension to more qubits in a common resonator for exploring dipolar spin models has been considered in [56].

Photons can typically only bet trapped for limited times in optical or microwave resonators. The achievable trapping times are moreover comparable or even below the experimental equilibration time-scales. In addition, the interaction of photons with matter involves absorbtion and emission processes. As a consequence of these properties, the number of photons is usually not conserved during an experiment and there is no equivalent to a chemical potential for photons.

In some experimental situations thermalisation of photons via scattering with phonons or repeated absorption-emission cycles was achieved on sufficiently fast time-scales so that a quasi-equilibrium builds up. Prominent examples for such situations are the condensation of exciton polaritons [57, 133, 135, 143, 237] in a semiconductor microcavity or the number-conserving thermalisation and Bose–Einstein condensation of a two-dimensional photon gas in a dye-filled optical microcavity [144]. Despite the observation of the characteristic features of Bose–Einstein condensation these non-equilibrium cases are distinct from the equilibrium situation [261] and creating a proper chemical potential for light remains an open question.

Concepts for generating such a chemical potential for photons have been proposed based on a parametric modulation of the system's coupling to its environment at a frequency that exceeds the highest spectral component of the environment [99]. Moreover, connections to regimes of ultra-strong light matter coupling have been discussed as their ground states can be viewed as containing a non-vanishing photon number which is however not observable without further manipulation [99, 226]. As the generation of a chemical potential for photons remains a challenging task, most discussions of equilibrium phase diagrams simply introduced such a potential by hand. We start our discussion of these works by reviewing those based on mean-field approaches that are expected to be accurate in high dimensional lattices.

The first study by Greentree et al [93] used a mean-field approach that is expected to become increasingly accurate with higher (three or more) lattice dimensions. Introducing $\psi =\langle {a}_{j}\rangle$ as a superfluid order parameter, a mean-field decoupling was performed in the Hamiltonian (14) and a chemical potential μ was added to get

$$\begin{eqnarray}\tilde{H} & = & \displaystyle \sum _{\vec{R}}{H}_{\vec{R}}^{{JC}}-\mu \displaystyle \sum _{\vec{R}}({a}_{\vec{R}}^{\dagger }{a}_{\vec{R}}+{\sigma }_{\vec{R}}^{+}{\sigma }_{\vec{R}})\\ & & -{{zJ}}_{{JC}}\displaystyle \sum _{\vec{R}}({a}_{\vec{R}}^{\dagger }\psi +{\rm{H.c.}})+{{zJ}}_{{JC}}| \psi {| }^{2},\end{eqnarray} \tag{ 28 }$$

where z denotes the number number of neighbouring lattice sites, i.e. z = 6 for three dimensions. Depending on whether the value of ψ that minimises the expectation value of $\tilde{H}$ vanishes or not, the system is in a Mott insulating state or supports a super-fluid component. The resulting phase diagram for the case where the transition frequency of the two-level system ${\omega }_{0}$ equals the cavity resonance frequency ${\omega }_{C}$, ${\rm{\Delta }}={\omega }_{0}-{\omega }_{C}=0$, is reproduced in figure 9. The analogous phase diagram for the Bose–Hubbard model had already been calculated in the late 1980s by Fisher et al [85]. Using a fieldtheoretic approach, Koch and Le Hur showed that the Jaynes–Cummings-Hubbard and Bose–Hubbard models are in the same universality class and that Jaynes–Cummings-Hubbard features multicritical curves, which parallel the presence of multicritical points in the Bose–Hubbard model [148].

Zoom In Zoom Out Reset image size

Comparisons between the mean-field results and exact calculations for small lattices have found that the Mott lobes shrink in finite systems, similar to the Bose–Hubbard model [177]. Moreover, the phase diagram of the Jaynes–Cummings-Hubbard model changes substantially for ultra-strong light matter coupling, where it maps to the transverse field Ising model and features a discrete parity symmetry-breaking transition [226].

Whereas mean-field approaches are expected to become increasingly accurate for higher and higher lattice dimensions, ground states and dynamics of one-dimensional systems can often be efficiently calculated using numerical approaches based on the Density Matrix Renormalisation Group (DMRG) [232]. For the approaches described in sections 2.1.2 and 2.1.5 the ground state phase diagrams (after a chemical potential term had been added) have been obtained with these techniques by Rossini et al [212, 214], see figure 10. Here, the microscopic origin of the nonlinearity in the explicit form of the atomic level structure has been taken into account for the approach to the Bose–Hubbard model, c.f. section 2.1.2, and deviations in the phase-diagram of the model from the original Bose–Hubbard model, c.f. equation (3), have been found due to small atom numbers. Moreover, the existence of a glassy phase due to non-vanishing disorder in the number of atoms per cavity has been predicted.

Zoom In Zoom Out Reset image size

Phase diagrams for one and two-dimensional systems have also been calculated with a variational cluster approach [2, 145, 146] or quantum Monte Carlo calculations [116, 117, 199], and analysed differences between the Bose–Hubbard and Jaynes-Cuminngs-Hubbard models due to the composite nature of the excitations in the latter. An analytic strong-coupling theory based on a linked-cluster expansion for the phase diagram of the Jaynes–Cummings-Hubbard model and its elementary excitations in the Mott phase has been derived by Schmidt et al [227, 228].

Recent work has also predicted quantum phase transitions in finite size and even single site systems of Rabi [124] and Jaynes–Cummings models [123].

Properties of the ground state phase diagrams of Bose–Hubbard and Jaynes–Cummings-Hubbard models can be investigated in non-equilibrium states using the following concept. The system parameters are initially tuned such that photon tunnelling between resonators is strongly suppressed and on-site interactions are very strong. Due to the nonlinearity of their spectra, an input pulse can generate a single excitation in each lattice site and thus prepare the system in a state very similar to a Mott insulating state. Despite not being the ground state, this initial state is an approximate eigenstate of the system for the chosen parameters. By changing the parameters of the system slowly enough to generate an adiabatic sweep, the system will stay in this energy eigenstate, which will however change its character and become very similar to the superfluid ground state when the tunnelling is increased and nonlinearities are decreased. This sequence has been analysed by Hartmann et al [108] in terms of the fluctuations of the polariton number in each lattice site and by Angelakis et al [4] for the excitation number fluctuations in a Jaynes–Cummings lattice. Particle number fluctuations are strongly suppressed in the Mott insulating regime where particles are localised but become finite in the superfluid regime where the particles become delocalised.

A scenario similar to sudden quenches has been explored by Tomadin et al [243], where, independently of the system parameters, the cavity array was assumed to be initialised into a direct product of single excitation Fock states for each cavity by a short but intense input pulse. By exploring the subsequent non-equilibrium dynamics, it was found that the rescaled superfluid order parameter $\psi /\sqrt{n}$ (n is the excitation density) decays to zero in the Mott insulating but remains finite in the superfluid regimes, see also [54] for analogue results for lattices with disorder. Due to the possibilities of local and single-site control offered by many setups for resonator arrays, one can also explore quenches that are not uniform across the lattice. For example if a bipartite lattice is initially prepared in a superfluid regime in one half and a Mott insulating regime in the other half, all excitations migrate to the superfluid half upon switching on a small tunneling between both parts [110]. This example shows that the tendency of non-equilibrium quantum systems to explore all the accessible Hilbert space, which typically results in the formation of local equilibria [73], can lead to states which are strongly inhomogeneous and not translation invariant provided the underlying system breaks such translation invariance.

Due to inevitable experimental imperfections, photons dissipate after a relatively short time from the quantum simulation devices described in this review. These photon losses can be compensated for by continuously loading new photons into the device via coherent or incoherent input fields. This approach is very feasible as photons are much easier and cheaper to produce as compared to for example ultra-cold atoms [67, 68]. The emulated quantum many-body systems are therefore most naturally and feasibly explored in driven-dissipative regimes where input drives continuously replace the dissipated excitations and the dynamical balance of loading and loss processes eventually leads to stationary states.

This mode of operation should not be viewed as merely a means of compensating for an imperfection of the technology. In fact, quantum many-body systems are much less explored in such non-equilibrium scenarios than in equilibrium regimes. Investigating driven-dissipative regimes of interacting photons thus leads onto largely unexplored territory and may lead to interesting discoveries. One may even search for non-equilibrium phase transitions, that is points where the properties of the stationary state of some driven-dissipative dynamics change abruptly as one of the system's parameters is varied [137].

For investigations of driven-dissipative regimes, coherent driving fields at each lattice site are often considered. To be able to perform calculations with a time-independent Hamiltonian, it is often useful to move to a frame that rotates at the frequencies ${\omega }_{d}$ of the input fields. In this frame, the frequency of a photon in each cavity (resonator), ${\omega }_{r}$, is replaced by the detuning between ${\omega }_{r}$ and the frequency of the drive ${\omega }_{r}\to {{\rm{\Delta }}}_{r}={\omega }_{r}-{\omega }_{d}$ and similarly the transition frequency of emitters in the resonators is replaced by their detuning from the drive frequency, ${\omega }_{e}\to {{\rm{\Delta }}}_{e}={\omega }_{e}-{\omega }_{d}$. A field that continuously drives the cavity modes is, in this rotating frame, described by an additional term in the Hamiltonian,

$$\begin{eqnarray}&&{H}_{d}=\displaystyle \sum _{j}\left(\displaystyle \frac{{{\rm{\Omega }}}_{j}}{2}{e}^{-i{\varphi }_{j}}{a}_{j}^{\dagger }+{\rm{H.c.}}\right),\end{eqnarray} \tag{ 29 }$$

where ${{\rm{\Omega }}}_{j}$ is the amplitude and ${\varphi }_{j}$ the phase of the driving field at lattice site j. Note that while a global phase of all diving fields can be gauged away, the assumption of coherent drives requires choosing a relative phase between each pair of drives. For the dissipation, local particle losses are typically assumed and modelled by Lindblad type damping terms [35]. Hence, the driven-dissipative models discussed here are described by master equations of the from

$$\begin{eqnarray}&&\dot{\rho }=-\displaystyle \frac{i}{{\hslash }}[H,\rho ]+\displaystyle \sum _{j}\displaystyle \frac{\gamma }{2}(2{a}_{j}\rho {a}_{j}^{\dagger }-{a}_{j}^{\dagger }{a}_{j}\rho -\rho {a}_{j}^{\dagger }{a}_{j}),\end{eqnarray} \tag{ 30 }$$

where γ is the rate at which photons are lost from the device and H the Hamiltonian of the considered model including driving terms such as in equation (29) and written in a frame that rotates at the frequencies of the driving fields. Note that the form of the damping terms is invariant under the transformation into this rotating frame. In our discussion of stationary regimes of driven-dissipative quantum many-body models, we first consider the Bose–Hubbard model.

The driven-dissipative regime of a single resonator that features a Kerr nonlinearity for the cavity mode and thus corresponds to a single lattice site of the Hamiltonian (3) has been studied by Drummond and Walls [69] long before quantum simulation applications have been considered. These initial studies aimed to explore the physics of nonlinear polarisability and where later extended to consider solid-state nanostructures for single photon sources [81].

For a Bose–Hubbard model describing tunnel-coupled resonators, that is driven by coherent input fields at all lattice sites, one would expect that the field in the cavity array should inherit the classical and coherence properties of the driving fields provided their intensity is strong enough to dominate over the effects caused by the nonlinearities. The boundaries of this semi-classical regime have been calculated by Hartmann [107] via a linearisation in the quantum fluctuations around the classical background component, see figure 11.

Zoom In Zoom Out Reset image size

In the regime of strong interactions or nonlinearities, the number of excitations per lattice site is at most 1/2. This bound becomes obvious in the limit of very strong interactions, where each lattice site can be approximated by a two-level system as higher excited states are far off resonance to the drive. Yet, as a coherent field cannot generate inversion [253] these two-level systems have an excitation probability below 1/2. Consequently, Mott insulating regimes with commensurate filling can in this way not be generated.

A lattice system with low particle density, more precisely a system where the average inter particle spacing greatly exceeds the lattice constant, can be viewed as an approximation to a continuum system. The properties of a one-dimensional driven-dissipative Bose–Hubbard model in the strongly interacting regime should thus rather be compared to a Lieb-Liniger model, c.f. equation (16).

To explore whether Lieb-Liniger physics can be observed in such driven-dissipative regimes, Carusotto et al [40] considered a five-site version of the Hamiltonian (3) and investigated whether the characteristic energies of the collective strongly correlated many-body states could be seen in a spectroscopy analysis. They scanned the frequency of the driving fields through the relevant range and found resonance peaks at the transition-frequencies of the Hamiltonian (3).

Given the expected relations to the Lieb-Liniger model, an interesting question is, whether a driven-dissipative Bose–Hubbard model exhibits similar density-density correlations, in particular whether a feature similar to Friedel oscillations [87] can be expected. Spatially resolved density-density correlations in the form of a ${g}^{(2)}$-function,

$$\begin{eqnarray}&&{g}_{j,l}^{(2)}=\displaystyle \frac{\langle {a}_{j}^{\dagger }{a}_{l}^{\dagger }{a}_{l}{a}_{j}\rangle }{\langle {a}_{j}^{\dagger }{a}_{j}\rangle \langle {a}_{l}^{\dagger }{a}_{l}\rangle }\end{eqnarray} \tag{ 31 }$$

where j and l label lattice sites, have been investigated by Hartmann [107] via a numerical integration of equation (30) with Matrix Product Operators [111, 232]. A spatical modulation similar to Friedel oscillations was indeed found for a non-vanishing relative phase between the driving fields of adjacent lattice sites, see figure 12. As this phase off-set generates a particle flux in the array, the strong anti-bunching on-site ${g}_{j,j}^{(2)}\ll 1$, slight bunching for neighbouring sites, ${g}_{j,j+1}^{(2)}\gt 1$, and anti bunching for further separated sites, ${g}_{j,l}^{(2)}\lt 1$ for $| j-l| \gt 1$, indicates that a flux of particle dimers extended over neighbouring sites flows through the lattice. A similar signature was later also found for the driven-dissipative Jaynes–Cummings-Hubbard model by Grujic et al [95].

Zoom In Zoom Out Reset image size

For a higher dimensional lattice, a mean-field approach based on the exact single site solution by Drummond et al [69] was used by Le Boité et al [153, 154] to predict mono — and bistable phases, which emerge as a consequence of the nonlinear nature of the mean-field equations, see figure 13. The related dynamical hysteresis was then investigated in [43]. Yet other calculation techniques rather show a first order phase transition in the regions where mean-field predicts a bistable phase [174, 256].

Zoom In Zoom Out Reset image size

In the limit of very large on-site interactions, the driven-dissipative Bose–Hubbard model maps to a driven-dissipative XY spin-1/2 model, the phase diagram of which has been investigated with a site-decoupled mean-field approximation [258]. Here, stationary state phases with canted antiferromagnetic order and limit cycle phases with persistent oscillatory dynamics together with bistabilities of these two phases, have been found.

Moreover, driven-dissipative phases of the Bose–Hubbard model with cross-Kerr interactions as discussed in section 3.2.1 were calculated with a mean-field technique for higher [129] and with Matrix Product Operators for one dimension [128]. For sufficiently strong cross-Kerr interactions and low enough excitation tunnelling, the model exhibits a density-wave ordering with ${g}^{(2)}(j,j+1)\lt {g}^{(2)}(j,j)$, see figure 14.

Zoom In Zoom Out Reset image size

In the studies discussed so far, the excitation dissipation has been assumed to be caused by the electromagnetic vacuum throughout. As has been investigated in [203], the scenario changes substantially if squeezed dissipation is considered, which does not obey a U(1) symmetry. We now turn to discuss driven-dissipative scenarios for the Jaynes–Cummings-Hubbard model.

The coherence and fluorescence properties of a coherently pumped driven-dissipative Jaynes–Cummings-Hubbard model were explored by Nissen et al [189]. For short arrays, the photon blockade regime was found to persist even up to large tunnelling rates, whereas there is a transition to a coherent regime for larger arrays as the tunnelling strength is increased. This size dependence is due to the fact that spectrally dense excitation bands only form in the limit of large system sizes whereas for a small lattice, say a dimer, the tunnelling causes a splitting of the spectral lines of single site spectra that leads to a collective spectrum which still remains anharmonic [189], see figure 15.

Zoom In Zoom Out Reset image size

A comparative study of the features of driven-dissipative Bose–Hubbard and Jaynes–Cummings-Hubbard models was conducted by Grujic et al [95] and found quantitative differences for the experimentally accessible observables of both models for realistic regimes of interactions even when the corresponding nonlinearities are of similar strength.

Further interesting effects appear for arrays where not every resonator couples to a superconducting qubit but nonlinear interactions only occur in regularly spaced lattice sites. This periodic arrangement leads to photonic flat bands, see also section 5, where the role of interactions is enhanced and polaritons can become incompressible [25, 42]. Experimentally, similar flat bands have been generated in coupled micro-pillars, each containing a cavity formed by two distributed Bragg reflectors [13], where condensation of polaritons in the flat band was observed. Such micro-pillars have also been coupled in a honey-comb lattice where they exhibit edge states similar to graphene structures [195], c.f. section 5. Whereas we have so far discussed scenarios with coherent driving fields, cavity arrays with incoherent input fields have been considered as well.

In the above discussed driven-dissipative systems of interacting photons, the driving fields were always assumed to be of a coherent nature. Therefore any coherence between lattices sites that is found in the system needs to be attributed at least in part to the coherent inputs. To explore whether coherence can spontaneously develop for non-equilibrium photonic systems, incoherent input fields or pumping of higher excited states and subsequent decay [179] should be considered. A study for an incoherently pumped Jaynes–Cummings-Hubbard model found that the scaling of the correlation length with the hopping rate J_JC changes as the hopping rate becomes larger than the light-matter coupling g [219]. See also the remarks in section 4.1.1 about engineering effective chemical potentials. For further work on non-equilibrium photon condensation, see [3, 144, 236], as well as the reviews [39, 235] and references therein.

After reviewing approaches that considered a uniform drive intensity for all lattice sites, we now turn to discuss scenarios where the input can be more intense at one end of a chain. This setup naturally leads to the analysis of transport properties.

Transport of photons in a waveguide that are scattered at a localised emitter has been theoretically explored via scattering theory [165, 187, 234], see also extensions to multiple emitters [156], and numerically using the Density Matrix Renormalization Group (DMRG) [170]. Experimentally such scattering effects have been explored with Rydberg atoms [46, 115, 171, 218, 221] and one [12, 118] or two [247] superconducting qubits coupled to an open coplanar waveguide resonator. The recent developments in this direction of research are summarised in the review by Roy et al [218]. In our discussion in relation to quantum simulation we here therefore focus on transport studies of the driven-dissipative regime of a quantum many-body system on a one-dimensional chain or two-dimensional band. Here, a particle flux can either be generated by pumping locally at one end of the chain (the band) or by implementing a phase off-set between the driving fields at adjacent resonator in the direction of transport [107].

A phase transition in the sense that two eigenvalues of the Liouvillian approach zero has been found for a spin chain with incoherent pump at one and losses at the opposite end by Prosen and Pižorn [200]. Hafezi et al studied the propagation of few photon pulses in the polaritonic Lieb-Lininger model described in section 2.2 by decomposing their wave-function in zero-, singe- and two-photon components [100, 102] and found that for an input that drives a single- or two-photon transition, the output will contain anti-bunched or bunched photons. A different scenario with a coherent input at one end of the chain and uniform dissipation in all lattice sites has been studied by Biella et al [24], where the transport was found to be strongly influenced by the many-body resonances related to extended eigen-states of the chain. These transport properties can be interpreted as a generalisation of photon blockade, c.f. section 2.1.3, to extended one-dimensional systems. More recently, Mertz et al [182] explored two scenarios, a source-drain setup with coherent drive at one end and enhanced dissipation at the opposite end of the chain, as well as the regime where relative phases between coherent inputs at neighbouring lattice sites generate a current in the presence of uniform dissipation. Employing a Gutzwiller mean-field approximation, the study considered two dimensional lattices with periodic boundary conditions in the direction perpendicular to the direction of transport. In addition to the dependencies of the current on the many-body spectrum it found that transport can be inhibited by strong dissipation at the 'drain'-end by the quantum Zeno effect.

A relative phase between the coherent input fields at different locations can already lead to interesting effects when considered for the two outer resonators of a three-site model, which we turn to discuss now.

The interplay of coherent tunnelling and on-site repulsion has been theoretically explored in three-cavity setups where the two outer cavities where driven by coherent fields the relative phase of which was varied. The central cavity in turn contained a Kerr nonlinearity, c.f. equation (3). The first study coined the setup 'Quantum Optical Josephson Interferometer' [90] and found that the destructive interference in the central cavity due to opposite phases of the driving fields is retained even for large nonlinearity. As photons enter the central cavity via tunnelling processes from the outer cavities, the transition from coherent to anti-bunched photon statistics in this cavity depends on the tunnelling rate. Similar behaviour was found for a case where two waveguides with a continuous spectrum replaced the two outer cavities. A later work explored the setup in a superconducting circuit context [130], c.f. section 3, and extended the study to higher excitation numbers where a different dependence of the coherent to anti-bunching transition on the tunnelling rate was found.

Due to the large dimension of their Hilbert spaces and the eventually long time scales on which stationary states are reached, modelling many-body systems of interacting photons is a demanding challenge. In the next section we thus review some existing powerful methods together with recent efforts to extend their application ranges and develop new ones.